Method of characteristics, recursion formulas

Question

I'm helping a friend out with a problem. He needs to plot out the solution to a water surge in a pipeline.

I have all the parameters set and I've put the equations in Mathematica, but I can't get the table that will allow me to use ListLinePlot to plot the data.

The problem goes like this.

The water is flowing through a horizontal pipeline from a reservoir. The pipeline is L = 1 km long and has a d = 500 mm internal diameter. If the valve closes in 4 seconds (and gives a linear retardation), estimate the pressure rise at the valve and at the mid point over an 8 second interval. The reservoir head is $H_1^0$ = 100 m and at the head at the valve before commencement of valve closure the head is $H_3^0$ = 4.75 m. Friction factor is $\lambda$ = 0.012, and the speed of sound is c = 1000 m/s.

I've set up the parameters like this:

L = 1000;
d = 0.5;
\[CapitalDelta]t = 0.5;
H10 = 100;
H30 = 4.75;
\[Lambda] = 0.012;
c = 1000;
g = 9.81;

Initial conditions: at t=0. Time interval is 0.5 s - the time necessary for the shock wave to get to the midpoint. We can find the water flow speed (u) at the valve:

$H_1^0-H_3^0=\frac{\lambda L (u_3^0)^2)}{2g d}$

In Mathematica the formula is:

u30 = Sqrt[(2 g d (H10 - H30))/(\[Lambda] L)]

8.82422

This is the speed of water flow at the valve. In the t=0, when the water is just flowing freely, the speed of the flow is the same along the pipeline and on the beginning. So $u_1^0=u_2^0=u_3^0$.

The numbers are indices in the formulas used to calculate speed and the head at the pipeline intake (i=1) midpoint (i=2) and at the valve (i=3). The upper indices are time intervals.

The solution is basically recursively finding the values for the speed and the head at the increasing time intervals at the midpoint and at the valve.

So the initial data, where H10 is the head at the pipeline beginning, H20 head at the midpoint, and H30 is the head at the valve is:

H10 = 100;
H20 = 52.38;
H30 = 4.75;

The head at the t=0 in the midpoint is

H20 = H30 + (\[Lambda]*(L/2)*u30^2)/(2 g d)

The characteristic coefficient C is c/g=101.937.

If we take the first step to be t = 0.5 s (and the steps can be altered - more steps, more points, and better the plot). We have valve area reduction (VAR) = 0.5/4 = 0.125, and the factor $A_R = 1- VAR = 0.875$. This factor changes as the time interval changes (so I can make a list of these values).

Now the equations look like this. For the valve the speed is calculated with this formula:

$u_{i+1}^{n+1}=-\frac{1}{2}\frac{C(u_{i+1}^0)^2A_R^2}{H_{i+1}^0}+\frac{1}{2} \bigg[\frac{C^2(u_{i+1}^0)^4 A_R^4}{(H_{i+1}^0)^2}-4\frac{(u_{i+1}^0)^2A_R^2}{H_{i+1}^0}\bigg(\frac{C\lambda u_i^n |u_i^n| }{2d}\delta t-Cu_i^n-H_i^n\bigg) \bigg]^{\frac{1}{2}}$

$i+1$ is the next point, so if the $i=2$ (midpoint) $i+1 = 3$ is the valve. The $n$ in the superscript denotes the time. So if the $n$ is the current time interval, $n+1$ will be the next time interval. The formula for the head is

$H_{i+1}^{n+1}-H_i^n+C(u_{i+1}^{n+1}-u_i^n)+\frac{C\lambda u_i^n|u_i^n|}{2d}\delta t=0$

Now if we substitute the valve point and midpoint (i=2, and i+1=3) we get the recursion formulas for the valve:

$u_{3}^{n+1}=-\frac{1}{2}\frac{C(u_{3}^0)^2A_R^2}{H_{3}^0}+\frac{1}{2} \bigg[\frac{C^2(u_{3}^0)^4 A_R^4}{(H_{3}^0)^2}-4\frac{(u_{3}^0)^2A_R^2}{H_{3}^0}\bigg(\frac{C\lambda u_2^n |u_2^n| }{2d}\delta t-Cu_2^n-H_2^n\bigg) \bigg]^{\frac{1}{2}}$

That is, if we take n=0 we can use the initial conditions

$u_{3}^{1}=-\frac{1}{2}\frac{C(u_{3}^0)^2A_R^2}{H_{3}^0}+\frac{1}{2} \bigg[\frac{C^2(u_{3}^0)^4 A_R^4}{(H_{3}^0)^2}-4\frac{(u_{3}^0)^2A_R^2}{H_{3}^0}\bigg(\frac{C\lambda u_2^0 |u_2^0| }{2d}\delta t-Cu_2^0-H_2^0\bigg) \bigg]^{\frac{1}{2}}$

That gives 8.81015 m/s which is correct. Also the head can be then calculated as

$H_{3}^{1}=H_2^0-C(u_{3}^{1}-u_2^0)-\frac{C\lambda u_2^0|u_2^0|}{2d}\delta t$

And that gives 6.18431 m, also correct.

I've managed to create a recursion formula for this in Mathematica like

u3[y_] := -((C' (u30)^2 (AR[[n]])^2)/(2*H30)) + 1/2 (((C')^2 (u30)^4 (AR[[n]])^4)/(H30)^2 - (4 (u30)^2 (AR[[n]])^2)/H30 ((C' \[Lambda] u2[y - 1] Abs[u2[y - 1]])/(2 d) \[Delta]t[[n]] - C' u2[y - 1] - H2[y - 1]))^(1/2)

I've created my time intervals and $A_R$ like

\[Delta]t = Table[i, {i, \[CapitalDelta]t, 8, \[CapitalDelta]t}]
VAR = Table[i/4, {i, \[CapitalDelta]t, 8, \[CapitalDelta]t}]
AR = 1 - VAR
C' = c/g

The problem I'm facing is that when ever I try to evaluate this, with u3[1]

I get the error:

The expression n cannot be used as a part specification. >>

-835.526 {0.875, 0.75, 0.625, 0.5, 0.375, 0.25, 0.125, 0., -0.125, -0.25, -0.375, -0.5, -0.625, -0.75, -0.875, -1.}[[ n]]^2 + 1/ 2 [Sqrt](2.79242*10^6 {0.875, 0.75, 0.625, 0.5, 0.375, 0.25, 0.125, 0., -0.125, -0.25, -0.375, -0.5, -0.625, -0.75, -0.875, \ -1.}[[n]]^4 - 65.5721 {0.875, 0.75, 0.625, 0.5, 0.375, 0.25, 0.125, 0., -0.125, -0.25, -0.375, -0.5, -0.625, -0.75, -0.875, -1.}[[ n]]^2 (-H2[0] - 101.937 u2[0] + 1.22324 Abs[ u2[0]] {0.5, 1., 1.5, 2., 2.5, 3., 3.5, 4., 4.5, 5., 5.5, 6., 6.5, 7., 7.5, 8.}[[n]] u2[0]))

When I evaluate this I get the correct result (a number), with the same error.

u3[1] /. n -> 1 /. u2[0] -> u30 /. H2[0] -> H20

Now I can create a table of values, and then I won't get an error, but I'll get a table with the u2[y-1] in them, and since I need to increase by 1 every step this isn't a good way to go about it. Also for the higher terms, I need the values of speed and head at the midpoint which I get from the fromulas:

$u_{i+1}^{n+1}=\frac{1}{2}\bigg[\frac{H_i^n-H_{i+2}^n}{C}+u_i^n+u_{i+2}^n-\frac{\lambda}{2d}(u_i^n |u_i^n|+u_{i+2}^n |u_{i+2}^n|)\delta t\bigg]$

$H^{n+1}_{i+1}=\frac{1}{2}\bigg[H^n_i+H^n_{i+2}-C(u_{i+2}^2-u_i^n)-\frac{C\lambda}{2d}(u_i^n |u_i^n|-u_{i+2}^n |u_{i+2}^n|)\delta t\bigg]$

So I need to input the values I got from the previous step in this, to get the values at the midpoint. These values I'll be able to input in the formula for the valve for the next time step, and so on, until I get to the end of my interval. I've created the recursive formulas for these two, but I get max recursion depth error, and everything falls apart :\

u2[y_] :=  1/2 ((H10 - H3[y - 1])/C' + u10 + u3[y - 1] - \[Lambda]/(2 d) (u10 Abs[u10] + u3[y - 1] Abs[u3[y - 1]]) \[Delta]t[[n]])

H2[y_] := 1/2 (H10 + H3[y - 1] - C' (u3[y - 1] - u10) - (\[Lambda] C')/(2 d) (u10 Abs[u10] - u3[y - 1] Abs[u3[y - 1]]) \[Delta]t[[n]])

I need to end up with a table of Head values for the midpoint and the valve, and speeds, and then I'd just plot that against my time interval. Also I'd be able to refine this, to a time step of 0.1 s instead of 0.5 s.

But I'm stuck. I've solve this by hand for time steps 0.5 s and 1 s, but calculating each step by hand is not practical.

I know this is a wall of text, but I tried to be as clear as possible. I haven't used Mathematica in a while. Any help is appreciated.