Solving differential equations for multi-component fluid flow

Specific Domains Modelling & Simulations
, ,
1

Hello everyone,

I would like to model a multi-component, isothermal gas flow in a 1D pipeline. To do this, I would like to solve the following simplified differential equations:

\quad \frac{\partial p}{\partial t} = - a \frac{\partial q_{m}}{\partial x} \\ \quad \frac{\partial p}{\partial x} = -b \frac{q_{m}|q_{m}|}{p} \\ \quad \frac{\partial c}{\partial t} = - a \frac{q_m}{p}\frac{\partial c}{\partial x}

The variables are the pressure p(x,t), the mass flow q_m(x,t) and the methane mass fraction c(x,t). a and b are constants. The initial state is the steady-state flow, i.e.

\quad q_m(x,0) = q_0 \\ \quad p(x,0) = \sqrt{p_0^2 - 2b\, q_0|q_0|x} \\ \quad c(x,0) = 1 \,.

For the boundary conditions I choose the outlet mass flow, the inlet pressure and the methane fraction,

\quad p(0,t) = p_{in}(t) \ \quad c(0,t) = c_{in}(t) \ \quad q_m(L,t) = q_{out}(t)\, .
In a first step, these boundary conditions can be chosen as constant. Later I would like to use time-dependent input.

How do I solve this problem with the Julia (DifferentialEquations)? In particular, what do I need to avoid numerical diffusion and get a numerically stable system?

Input

Here are some realistic input values:

L = 10e3
x = [0,L]
t = [0,43200]  
p_0 = 50e5
q_0 = 18 
a = 1.3e6
b = 1.9e5
2

Not sure if this is your complete problem, or just a small starting point. I would suggest to consider using Home · ModelingToolkit.jl , because it allows you just to write down your differential equations in symbolic form and it simplifies your system of differential equations symbolically before solving it numerically.

If this is is your complete problem, using ModelingToolkit might be overkill, though.

1 Like
3 

using MethodOfLines, ModelingToolkit, DomainSets

@parameters t, x, a, b
@variables p(..), q_m(..), c(..)  

Dt = Differential(t); Dx = Differential(x)

a, b = 1.3e6, 1.9e5 
p0, q0 = 5e6, 18

eq = [
    Dt(p(t, x)) ~ -a * Dx(q_m(t, x)), 
    Dx(p(t, x)) ~ -b * q_m(t, x) * abs(q_m(t, x)) / p(t, x),
    Dt(c(t, x)) ~ -a * q_m(t, x) / p(t, x) * Dx(c(t, x)), 
]

L = 1000.0

domain = [x ∈ Interval(0.0, L),
          t ∈ Interval(0.0, 43200.0)
]


ic_bc = [q_m(0.0, x) ~ q0,
         p(0, x) ~ sqrt(p0^2 - 2b * q0 * abs(q0) * x),
         c(0, x) ~ 1,
         q_m(t, L) ~ 1, # boundary condition be chosen as constant
         p(t, 0) ~ p0,
         c(t, 0) ~ 1    # boundary condition be chosen as constant
]

@named sys = PDESystem(eq, ic_bc, domain, [t, x], [p(t, x), q_m(t, x), c(t, x)])

dx = 2.0

discretization = MOLFiniteDifference([x => dx], t, approx_order = 2)

prob = discretize(sys, discretization)

using DifferentialEquations

sol = solve(prob, saveat = 1.0)

which throws error:

ERROR: MethodError: no method matching SciMLBase.PDETimeSeriesSolution(::SciMLBase.PDETimeSeriesSolution{Float64, 3, Dict{Num, Matrix{Float64}}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, Nothing, MethodOfLines.ScalarizedDiscretization}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, ModelingToolkit.MTKParameters{Tuple{}, Tuple{}, Tuple{}, Tuple{}, Tuple{}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out,), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x856c221d, 0x2ab738c3, 0x90c8aeff, 0x564502a8, 0x9bbfff35), Expr}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xc9c2ab17, 0x3395fe14, 0x86c47b71, 0xce7f7866, 0xb9d3bbeb), Expr}}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#740"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x244bcc7c, 0x4ad8c85d, 0xd984c452, 0xe91bc568, 0x17bbcbce), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x29d5008a, 0x0fae669a, 0xc0886d32, 0xb3a1ef51, 0x6483659f), Nothing}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, Nothing, MethodOfLines.ScalarizedDiscretization}}, CompositeAlgorithm{1, Tuple{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Rosenbrock23{1, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, Rodas5P{1, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, FBDF{5, 1, false, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, typeof(OrdinaryDiffEq.trivial_limiter!)}, FBDF{5, 1, false, KrylovJL{typeof(Krylov.gmres!), Int64, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, typeof(OrdinaryDiffEq.trivial_limiter!)}}, OrdinaryDiffEq.AutoSwitchCache{Tuple{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Vern7{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, Tuple{Rosenbrock23{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, Rodas5P{0, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, FBDF{5, 0, false, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, typeof(OrdinaryDiffEq.trivial_limiter!)}, FBDF{5, 0, false, KrylovJL{typeof(Krylov.gmres!), Int64, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}}, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, typeof(OrdinaryDiffEq.trivial_limiter!)}}, Rational{Int64}, Int64}}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#740"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x244bcc7c, 0x4ad8c85d, 0xd984c452, 0xe91bc568, 0x17bbcbce), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x29d5008a, 0x0fae669a, 0xc0886d32, 0xb3a1ef51, 0x6483659f), Nothing}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, Vector{Int64}, OrdinaryDiffEq.DefaultCache{OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, uNoUnitsType, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False} where uNoUnitsType<:(Vector), OrdinaryDiffEq.Vern7Cache{Vector{Float64}, Vector{Float64}, uNoUnitsType, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False} where uNoUnitsType<:(Vector), OrdinaryDiffEq.Rosenbrock23Cache{Vector{Float64}, Vector{Float64}, uNoUnitsType, Matrix{Float64}, Matrix{Float64}, TabType, SciMLBase.TimeGradientWrapper{true, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#740"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x244bcc7c, 0x4ad8c85d, 0xd984c452, 0xe91bc568, 0x17bbcbce), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x29d5008a, 0x0fae669a, 0xc0886d32, 0xb3a1ef51, 0x6483659f), Nothing}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, Vector{Float64}, ModelingToolkit.MTKParameters{Tuple{}, Tuple{}, Tuple{}, Tuple{}, Tuple{}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out,), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x856c221d, 0x2ab738c3, 0x90c8aeff, 0x564502a8, 0x9bbfff35), Expr}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xc9c2ab17, 0x3395fe14, 0x86c47b71, 0xce7f7866, 0xb9d3bbeb), Expr}}}, SciMLBase.UJacobianWrapper{true, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#740"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x244bcc7c, 0x4ad8c85d, 0xd984c452, 0xe91bc568, 0x17bbcbce), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x29d5008a, 0x0fae669a, 0xc0886d32, 0xb3a1ef51, 0x6483659f), Nothing}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, Float64, ModelingToolkit.MTKParameters{Tuple{}, Tuple{}, Tuple{}, Tuple{}, Tuple{}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out,), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x856c221d, 0x2ab738c3, 0x90c8aeff, 0x564502a8, 0x9bbfff35), Expr}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xc9c2ab17, 0x3395fe14, 0x86c47b71, 0xce7f7866, 0xb9d3bbeb), Expr}}}, F, FiniteDiff.JacobianCache{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, UnitRange{Int64}, Nothing, Val{:forward}(), Float64}, FiniteDiff.GradientCache{Nothing, Vector{Float64}, Vector{Float64}, Float64, Val{:forward}(), Float64, Val{true}()}, Float64, Rosenbrock23{1, false, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)}, Vector{Bool}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!)} where {uNoUnitsType<:(Vector), TabType<:OrdinaryDiffEq.Rosenbrock23Tableau, F<:(LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, LinearSolve.DefaultLinearSolver})}, OrdinaryDiffEq.Rosenbrock5Cache{Vector{Float64}, Vector{Float64}, uNoUnitsType, Matrix{Float64}, Matrix{Float64}, TabType, SciMLBase.TimeGradientWrapper{true, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#740"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x244bcc7c, 0x4ad8c85d, 0xd984c452, 0xe91bc568, 0x17bbcbce), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x29d5008a, 0x0fae669a, 0xc0886d32, 0xb3a1ef51, 0x6483659f), Nothing}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ODESystem}, Nothing, ODESystem, Nothing, Nothing}, Vector{Float64}, ModelingTool…
Stacktrace:

wrap_sol(sol::SciMLBase.PDETimeSeriesSolution{…}, metadata::MethodOfLines.MOLMetadata{…}) at pde_solutions.jl

wrap_sol(sol::SciMLBase.PDETimeSeriesSolution{…}) at basic_solutions.jl

solve(::ODEProblem{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::Base.Pairs{…}) at solve.jl

kwcall(::@NamedTuple{…}, ::typeof(solve), ::ODEProblem{…}) at solve.jl

top-level scope at multi-component fluid flow.jl
4

This is just a starting point for a larger problem. I’ve tried using the ModelingToolkit, but I haven’t managed to get a stable solution and the propagation stops early. The mass flow of the last time step shows some unphysical oscillations, it seems to me that the boundary conditions might be the problem, but I have not managed to fix it.

Code

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets, Plots

function test2(L::Float64, T::Float64)
    a = 1.277e6
    b = 1.865e5
    q = 15.0
    p_in = 5e6
 
    @parameters x t
    @variables p(..) q_m(..) c(..)
    ∂t = Differential(t)
    ∂x = Differential(x)
 
    domains = [
        x ∈ Interval(0, L),
        t ∈ Interval(0, T),
    ]
 
    eq = [
        ∂t(p(x, t)) ~ (-a) * ∂x(q_m(x, t)),
        ∂x(p(x, t)) ~ - b * q_m(x, t) * abs(q_m(x, t)) / p(x, t),
        ∂t(c(x, t)) ~ - a * q_m(x, t) / p(x, t) * ∂x(c(x, t)),
    ]
 
    p_0(x) = sqrt((p_in)^2 - 2 * b * q * q * x)

    bcs = [
        # -------------------
        # Boundary conditions
        # -------------------
        p(0, t)   ~ p_in, 
        q_m(L, t) ~ q, 
        c(0, t)   ~ 1.0,
        # -------------------
        # Initial conditions
        # -------------------
        p(x, 0)   ~ p_0(x),
        q_m(x, 0) ~ q_m0(x),
        c(x, 0)   ~ 1.0,
    ]
 
    @named pde_system = PDESystem(
        eq, bcs, domains,
        [x, t],
        [p(x, t), q_m(x, t), c(x, t)],
    )
 
    return pde_system
end
 
pde_system = test2(500.0, 3600.0 * 1)
x, t = pde_system.ivs
p, q_m, c = pde_system.dvs
 
discretization = MOLFiniteDifference([x => 50], t; approx_order=2)
 
prob = discretize(pde_system, discretization)
sol = solve(prob, Rosenbrock23())

Output/Error

PDESystem
Equations: Equation[Differential(t)(p(x, t)) ~ -1.277e6Differential(x)(q_m(x, t)), Differential(x)(p(x, t)) ~ (-186500.0abs(q_m(x, t))*q_m(x, t)) / p(x, t), Differential(t)(c(x, t)) ~ (-1.277e6q_m(x, t)*Differential(x)(c(x, t))) / p(x, t)]
Boundary Conditions: Equation[p(0, t) ~ 5.0e6, q_m(500.0, t) ~ 15.0, c(0, t) ~ 1.0, p(x, 0) ~ sqrt(2.5e13 - 8.3925e7x), q_m(x, 0) ~ 15.0, c(x, 0) ~ 1.0]
Domain: Symbolics.VarDomainPairing[Symbolics.VarDomainPairing(x, 0.0 .. 500.0), Symbolics.VarDomainPairing(t, 0.0 .. 3600.0)]
Dependent Variables: Num[p(x, t), q_m(x, t), c(x, t)]
Independent Variables: Num[x, t]
Parameters: SciMLBase.NullParameters()
Default Parameter ValuesDict{Any, Any}()

2-element Vector{Num}:
 x
 t

3-element Vector{Num}:
   p(x, t)
 q_m(x, t)
   c(x, t)

MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}(Dict{Num, Int64}(x => 50), t, 2, UpwindScheme(1), MethodOfLines.CenterAlignedGrid(), true, false, MethodOfLines.ScalarizedDiscretization(), true, Any[], Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}())

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 3600.0)
u0: 147-element Vector{Float64}:
  4.999914361511503e6
  4.999828721556165e6
  4.999743080133912e6
  4.999657437244668e6
  4.999571792888356e6
  4.999486147064904e6
  4.999400499774233e6
  4.999314851016271e6
  4.99922920079094e6
  4.999143549098165e6
  4.999057895937871e6
  4.998972241309983e6
  4.998886585214425e6
  ⋮
 15.0
 15.0
 15.0
 15.0
 15.0
 15.0
 15.0
 15.0
 15.0
 15.0
 15.0
 15.0
 15.0

┌ Warning: At t=0.00031538038626707164, dt was forced below floating point epsilon 5.421010862427522e-20, and step error estimate = 3.677287959149432. Aborting. There is either an error in your model specification or the true solution is unstable (or the true solution can not be represented in the precision of Float64).
└ @ SciMLBase C:\Users\OrsolitsD\.julia\packages\SciMLBase\sakPO\src\integrator_interface.jl:633
PDESolution:
  Return Code:
    Unstable
  Dependent variables:
    p(x, t): (50, 533) sized solution
    q_m(x, t): (50, 533) sized solution
    c(x, t): (50, 533) sized solution
  Domain:
    t ∈ (0.0, 0.00031538038626707164) with 533 non-uniform points. average step size 5.928202749381046e-7
    x ∈ (0.0, 500.0) with 50 points, step size 10.204081632653061
  From system:
    Equations:
L"\begin{align}
1.277 \cdot 10^{6} \frac{\mathrm{d}}{\mathrm{d}x} q_{m}\left( x, t \right) + \frac{\mathrm{d}}{\mathrm{d}t} p\left( x, t \right) =& 0 \\
\frac{1.865 \cdot 10^{5} q_{m}\left( x, t \right) \left|q_{m}\left( x, t \right)\right|}{p\left( x, t \right)} + \frac{\mathrm{d}}{\mathrm{d}x} p\left( x, t \right) =& 0 \\
\frac{1.277 \cdot 10^{6} q_{m}\left( x, t \right) \frac{\mathrm{d}}{\mathrm{d}x} c\left( x, t \right)}{p\left( x, t \right)} + \frac{\mathrm{d}}{\mathrm{d}t} c\left( x, t \right) =& 0
\end{align}
"
    Boundary/Initial Conditions:
L"\begin{align}
p\left( 0, t \right) =& 5 \cdot 10^{6} \\
q_{m}\left( 500, t \right) =& 15 \\
c\left( 0, t \right) =& 1 \\
p\left( x, 0 \right) =& \sqrt{2.5 \cdot 10^{13} - 8.3925 \cdot 10^{7} x} \\
q_{m}\left( x, 0 \right) =& 15 \\
c\left( x, 0 \right) =& 1
\end{align}

Plot Output

plot(sol[x][:], sol[q_m][:,end], label="mass flow")
5

Some remarks:
DifferentialEquations.jl supports supports hundreds of solvers, you might need one or two days to find the one that is best suited for your problem. In addition it might be important to set the absolute and relative accuracy of the solvers to a lower value, e.g. 1e-7.

I use (for a very different problem) code like this:

    prob = ODEProblem(m.sys_simple, u0, tspan, p)
    tol=1e-7
    sol = solve(prob, Rodas5(), dt=dt, tstops=tstops, abstol=tol, dtmax=1, reltol = tol, saveat=ts)
6

Rodas5P or FBDF would be a better choice. Rosenbrock23 isn’t great at DAEs.

2 Likes
7

It is possible that swapping the boundary conditions will lead to more stable results: The flow should be specified upstream and the pressure downstream.
A weno discretization of the space derivatives could also be helpful.

And
0 ~ ∂x(p(x, t)) + b * q_m(x, t) * abs(q_m(x, t)) / p(x, t),
instead of
∂x(p(x, t)) ~ - b * q_m(x, t) * abs(q_m(x, t)) / p(x, t),
seems preferable.

1 Like
8

I didn’t change the BC’s as they wouldn’t fit into the structure of my problem. However, rewriting the second equation and using FBDF()/Rodas5P() leads to a stable solution. I could not observe any difference for different advection schemes.