B-LIE
commented
1 year ago In case you are interested in my code [note! I have successfully solved at least 3 PDEs using MethodsOfLine.jl, so I have a little experience with it]:
# PACKAGES
using ModelingToolkit
using DifferentialEquations
using Plots, LaTeXStrings
using MethodOfLines, DomainSets
# MODEL
parameters RdM T A_p g H L ℘ f
@variables t x p(..) ṁ(..) # p_0(t) ṁ_L(t)
Dt = Differential(t)
Dx = Differential(x)
#
eqs = [Dt(p(t,x)) ~ -(RdM*T)/(A_p)*Dx(ṁ(t,x)),
Dt(ṁ(t,x)) ~ -(RdM*T)/(A_p)*Dx((ṁ(t,x))^2/p(t,x)) - A_p*Dx(p(t,x)) +
(g*H)/L*A_p/(RdM*T)*p(t,x) - RdM*T*℘/A_p^2/2*f*(ṁ(t,x))^2/p(t,x)]
#= In case term needs to be expanded out:
eqs = [Dt(p(t,x)) ~ -(RdM*T)/(A_p)*Dx(ṁ(t,x)),
Dt(ṁ(t,x)) ~ -(RdM*T)/(A_p)*( 2*ṁ(t,x)/p(t,x)*Dx(ṁ(t,x)) - (ṁ(t,x))^2/p(t,x)*Dx(p(t,x)) ) - A_p*Dx(p(t,x)) +
(g*H)/L*A_p/(RdM*T)*p(t,x) - RdM*T*℘/A_p^2/2*f*(ṁ(t,x))^2/p(t,x)]
=#
# NUMERIC VALUES
# From A. Carlson and B. Lie, 1992
D0 = 1.422 # m
Ap0 =pi*D0^2/4
℘0 = pi*D0
L0 = 122e3 # m
RdM0 = 493 # J/(kg.K)
T0 = 12 + 273.15 # K
f0 = 0.28e-3 # friction factor
g0 = 9.81
H0 = 0
p0 = 8.4e6 # Pa
ṁ0 = 1.5e6/3600 # kg/s
# BOUNDARY FUNCTIONS
@register_symbolic p_0(t)
@register_symbolic ṁ_L(t)
p_0(t) = p0
ṁ_L(t) = ṁ0
# DOMAINS
L = L0
doms = [t ∈ Interval(0.0, 30*3600), x ∈ Interval(0, L)]
# (INITIAL &) BOUNDARY CONDITIONS
bcs = [ p(0,x) ~ p0,
ṁ(0,x) ~ ṁ0,
p(t,0) ~ p_0(t),
ṁ(t,L) ~ ṁ_L(t)]
# SYMBOLIC MODEL
@named pipe = PDESystem(eqs,bcs,doms,[t,x],[p(t,x), ṁ(t,x)],
[RdM=>RdM0, T=>T0, A_p=>Ap0, g=>g0, H=>H0, L=>L0, ℘=>℘0, f=>f0])
# The above leads to warnings, but nothing that stops the execution
# DISCRETIZATION INTO NUMERIC MODEL
Nmol = 10
order = 2
disc = MOLFiniteDifference([x=>Nmol], t, approx_order=order)
prob = discretize(pipe,disc)
# The above `discretize` function runs for 1 minute ++ on my 12th gen intel laptop, then crashes
# SOLVING MODEL AND EXTRACTING RESULTS
sol = solve(prob)
solp = sol[p(t,x)]
solṁ = sol[ṁ(t,x)]
solx = sol[x]
solt = sol[t]
# Code crashes before this headline
xtalax
bgctw
ChrisRackauckas
bolognam
sathvikbhagavan
Here is my model:
This 1D model describes the flow of ideal gas in a pipe. I assume that inlet pressure $p_0$ is a known function of time, and that the outlet mass flow rate $\dot{m}_L$ is a known function of time.
I get a few
warningsafter the statement:but the code still runs. However, when I then do:
I get more warnings, and eventually the running crashes.
Question: Is there a problem with the term:
?? Is only the differentiation of a simple function allowed?
If so, how do I get around it?
If you guys need it, I can provide the entire code.