# Implementing Maple water hammer example (PDE's)

**URL:** <https://discourse.julialang.org/t/implementing-maple-water-hammer-example-pdes/103501>\
**Category:** Modelling & Simulations\
**Tags:** pde\
**Created:** [September 4, 2023, 9:58am UTC](https://discourse.julialang.org/t/implementing-maple-water-hammer-example-pdes/103501 "2023-09-04T09:58:24Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![victor](https://avatars.discourse-cdn.com/v4/letter/v/3e96dc/32.png) [@victor](https://discourse.julialang.org/u/victor)\
**Post date:** [September 4, 2023, 9:58am UTC](https://discourse.julialang.org/t/implementing-maple-water-hammer-example-pdes/103501/1 "2023-09-04T09:58:24Z")

</div>

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.

> **[Water Hammer - Maple Help](https://www.maplesoft.com/support/help/maple/view.aspx?path=applications%2FWaterHammer)**
>
> Water Hammer When a valve at the end of a pipeline suddenly closes, a pressure surge hits the valve and travels along the pipeline. This is known as water hammer. This process is modeled by two partial differential equations (PDEs). The PDEs can be...

```julia
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

```

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [September 4, 2023, 12:20pm UTC](https://discourse.julialang.org/t/implementing-maple-water-hammer-example-pdes/103501/2 "2023-09-04T12:20:49Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![victor](https://avatars.discourse-cdn.com/v4/letter/v/3e96dc/32.png) [@victor](https://discourse.julialang.org/u/victor)\
**Post date:** [September 5, 2023, 3:16pm UTC](https://discourse.julialang.org/t/implementing-maple-water-hammer-example-pdes/103501/3 "2023-09-05T15:16:51Z")

</div>

> [@xtalax](#):
>
> ENO scheme, and/or an

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.

---

<div class="post-metadata">

**Author:** ![stephans3](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stephans3/32/34552_2.png) [@stephans3](https://discourse.julialang.org/u/stephans3)\
**Post date:** [September 14, 2023, 10:00am UTC](https://discourse.julialang.org/t/implementing-maple-water-hammer-example-pdes/103501/4 "2023-09-14T10:00:08Z")

</div>

Hi victor,  
your example seems to be interesting 🙂

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 🙂