Using velocity profiles with DiffEqOperators.jl for 2D convection-diffusion problem

Specific Domains Modelling & Simulations
1

Hello all. I am trying to discretize and solve the following PDE using DiffEqOperators.jl

\frac{\partial C}{\partial t} = D\nabla^2 C - \vec{v} \cdot \nabla C + R

for solving non-steady convection diffusion of a generic chemical species with a generic velocity profile v_x(y).

So far I have been able to generate and solve the problem using scalar values for the velocity field in each direction in both 2D and 3D using the following code:


# Loading Dependencies

using DifferentialEquations

using DiffEqOperators

using Plots

# Physical Constants

α = 0.05 # Diffusion Coefficient

v_x = 0.5; # Advection speed in x [m/s]

v_y = 0.25; # Advection speed in y [m/s]

γ=0.05      # Source Term

# Discretization Parameters

s = x, y = (-5:0.2:5, -5:0.2:5)

dx = dy = x[2] - x[1]

t0 = 0.0

t1 = 20.0

# Initial Distribution

u_low = convert(Int64, round(2 / 10 * length(x)))

u_high = convert(Int64, round(4 / 10 * length(x)))

u0 = fill(0.0, (length(x), length(y)));

u0[u_low:u_high,u_low:u_high] .= 5.0;   # a small square in the middle of the domain

# Generate Coefficient Matrices

laplacianx = CenteredDifference{1}(2, 4, dx, length(x))

laplaciany = CenteredDifference{2}(2, 4, dy, length(y))

Δ = laplacianx + laplaciany

divx = CenteredDifference{1}(1, 4, dx, length(x))

divy = CenteredDifference{2}(1, 4, dy, length(y))

# Generate Boundary Conditions

Qx, Qy = Neumann0BC(Float64, (dx, dy), 4, (length(x), length(y)))

Q = compose(Qx, Qy)

# Define matrix system and solve

function f!(du, u, p, t)

    du .= α .* (Δ * Q * u) + (((-v_x * divx) + (-v_y * divy)) * Q * u) + γ .* u

###        #Diffusion         #Advection in x   #Advection in y           #Source

end

prob = ODEProblem(f!, u0, (t0, t1))

alg = Tsit5()

sol2D = solve(prob, alg)

My question is whether or not it is possible to use a velocity profile for the v_x, v_y, v_w terms and some recommendations on how to go about it. Likewise I was wondering if there is a more compact form of the \vec{v} \cdot \nabla C expression that can be expressed using DiffEqOperators.jl. Thanks!

2

I don’t quite get the question. Try it again?

3

Alright so the \vec{v} \cdot \nabla C portion of my governing equation models advection right. What this looks like expanded for an incompressible flow is v_x \frac{\partial C}{\partial x} \vec{i} + v_y \frac{\partial C}{\partial y} \vec{j} + v_z \frac{\partial C}{\partial z} \vec{k}. Now when I say velocity profile, I mean using a function to give a varying velocity speed dependent on the position. So because they are not allowed to vary along its prinicipal axis, this would look like v_x = v_x(y,z), v_y = v_y(x,z), v_z(y,z) etc.
Now for the application, I was able to model a fixed scalar advection speed using the following code:

function f!(du, u, p, t)
    du .= α .* (Δ * Q * u) + (((-v_x * divx) + (-v_y * divy)) * Q * u) + γ .* u
###        #Diffusion         #Advection in x   #Advection in y           #Source
end

Where v_x, v_y are constant scalars and then divx is the divergence in x and so on. So what I would like to do is instead use a parabolic flow field given by v_x(y) = v_{max} \left[ 1 - \left( \frac{y - H/2}{H/2} \right) \right], for example but this could be any arbitrary function for the flow. So I’m thinking that if I can generate v_x to be an array of the appropriate shape then I could do .* with the divergence operator divx for each direction and this will give me the behavior I want. I just don’t know what the ordering/shape would be like so that it corresponds to each line of y. Then lastly, I was wondering if there were any functions I could use to just do the dot multiply operation of \vec{v} \cdot \nabla C instead of having to expand it like this (((-v_x * divx) + (-v_y * divy)) * Q * u) (this one is more just wishful thinking). Hopefully this clears things up.

TLDR cross or dot product of a vector function by a DiffEqOperator.jl derivative operator (in this case divx = CenteredDifference{1}(1, 4, dx, length(x)) in order to perform operations such as \vec{v} \cdot \nabla C which come up often in PDEs.

Thanks again.

4

Make v_x a vector and use -v_x .* divx

5

Hi Chris, I tried your suggestion and I got an error about a missing method regarding the .* when doing v_x .* divx. However I got a solution for what I’m trying to do that I want to document for other people trying to do dot products in DiffEqOperators.jl. The trick was actually changing the order of the operations for the governing equation. Essentially we went from


function f!(du, u, p, t)

    du .= α .* (Δ * Q * u) + (((-v_x .* divx) + (-v_y .* divy)) * Q * u) + γ .* u

end

to this


function f!(du, u, p, t)

    du .= α .* (Δ * Q * u) + ((-v_x' .* (divx * Q * u) + -v_y' .* (divy * Q * u)) ) + γ .* u

end

Following this change, v_x and v_y and even \alpha can all be matricies and I can functionally perform the dot product \vec{v} \cdot \nabla C since I can multiply each part of the gradient by the correct velocity component a la
v_x \frac{\partial C}{\partial x} \vec{i} + v_y \frac{\partial C}{\partial y} \vec{j} + v_z \frac{\partial C}{\partial z} \vec{k}.

So now then, the reason all of this was important to me was so that now I can have arbitrary flow fields (or more generally nonconstant properties that vary along the domain) being modeled which opens up the possibilities of the work I can do with this package. Thanks for the help.

1 Like