103: 1D Transient convection-diffusion equation

(source code)

Solve the time-dependent convection-diffusion equation

\[\partial_t u -\nabla ( D \nabla u - v u) = 0\]

in $\Omega=(0,1)$ with homogeneous Neumann boundary condition at $x=0$ and outflow boundary condition at $x=1$.

This is the time-dependent version of the convection-diffusion problem from Example102. The equation models the evolution of a scalar quantity $u$ (e.g., concentration, temperature) under the combined effects of diffusion (with coefficient $D$) and advection (with velocity $v$).

Physical Interpretation

• Diffusion term: $D \nabla u$ represents spreading due to random motion
• Convection term: $v u$ represents transport by bulk motion of the medium
• Time evolution: $\partial_t u$ describes how the quantity changes over time

Boundary Conditions

• Left boundary ($x=0$): Homogeneous Neumann condition $\partial_n u = 0$ (no flux)
• Right boundary ($x=1$): Outflow condition $D\partial_n u = 0$, equivalent to $ ( D \nabla u - v u)\cdot \vec n = v u \cdot \vec n $

allowing material transported by convection towards the boundary to leave the domain.

Discretization

The spatial discretization uses the exponential fitting scheme (Scharfetter-Gummel method) which maintains monotonicity properties and provides accurate solutions even for convection-dominated transport (high Peclet numbers). This method uses the Bernoulli function:

\[B(x) = \frac{x}{e^x - 1}\]

to construct fluxes that exactly solve the local two-point boundary value problem on each edge.

Initial Condition

The simulation starts with a linear profile $u(x,0) = 1 - 2x$, which evolves under the combined effects of diffusion and convection until it reaches a steady state or is advected out of the domain.

module Example103_ConvectionDiffusion1Dusing Printfusing VoronoiFVMusing ExtendableGridsusing GridVisualize# Mutable struct to hold problem parameters# This encapsulates all physical and numerical parametersmutable struct ProblemData    D::Float64      ## Diffusion coefficient    v::Vector{Float64}  ## Velocity vectorend# Bernoulli function used in the exponential fitting discretizationfunction bernoulli(x)    if abs(x) < nextfloat(eps(typeof(x)))        return 1    end    return x / (exp(x) - 1)endfunction exponential_flux!(f, u, edge, data)    vh = project(edge, data.v)    Bplus = data.D * bernoulli(vh / data.D)    Bminus = data.D * bernoulli(-vh / data.D)    f[1] = Bminus * u[1, 1] - Bplus * u[1, 2]    return nothingendfunction outflow!(f, u, node, data)    if node.region == 2        f[1] = data.v[1] * u[1]    end    return nothingendfunction main(; n = 10, Plotter = nothing, D = 0.01, v = 1.0, tend = 100)    # Create a one-dimensional discretization    h = 1.0 / n    grid = simplexgrid(0:h:1)    # Initialize problem parameters in data structure    problem_data = ProblemData(D, [v])    sys = VoronoiFVM.System(        grid,        VoronoiFVM.Physics(;            flux = exponential_flux!,            breaction = outflow!,            data = problem_data        )    )    # Add species 1 to region 1    enable_species!(sys, 1, [1])    # Set boundary conditions    boundary_neumann!(sys, 1, 1, 0.0)    # Create a solution array    inival = unknowns(sys)    inival[1, :] .= map(x -> 1 - 2x, grid)    # Transient solution of the problem    control = VoronoiFVM.SolverControl()    control.Δt = 0.01 * h    control.Δt_min = 0.01 * h    control.Δt_max = 0.1 * tend    tsol = solve(sys; inival, times = [0, tend], control)    vis = GridVisualizer(; Plotter = Plotter)    for i in 1:length(tsol.t)        scalarplot!(            vis[1, 1], grid, tsol[1, :, i]; flimits = (0, 1),            title = "t=$(tsol.t[i])", show = true        )        sleep(0.01)    end    return tsolendusing Testfunction runtests()    tsol = main()    @test maximum(tsol) <= 1.0 && maximum(tsol.u[end]) < 1.0e-20    return nothingendend

This page was generated using Literate.jl.