Implementing Maple water hammer example (PDE's)

Specific Domains Modelling & Simulations
1

Hi

I’m trying to recreate this Maple water hammer example in Julia. Link to Maple example below. I’m not able to get anything close to stable results. Some help here would be greatly appreciated. (I know I’ve made some simplifications, but can’t see that this should effect the results). Pretty new to Julia so there is probably a lot of improvments to be made.

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets, Plots

@parameters x y t
@variables P(..) V(..)
Dt = Differential(t)
Dx = Differential(x)

friction(V) = 64/((0.1*V*1000)/0.001)

x_min = y_min = t_min = 0.0
x_max = 25.0
t_max = 3.0

α = 10.

Ks = 200000000
p = 1000
Psource = 755147.8914031894

u0(x,y,t) = 22(y*(1-y))^(3/2)
v0(x,y,t) = 27(x*(1-x))^(3/2)
P0(x,t) = ((25-x)/25)*Psource
V0(x,t) = (14.19058741*(t<=2)) +(14.19058741^(-70*(t-2)) * (t>2))

eq = [Dt(V(x,t)) + (1/p)*Dx(P(x,t)) + ((0.03*abs(V(x,t))*V(x,t))/(2*0.1)) ~ 0, 
    Dx(V(x,t)) + (1/Ks)+Dt(P(x,t)) ~ 0]

domains = [x ∈ Interval(x_min, x_max),
              t ∈ Interval(t_min, t_max)]

# Periodic BCs
bcs = [V(x,0) ~ 14.19058741, # t=0
       V(x_max,t) ~ V0(x,t),
     
       P(x,0) ~ P0(x,t),
       P(0,t) ~ Psource]
    #   P(x_max,t) ~ 0] 

@named pdesys = PDESystem(eq,bcs,domains,[x,t],[P(x,t),V(x,t)])

N = 32

order = 2 # This may be increased to improve accuracy of some schemes

# Integers for x and y are interpreted as number of points. Use a Float to directtly specify stepsizes dx and dy.
discretization = MOLFiniteDifference([x=>N], t, approx_order=order)

# Convert the PDE problem into an ODE problem
println("Discretization:")
@time prob = discretize(pdesys,discretization)

println("Solve:")
@time sol = solve(prob, TRBDF2(), saveat=0.1)

plot(sol[V(x,t)][:,26]) ## Plotting V at time 26
2

Have you tried using the WENO scheme, and/or an implicit solver such as FBDF or QBDF?

3

Yes, I’ve tried all combination of it. But with the WENO scheme it trows lots of errors which I’m not able to solve.

4

Hi victor,
your example seems to be interesting :slight_smile:

I am not familiar with it and have two questions:

1. Boundary Conditions

In the example description the math formulas of the boundary conditions (BCs) are missing.

As we have two states (P, V) and a 1D geometry we need 2x2=4 BCs (2 BC per side: left/right) - or only 2 if we assume periodic BCs.

I try to decode them from the Maple code:

P(x=0,t)=Psource
P(x=x_max,t) = 0

V(x=x_max,t) = Vsteady if t<2
V(x=x_max,t) = Vsteady * exp(-70(t-2))* if t>=2

About V(x=0,t) I guess they mean a (zero) Neumann BC:

dV(x=0,t)/dx = 0

These BCs do not seem to be periodic. In your code you only have 1 BC for P and S.
So you need to implement all 4 BCs in your code.

2. Friction term

In your code you implemented the friction term but you do not use it in the PDE equation.

Furthermore, your friction term seems to differ from the Maple code (-> Friction Factor) which as several (if-else) cases.

If you find your issues and you can run your simulation, can you please post it here :slight_smile: