MAE 103: Differential Equations in Julia

Solve differential equations in Julia and pathlines

In today's lecture, we will learn how to solve ordinary differential equations. We will then use it to plot pathlines of fluid particles

A 1D differential equation

In this first example, we will solve the equation: dudt=−u \frac{d u}{dt} = -u with u(t=0)=1u(t=0) = 1

We know the analytical solution for this problem : u(t)=exp(−t)u(t) = exp(-t)

In Julia, we can call the package OrdinaryDiffEq to solve numerically ODEs

using DifferentialEquations
Shift+Enter to run
Julia
Julia 1.2

The first step is to define the right hand side of the differential equation, here f(u,p,t) = -u

f(u,p,t) = -u
Shift+Enter to run
Julia
Julia 1.2
nil

The initial condition is u0 = 1

u0=1.0
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)
Shift+Enter to run
Julia
Julia 1.2
nil

The variable prob holds all the useful informations about the differential problem : right-hand side, initial condition, time span...

We will solve this ODe with a Runge-Kutta 4 scheme

sol = solve(prob,RK4())
Shift+Enter to run
Julia
Julia 1.2
nil

Plotting the solution

using Plots
Shift+Enter to run
Julia
Julia 1.2
plot(sol.t, exp.(-sol.t), linewidth = 4, label = "Analytical solution")
plot!(sol,linewidth=4, xaxis="t",yaxis="u", linestyle =:dash, label = "Numerical solution" ) # legend=false
Shift+Enter to run
Julia
Julia 1.2
image/svg+xml
using Interpolations
Shift+Enter to run
Julia
Julia 1.2

Solve for the pathline of a fluid particle in a Couette Flow

For a Couette flow, the velocity profile is: u=u(y)ex=UyHex \boldsymbol{u} = u(y)\boldsymbol{e}_x = U\frac{y}{H} \boldsymbol{e}_x with UU the velocity of the upper plate and HH the depth of the gap

The trajectory of a fluid particle A is denoted xA\boldsymbol{x}_A and a solution of:

dxAdt=uA(t)=u(xA,yA,zA,t)\frac{d \boldsymbol{x}_A}{dt} = \boldsymbol{u}_A(t) = \boldsymbol{u}(x_A, y_A, z_A, t)

with initial condition xA(t=0)=x0\boldsymbol{x}_A(t= 0) = \boldsymbol{x}_0

Let's solve for the pathline of a fluid particle. For this problem we have :

dxAdt=u(xA,yA,t)=u(yA)=UyAH \frac{d x_A}{dt} = u(x_A, y_A,t) = u(y_A) = U\frac{y_A}{H}
dyAdt=v(xA,yA,t)=0 \frac{d y_A}{dt} = v(x_A, y_A,t) = 0

Define the parameters for the Couette flow

U = 1
H = 1
couette_velocity(y) = U*y/H
Shift+Enter to run
Julia
Julia 1.2
nil

Let's define the right-hand side of the equation

function couette(du,u,p,t)
 du[1] = couette_velocity(u[2]) # U * (yₐ/H)
 du[2] = 0.0
end
Shift+Enter to run
Julia
Julia 1.2
nil
## Initial position of the fluid particle
x0=[0.0, 0.5]
# Time span
tspan = (0.0,1.0)
prob_couette = ODEProblem(couette, x0, tspan)
Shift+Enter to run
Julia
Julia 1.2
nil

Now, we can solve the differential equation and get the pathline of this fluid particle

couette_sol  = solve(prob_couette, RK4())
Shift+Enter to run
Julia
Julia 1.2
nil
plot(couette_sol, vars = 1, legend = false, xlabel = "t", ylabel = "xA(t)")
Shift+Enter to run
Julia
Julia 1.2
image/svg+xml
plot(couette_sol, vars = 2, legend = false, xlabel = "t", ylabel = "yA(t)")
Shift+Enter to run
Julia
Julia 1.2
image/svg+xml
plot(couette_sol, vars = (1,2), ylim =(0, H))
plot!([1.0], line = :hline, linewidth = 8, label = "Upper wall")
plot!([0.0], line = :hline, linewidth = 8, label = "Lower wall")
scatter!([x0[1]], [x0[2]])
Shift+Enter to run
Julia
Julia 1.2
image/svg+xml

Trajectory of a fluid particle in a cellular flow

A = 32.0
B = 3.2
ω = 9.2
ψ(x, y, t) = A*cos(x + B*sin(ω*t))*cos(y)
ucell(x, y, t) =  A*cos(x + B*sin(ω*t))*sin(y)
vcell(x, y, t) = -A*sin(x + B*sin(ω*t))*cos(y)
Shift+Enter to run
Julia
Julia 1.2
nil
xrange = -4.0:0.05:12.0
yrange = -4.0:0.05:12.0
function ψtab(t) 
    out  = zeros(xrange.len ,yrange.len)
for (i, xi) in enumerate(xrange)
    for (j, yj) in enumerate(yrange)
    
        out[i, j] = ψ(xi, yj , t)
        
    end
end
    
    return out
end
Shift+Enter to run
Julia
Julia 1.2
nil
plot(xrange, yrange, ψtab(0.0)', ratio = 1, xlim = (xrange[1], xrange[end]), ylim = (yrange[1], yrange[end]), colorbar = false,
            xlabel ="x", ylabel = "y", title = "Streamlines at t = "*string(0.0))
Shift+Enter to run
Julia
Julia 1.2
image/svg+xml
plot(xrange, yrange, ψtab(1.0)', ratio = 1, xlim = (xrange[1], xrange[end]), ylim = (yrange[1], yrange[end]), colorbar = false,
            xlabel ="x", ylabel = "y", title = "Streamlines at t = "*string(1.0))
Shift+Enter to run
Julia
Julia 1.2
image/svg+xml
plot(xrange, yrange, ψtab(2.0)', ratio = 1, xlim = (xrange[1], xrange[end]), ylim = (yrange[1], yrange[end]), colorbar = false,
            xlabel ="x", ylabel = "y", title = "Streamlines at t = "*string(2.0))
Shift+Enter to run
Julia
Julia 1.2
image/svg+xml
function cellular_flow(du,u,p,t)
 du[1] = ucell(u[1], u[2], t)
 du[2] = vcell(u[1], u[2], t)
end
Shift+Enter to run
Julia
Julia 1.2
nil
## Initial position of the fluid particle
x0_cell=[3.20, 1.88]
# Time span
tspan = (0.0,2.0)
prob_cellular = ODEProblem(cellular_flow, x0_cell, tspan)
Shift+Enter to run
Julia
Julia 1.2
nil
cellular_sol  = solve(prob_cellular, RK4())
Shift+Enter to run
Julia
Julia 1.2
nil
plot(xrange, yrange, ψtab(0.0)', ratio = 1, xlim = (xrange[1], xrange[end]), ylim = (0.0, 5.0), colorbar = false,
            xlabel ="x", ylabel = "y", title = "Streamlines at t = "*string(0.0))
plot!(cellular_sol, vars = (1,2))
scatter!([x0_cell[1]], [x0_cell[2]])
Shift+Enter to run
Julia
Julia 1.2
image/svg+xml
Shift+Enter to run
Julia
Julia 1.2
Shift+Enter to run
Julia
Julia 1.2
Runtimes (1)