# Water Hammer - Numerically solving system of PDEs

I'm trying to use Mathematica to solve the water hammer effect.

g = 9.81;
a = 1350;
L = 3500;
h0 = 4;
v0 = Sqrt[2 g h0];
R = 0.003;

sol = NDSolve[{
D[h[x, t], x] - R*v[x, t]*Abs[v[x, t]] == 1/g D[v[x, t], t],
D[v[x, t], x] == g/a^2*D[h[x, t], t],

v[x, 0] == v0,
v[0, t] == v0 Exp[-t^2/0.4],
h[L, t] == h0,
h[x, 0] == h0},

{h, v},
{x, 0, L}, {t, 0, 10}
];

Manipulate[
Plot[Evaluate[v[x, t] /. sol], {x, 0, L}, PlotRange -> {-2 v0, 2 v0}],
{t, 0, 10}]


What I get near the end of the time interval is something I'm not expecting:

The documentation tells me to use the option:

Method -> {"MethodOfLines","SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 750}}


but it just makes it worse.

Can somebody help me out with this one?

PS: Take R=0 and you get a lossless system, and the solution should be a wave traveling and reflecting for h and v.

You need the magic of "Pseudospectral" or a dense enough 2nd order spatial difference grid:

mol[n_Integer, o_:"Pseudospectral"] := {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n,
"MinPoints" -> n, "DifferenceOrder" -> o}}

g = 9.81;
a = 1350;
L = 3500;
T = 30;
h0 = 4;
v0 = Sqrt[2 g h0];
R = 0.003;

(* Solution 1 *)
sol = NDSolve[{D[h[x, t], x] - R v[x, t] Abs[v[x, t]] == D[v[x, t], t]/g,
D[v[x, t], x] == g D[h[x, t], t]/a^2, v[x, 0] == v0, v[0, t] == v0 Exp[-(t^2/0.4)],
h[L, t] == h0, h[x, 0] == h0}, {h, v}, {x, 0, L}, {t, 0, T}, Method -> mol[45]];

(* Solution 2 *)
sol2 = NDSolve[{D[h[x, t], x] - R v[x, t] Abs[v[x, t]] == D[v[x, t], t]/g,
D[v[x, t], x] == g D[h[x, t], t]/a^2, v[x, 0] == v0, v[0, t] == v0 Exp[-(t^2/0.4)],
h[L, t] == h0, h[x, 0] == h0}, {h, v}, {x, 0, L}, {t, 0, T}, Method -> mol[200, 2]];

(* Use sol2 inside Plot if you like *)
Manipulate[
Plot[Evaluate[v[x, t] /. sol], {x, 0, L}, PlotRange -> {-2 v0, 2 v0}], {t, 0, T}]


Velocity at the end of the pipe:

(* Use sol2 inside Plot if you like *)
Plot[Evaluate[v[L, t] /. sol], {t, 0, T}, PlotRange -> {-2 v0, 2 v0}]


• Awesome! Thank you! This works great and fast enough. Could you elaborate on the "magic" part? – Ivan Dec 12 '14 at 19:42
• @Ivan To be honest, I'm not able to… the only thing I know is that "Pseudospectral" seems to be quite effective on certain kind of PDEs (specifically speaking, PDEs related to fluid dynamics and quantum mechanics, according to the examples appearing in this site so far), but I never studied why. Maybe you can consider asking it in scicomp.stackexchange.com ? – xzczd Dec 15 '14 at 5:25

With the option MaxStepSize -> 1., it seems to work. (A little bit magic)

sol = NDSolve[{
D[h[x, t], x] - R*v[x, t]*Abs[v[x, t]] == 1/g D[v[x, t], t],
D[v[x, t], x] == g/a^2*D[h[x, t], t],

v[x, 0] == v0,
v[0, t] == v0 Exp[-t^2/0.4],
h[L, t] == h0,
h[x, 0] == h0},

{h, v},
{x, 0, L}, {t, 0, 10},
MaxStepSize -> 1.
];


Here are the step sizes with the option MaxStepSize -> 1. :

Of course it is a way below 1...

... and It is finer that with the default MaxStepSize :

In the comments, @belisarius ask a plot of the velocity at the end of the pipe. Here it is :

• @belisarius This is consistent with the simulation above. At the end of the line (x=L) there is a oscillation. The caracteristics of the oscillation depend of the length of the line. – andre314 Dec 10 '14 at 20:23
• I'm not sure, though. I think the fluid should come to a stop due to the dissipation. – Dr. belisarius Dec 10 '14 at 20:49
• @belisarius This is precisely what happens : when you move the cursor you see that the wave is vanishing. – andre314 Dec 10 '14 at 20:55
• Ok, I can't run this in my machine (takes too much time). Perhaps you could include a plot of the velocity as a function of time at the end of the pipe. – Dr. belisarius Dec 10 '14 at 21:00
• @belisarius I have added a plot of the velocity. – andre314 Dec 10 '14 at 21:09