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