Exact computation of the axial vibration of two coupled liquid-filled pipes

Exact computation of the axial vibration of two coupled

liquid-filled pipes

Citation for published version (APA):

Tijsseling, A. S. (2009). Exact computation of the axial vibration of two coupled liquid-filled pipes. (CASA-report; Vol. 0914). Technische Universiteit Eindhoven.

Document status and date: Published: 01/01/2009

Document Version:

Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers)

Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne

Take down policy

If you believe that this document breaches copyright please contact us at:

openaccess@tue.nl

providing details and we will investigate your claim.

EINDHOVEN UNIVERSITY OF TECHNOLOGY

Department of Mathematics and Computer Science

CASA-Report 09-14

May 2009

Exact computation of the axial vibration

of two coupled liquid-filled pipes

by

A.S. Tijsseling

Centre for Analysis, Scientific computing and Applications

Department of Mathematics and Computer Science

Eindhoven University of Technology

P.O. Box 513

5600 MB Eindhoven, The Netherlands

ISSN: 0926-4507

1

EXACT COMPUTATION OF THE AXIAL VIBRATION OF TWO COUPLED

LIQUID-FILLED PIPES

Arris S. Tijsseling

Department of Mathematics and Computer Science

Eindhoven University of Technology

P.O. Box 513, 5600 MB Eindhoven

The Netherlands

Phone: +31 40 247 2755

Fax: +31 40 244 2489

E-Mail:

a.s.tijsseling@tue.nl

ABSTRACT

A simple recursion is presented that finds exact solutions to the problem of two coupled axially-vibrating liquid-filled pipes. Fluid-structure interaction at pipe ends and junction, and along the pipe because of axial-radial Poisson contraction, is taken into account. The solutions obtained for a waterhammer problem show unprecedented details that resemble noise.

Keywords waterhammer, pipe vibration, fluid-structure interaction, travelling wave, method of characteristics INTRODUCTION

Strong fluid-structure interaction (FSI) may happen in axially vibrating liquid-filled pipes with closed ends. At the pipe ends, pressure waves in the liquid interact with stress waves in the pipe walls. Such interaction takes also place at unrestrained junctions like pipe diameter changes and elbows (junction coupling). Second and third types of interaction occur along the entire pipe as a result of axial-radial contraction effects (Poisson coupling) and skin friction (friction coupling). The significance of the subject lies in the study of severe waterhammer events and in the altered resonant vibration of liquid-filled pipes [1-3].

A general method is presented that finds exact solutions of the sketched problem if the wave propagations are one-dimensional and non-dispersive. Wave fronts are automatically traced back in time to a given initial situation with the aid of a simple recursion. There are no interpolations, no adjustments of system parameters and no numerical approximations. The strength of the method is the simplicity of the algorithm (recursion) and the accurate (exact) computation of transient events. Its weakness is that the

computation time increases exponentially for events of longer duration.

The test problem concerns the axial vibration of two connected, different pipes, where at the junction one incident wave results in two reflected and two transmitted waves. Instantaneous valve closure generates a waterhammer impact load of the pipes. The new exact solution of the transient vibration serves as a benchmark in the verification of numerical results and schemes. Furthermore, the method can be used to perform parameter variation studies without parameter changes generated by the numerical scheme itself (e.g. changed wave speeds or pipe lengths). The present study is a combination and extension of earlier work [4-5]. It gives for the first time exact solutions for FSI in a two-pipe system and as such is a valuable addition to the single-pipe studies [5-6].

MATHEMATICAL MODEL

Skalak [7] presented a four-equation model for the low-frequency acoustic vibration of straight, thin-walled, linearly elastic, prismatic pipes of circular cross-section, filled with inviscid liquid. The four equations, governing fluid velocity, V, fluid pressure, P, axial pipe velocity, u , and axial pipe stress, _z

z σ , are 1 0 f V P t ρ z ∂ ∂ + = ∂ ∂ , (1)

2 1 2 2 ( ) z 0 V R P z K E e t E t ν σ ∂ ∂ ∂ + + − = ∂ ∂ ∂ , (2) 1 0 z z s u t z σ ρ ∂ ∂ − = ∂ ∂  , (3) 1 0 z z R P u z E t E e t ν σ ∂ ∂ ∂ − + = ∂ ∂ ∂  , (4) where E = Young modulus, e = wall thickness, K = bulk modulus,

R = inner pipe radius, t = time, z = distance along pipe, ν =

Poisson ratio, ρ = mass density; the subscripts f and s refer to fluid and structure, respectively.

Initial conditions for all four variables, and two boundary conditions at each end of the pipe, complete the mathematical model.

If the Poisson ratio is set equal to zero (ν = ) the Eqs (1-2) 0 and (3-4) are reduced to classical waterhammer and longitudinal beam equations, respectively. Still, they may be coupled via mutual boundary conditions.

The Eqs (1-4) themselves are energy conserving, but damping can be introduced through the boundary conditions, for example an orifice for the fluid and a dashpot for the structure.

RIEMANN INVARIANTS

Method-of-characteristics (MOC) analysis of the basic Eqs (1-4) reveals the four Riemann invariants

2 1 1 1 _{2 2 2} 1 2 1 1 2 2 2 2 1 1 2 2 2 , s f f s z z s s s R V P e c c u c c λ _ν λ η ρ ρ λ λ _ν λ ν σ ρ λ λ ⎛ ⎞ ⎜ ⎟ = + + + ⎜ ₋ ⎟ ⎝ ⎠ ⎛ ⎞ ⎛ ⎞ ⎜ ⎟ ⎜ ⎟ + − ⎜ ₋ ⎟ ⎜ ₋ ⎟ ⎝ ⎠ ⎝ ⎠  (5) 2 2 2 2 _{2 2 2} 2 2 2 2 2 2 2 2 2 2 2 2 2 , s f f s z z s s s R V P e c c u c c λ _ν λ η ρ ρ λ λ _ν λ ν σ ρ λ λ ⎛ ⎞ ⎜ ⎟ = + + + ⎜ ₋ ⎟ ⎝ ⎠ ⎛ ⎞ ⎛ ⎞ ⎜ ⎟ ⎜ ⎟ + − ⎜ ₋ ⎟ ⎜ ₋ ⎟ ⎝ ⎠ ⎝ ⎠  (6) 2 2 2 3 3 3 2 2 2 2 3 3 2 3 3 2 2 , f f f f f z z s s s c c R R V P Eec Eec u c c λ λ ν ν η ρ λ λ λ λ σ ρ ⎛ ⎞ ⎛ ⎞ ⎜ ⎟ ⎜ ⎟ = _{⎜ ⎟} +_{⎜ ⎟} + − − ⎝ ⎠ ⎝ ⎠ ⎛ ⎞ _{⎛ ⎞} ⎜ ⎟ +_{⎜ ⎟} − ⎜_⎜ ⎟_⎟ ⎝ ⎠ ⎝ ⎠  (7) 2 2 2 4 4 4 _{2 2 2 2} 4 4 2 4 4 2 2 , f f f f f z z s s s c c R R V P Eec Eec u c c λ λ ν ν η ρ λ λ λ λ σ ρ ⎛ ⎞ ⎛ ⎞ ⎜ ⎟ ⎜ ⎟ = + + ⎜ ₋ ⎟ ⎜ ₋ ⎟ ⎝ ⎠ ⎝ ⎠ ⎛ ⎞ ⎛ ⎞ ⎜ ⎟ +_{⎜ ⎟} − ⎜_⎜ ⎟_⎟ ⎝ ⎠ ⎝ ⎠  (8)

which are constant along the characteristic lines in the distance-time plane with directions d /dz t=λ₁, d /dz t=λ₂, d /dz t=λ₃

and d /dz t=λ₄, respectively. The wave propagation speeds are 1/2 1/2 2 4 2 2 1, 2 1 ( 4 ) 2 f s c c λ = ± ⎡_⎣γ − γ − ⎤_⎦ and 1/2 1/2 2 4 2 2 3, 4 1 ( 4 ) 2 c c f s λ _{= ±} ⎡γ ₊ γ ₋ ⎤ ⎣ ⎦ , (9)

where the constants

1/2 2 2 (1 ) f f f R c K Ee ρ ρ ν − ⎡ ⎤ = _⎢ + − _⎥ ⎣ ⎦ and 1/2 s s E c ρ ⎛ ⎞ = ⎜_⎜ ⎟_⎟ ⎝ ⎠ (10) are the classical pressure and axial-stress wave speeds, and

2 _{( 1 2} 2 f ₎ 2 2 f s s R c c e ρ γ ν ρ = + + . (11)

A full derivation of the wave speeds and the Riemann invariants can be found in [5 and 8].

Algebraic inversion of the Eqs (5-8) yields

2 2 1 1 2 2 2 2 2 2 3 1 2 2 1 3 4 2 2 2 1 3 1 2 2 2 , 2 f f s f s s s R c V e c c c c λ η η ρ ν ρ λ λ λ η η ν λ λ ⎛ ⎞ ⎛ + ⎞ ⎜ − ⎟ = ⎜ ⎟+ ⎜ _{− −} ⎟ _{⎝ ⎠} ⎝ ⎠ + ⎛ ⎞ − ⎜ ⎟ − _{⎝ ⎠} (12) 2 2 2 3 1 2 2 2 2 2 2 1 3 1 2 2 3 4 2 2 3 1 1 2 2 2 , 2 f f f f s f s f f s s R c c P e c c c c c λ η η ρ ρ ν ρ λ λ λ η η ρ ν λ λ ⎛ ⎞ ⎛ − ⎞ ⎜ − ⎟ = ⎜ ⎟+ ⎜ _{− −} ⎟ _{⎝ ⎠} ⎝ ⎠ − ⎛ ⎞ − _{⎜ ⎟} − _{⎝ ⎠} (13) 2 2 2 1 1 2 2 2 2 2 2 2 2 3 1 3 2 3 4 2 3 1 2 2 , (14) 2 f f f f z s f s s f s R c R c u e c c e c c λ η η ρ ρ ν ν ρ λ λ ρ λ η η λ ⎛ ⎞ ⎛ + ⎞ ⎜ − ⎟ = − ⎜ ⎟+ ⎜ _{− −} ⎟ _{− ⎝ ⎠} ⎝ ⎠ + ⎛ ⎞ + ⎜ ⎟ ⎝ ⎠  2 2 2 2 3 1 2 2 2 2 2 2 2 2 1 3 1 3 2 2 3 4 2 2 2 3 1 1 2 2 1 2 . (15) 2 f f f f f z s f s f f f s s s s R c R c c e c c e c R c c e c λ η η ρ ρ ν σ ν ρ λ λ λ λ η η ρ ρ ν ρ λ λ ⎛ ⎞ ⎛ − ⎞ ⎜ − ⎟ = ⎜ ⎟+ ⎜ _{− −} ⎟ _{− ⎝ ⎠} ⎝ ⎠ ⎛ ⎞ ⎛ − ⎞ ⎜ ⎟ − +_{⎜ ⎟} ⎜ ⎟ − _{⎝ ⎠} ⎝ ⎠

The Eqs (5-9) and (12-15) are summarised in matrix-vector notation by 1 ( , ) z t = S− ( , )z t η φ and ( , ) φ z t = S

η

( , )z t , (16) where T 1 2 3 4 ( , , , )η η η η = η and _{( , , , )}T z z V P u σ =  φ .

3 SOLUTION ALGORITHM FOR A SINGLE PIPE

The analysis is in terms of the Riemann invariants η 1, η 2, η 3

and η 4. Therefore, the initial and boundary conditions for V, P,

z

u and σ_z are translated into Riemann invariants by means of the

Eqs (5-8). The solution T 1 2 3 4 ( , , , )η η η η =

η in interior point P = (z, t), 0 < z < L, 0 < t, in the distance-time plane, is determined by the values of η 1, η 2, η 3 and η 4 generated at the boundary points

A 1, A 2, A 3 and A 4, respectively, as displayed in Fig. 1. The value

of η in the boundary points, at z = 0 or z = L, is calculated from two compatibility equations, (5, 7) or (6, 8), and two given linear boundary conditions. Because each boundary calculation is mathematically the same, the entire solution can be condensed into a simple recursion. See Appendix A. Tracing characteristic lines backwards in time, the recursion stops at time level t = 0, where η is assigned a known initial value. The general method is fully described and verified in Ref. [5]. At last, solution ( , )η z t is

decomposed into the original variables through the Eqs (12-15).

Figure 1. Distance-time plane for single pipe.

SOLUTION ALGORITHM FOR A DOUBLE PIPE

If two different pipes (1 and 2) are connected at z = zj , the

local solution vector η (P) is discontinuous and split into η1(P+_{) and η2 (P}−_{). The compatibility relations (5-8) provide}

four equations for the eight unknowns in η (P±_{), as seen from}

Fig. 2, and four additional equations are given by the linear junction condition j j ( )t (z−, )t = ( )t (z+, )t + ( )t J1 φ J2 φ r or j j ( )t (z t, ) = ( )t (z t, ) + ( )t J1 S1 η1 J2 S2 η2 r , (17) where J1 and J2 are 4 by 4 coefficient matrices and the vector r represents for instance excitation or leakage. The pipe junction is without wave reflection only if J1( )t S1 = J2( )t S2 and

( )t =

r 0 .

Conservation of volume and equilibrium of forces at an axial junction of two pipes, with cross-sectional areas A1 and A2 of fluid (f) and solid (s), is described herein by the two matrices:

1 0 1 0 0 1 0 0 0 0 1 0 0 1 0 1 f f f s A A A A − ⎛ ⎞ ⎜ ⎟ ⎜ ⎟ = ⎜ ⎟ ⎜ ⎟ ⎜ ₋ ⎟ ⎝ ⎠ J1 , 2 0 2 0 0 1 0 0 0 0 1 0 0 2 0 2 f f f s A A A A − ⎛ ⎞ ⎜ ⎟ ⎜ ⎟ = ⎜ ⎟ ⎜ ⎟ ⎜ ₋ ⎟ ⎝ ⎠ J2 and ( )r t =0 . (18)

Separating the unknown values T

1 2 3 4 ( 2 , 1 , 2 , 1 )η η η η =

u η

from the known values T

1 2 3 4 ( 1 , 2 , 1 , 2 )η η η η =

k

η at the junction,

see Fig. 2, gives the following relation: 1 − = u JSU JSK k η η (19) with 11 12 13 14 21 22 23 24 31 32 33 34 41 42 43 44 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) − − ⎛ ⎞ ⎜_{− −} ⎟ ⎜ ⎟ =_{⎜ ⎟} − − ⎜ ⎟ − − ⎝ ⎠ J1 S1 J2 S2 J1 S1 J2 S2 J1 S1 J2 S2 J1 S1 J2 S2 JSK J1 S1 J2 S2 J1 S1 J2 S2 J1 S1 J2 S2 J1 S1 J2 S2 and 11 12 13 14 21 22 23 24 31 32 33 34 41 42 43 44 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) − − ⎛ ⎞ ⎜_{− −} ⎟ ⎜ ⎟ = ⎜− − ⎟ ⎜ ⎟ − − ⎝ ⎠ J2 S2 J1 S1 J2 S2 J1 S1 J2 S2 J1 S1 J2 S2 J1 S1 JSU J2 S2 J1 S1 J2 S2 J1 S1 J2 S2 J1 S1 J2 S2 J1 S1 . (20)

Symbolic matrix inversion of the 4 by 4 matrix JSU is possible, but not carried out herein because there is no clear advantage over a numerical inversion. Because J1 and J2 in Eq. (18) are constant matrices, JSU has to be inverted numerically only once.

The solution procedure is tracking wave paths back in time and solving the boundary and junction equations in a recursive way as described for the single pipe. The junction coupling equation (19) is now included in the recursion. The algorithm is given in Appendix A.

4 Figure 3. Sketch of test problem; PT = pressure transducer; distances in metres.

Table 1. Input data for the calculation: material and geometrical properties of liquid and pipes.

Liquid Pipe 1 Pipe 2

K = 2.1 MPa ρf = 1000 kg/m 3 Q0 = 0.499 m3/s L1 = 10 m R1 = 398.5 mm e1 = 16 mm E1 = 2.1 GPa ρ1 s = 7900 kg/m 3 ν1 = 0.30 L2 = 10 m R2 = 398.5 mm e2 = 8 mm E2 = 2.1 GPa ρ2 s = 7900 kg/m 3 ν2 = 0.30

Table 2. Initial liquid velocities (V0), classical wave speeds (c), coupled wave speeds (λ), and wave travel times (Δ t).

Pipe 1 Pipe 2 Relations Wave travel times

V10 = 1 m/s V20 = 1 m/s A1 f V10 = A2 f V20 = Q0 L1/ V10 = L2/ V20 = 10 s c1 f = 1202.1 m/s c1 s = 5155.8 m/s λ11 = 1183.0 m/s λ13 = 5239.1 m/s c2 f = 1049.5 m/s c2 s = 5155.8 m/s λ21 = 1024.7 m/s λ23 = 5280.5 m/s c1 f > c2 f c1 s = c2 s λ11 > λ21 λ13 < λ23 c1 f > λ1 1 c1 s < λ1 3 c2 f > λ2 1 c2 s < λ2 3 Δ t11 = L1/λ11 = 8.453 ms Δ t21 = L2/λ21 = 9.759 ms Δ t13 = L1/λ13 = 1.909 ms Δ t23 = L2/λ23 = 1.894 ms CALCULATED RESULTS

The test problem is instantaneous valve closure in the reservoir-pipe-valve system shown in Figure 3. The pipe is 20 m long and divided into two 10 m sections (1 and 2) with different wall thicknesses as specified in Table 1. The wave speeds in each section are different as listed in Table 2 and the used boundary conditions are given in Table 3. Wave reflection will occur at the junction of the two sections. Liquid waves will reflect because of the different wave speeds in the sections and solid waves will primarily reflect because of the difference in cross-sectional wall area. As a result of FSI (junction coupling), one incident wave produces two reflected and two transmitted waves. At the pipe ends, one incident wave generates two reflected waves: at z = 0 (fixed reservoir) because of Poisson coupling and at z = L (moving massless valve) because of junction coupling. The continuous double reflection brings about an exponentially growing number of wave fronts travelling in the system.

The double-pipe algorithm has been tested and successfully verified against two-pipe waterhammer (no junction and Poisson coupling [4]), one-pipe FSI (no junction [5]), and Fortran MOC results obtained with slightly modified wave speeds [9].

Table 3. Boundary conditions. boundary z = 0 z = L no coupling z = L FSI coupling fluid P1 = 0 V2 = 0 pipe u1z =0 u2z =0 2 2_z V =  u 2_f 2 2_s 2_z A P =A σ

Time histories for pressure (P), fluid velocity (V), axial stress (σ_z) and pipe velocity (u_z) calculated at z = 11.15 m are

shown in the Figs 4-7. Each graph consists of 1000 points, and results without FSI coupling (frictionless classical waterhammer) are given as a reference.

Figure 4 shows that the maximum pressure predicted by classical waterhammer theory (dashed line) is exceeded considerably because of FSI (continuous line). Furthermore, a precursor travelling roughly at speed cs is running ahead of the

main waterhammer wave which travels roughly at speed cf. The

variations in fluid velocity in Fig. 5, in particular the times of the jumps, are consistent with the variations in pressure in Fig. 4.

5 The axial pipe stress in Fig. 6 has a large positive value as long as

the pressure acting on the free closed end at z = L is large and positive. In principle, the entire solution can be calculated and explained from a detailed wave paths diagram and corresponding jump conditions [10-12]. This tremendous task is carried out automatically by the simple recursion given in Appendix A.

0 0.035 0.07 1.5 .106 0 1.5 .106 time (s) pressure (Pa)

Fig. 4. Pressure history at z = 11.15 m. Continuous line: with FSI. Dashed line: classical waterhammer.

0 0.035 0.07 1.5 0 1.5 time (s) fluid velocity (m/s)

Fig. 5. Fluid velocity history at z = 11.15 m. Continuous line: with FSI. Dashed line: classical waterhammer.

0 0.035 0.07

5 .107 0 5 .107

time (s)

axial stress (Pa)

Fig. 6. Stress history at z = 11.15 m. Calculation with FSI.

0 0.035 0.07 0.5 0 0.5 time (s) pipe velocity (m/s)

Fig. 7. Pipe velocity history at z = 11.15 m. Calculation with FSI.

The Figs 4-7 show an amazing amount of detail. When time increases the calculated traces give the impression to contain noise, but that is not the case: the solutions are exact. In fact, to detect all tiny details, the resolution (the number of points in the graphs) should increase when time advances. This is easy to do, because the method is mesh-free: z and t can be chosen arbitrarily.

The Figs 8-11 show spatial distributions at time t = 0.03 s. The main pressure wave, generated by the valve closure at t = 0, has reflected from the reservoir and is moving back to the valve (to the right in Figs 8-9). The axial stress in Fig. 10 jumps by a factor two at the junction at z = L/2, because the wall areas in pipes 1 and 2 differ by a factor two. The pipe velocity in Fig. 11 is continuous at z = L/2. 0 10 20 5 .105 0 5 .105 1 .106 1.5 .106 distance (m) pressure (Pa)

Fig. 8. Pressure distribution at t = 0.03 s. Continuous line: with FSI. Dashed line: classical waterhammer.

6 0 10 20 1.5 1 0.5 0 0.5 distance (m) fluid velocity (m/s)

Fig. 9. Fluid velocity distribution at t = 0.03 s. Continuous line: with FSI. Dashed line: classical waterhammer.

0 10 20

1 .107 2 .107 3 .107

distance (m)

axial stress (Pa)

Fig. 10. Stress distribution at t = 0.03 s. Calculation with FSI.

0 10 20 0 0.1 0.2 distance (m) pipe velocity (m/s)

Fig. 11. Pipe velocity distribution at t = 0.03 s. Calculation with FSI.

DISCUSSION

The mathematical challenge was to obtain exact solutions in a system with an exponentially growing number of travelling discontinuities. This has been achieved by means of a simple recursion. The solutions presented for waterhammer with FSI in

a two-pipe system show an amazing amount of detail, which looks like noise, but is not. In fact, the solutions mimic reality, where small imperfections in otherwise clean laboratory experiments can cause noise-like signals.

The amount of detail in the solutions increases with time, and this calls for non-uniform resolution in the graphs. This can easily be accomplished, even during the calculation, because the method is mesh-free and ideal for adaptivity.

The method is not suitable to simulate events of longer duration; the calculation time is prohibitive. A restart from a highly accurate solution obtained at time t = T may then be an option, that is, restart with initial value η0 = η(z, T), where η(z, T)

has high spatial resolution. This procedure is to the detriment of some smearing, because interpolations are needed to find the restart's initial value for arbitrary z. The restart procedure is expected to work well for gradual valve closure, but for instantaneous valve closure (travelling contact discontinuities) the interpolation error might become unacceptable.

CONCLUSION

The presented method for simulating waterhammer with FSI is able to reveal subtle differences otherwise lost in numerical error, for example deviations due to marginally thicker sections, pipe clamps, small leaks, trapped air pockets, etc. These are minor effects that may accumulate in systems with low damping rates.

The exact solutions are to be used in the verification of numerical schemes.

REFERENCES

[1] Wiggert, D.C., 1996, "Fluid transients in flexible piping systems (a perspective on recent developments)", Proceedings of the 18th _{IAHR Symposium on Hydraulic Machinery and}

Cavitation, E. Cabrera, V. Espert, and F. Martínez, eds., Valencia, Spain; Kluwer Academic Publishers, Dordrecht, The Netherlands, ISBN 0-7923-4210-0, pp. 58-67.

[2] Tijsseling, A.S., 1996, "Fluid-structure interaction in liquid-filled pipe systems: a review", Journal of Fluids and Structures 10 (2) 109-146.

[3] Wiggert, D.C., and Tijsseling, A.S., 2001, "Fluid transients and fluid-structure interaction in flexible liquid-filled piping", ASME Applied Mechanics Reviews 54 (5) 455-481.

[4] Tijsseling, A.S., and Bergant, A., 2007, "Meshless computation of water hammer", Scientific Bulletin of the “Politehnica” University of Timişoara, Transactions on Mechanics 52(66) (6) 65-76, ISSN 1224-6077.

[5] Tijsseling, A.S., 2003, "Exact solution of linear hyperbolic four-equation systems in axial liquid-pipe vibration", Journal of Fluids and Structures 18 (2) 179-196.

[6] Li, Q.S., Yang, K., and Zhang, L., 2003, "Analytical solution for fluid-structure interaction in liquid-filled pipes subjected to impact-induced water hammer", ASCE Journal of Engineering Mechanics 129 (12) 1408-1417.

[7] Skalak, R., 1956, "An extension of the theory of waterhammer", Transactions of the ASME 78 105-116.

7 [8] Tijsseling, A.S., 1993, "Fluid-structure interaction in case of

waterhammer with cavitation", Ph.D. thesis, Delft University of Technology, Faculty of Civil Engineering, Communications on Hydraulic and Geotechnical Engineering, Report No. 93-6, ISSN 0169-6548, Delft, The Netherlands;

http://repository.tudelft.nl/file/264337/201268.

[9] Lavooij, C.S.W., and Tijsseling, A.S., 1991, "Fluid-structure interaction in liquid-filled piping systems", Journal of Fluids and Structures 5 (5) 573-595.

[10] Bürmann, W., 1975, "Water hammer in coaxial pipe systems", ASCE Journal of the Hydraulics Division 101 (6) 699-715. (Errata on p. 1554.)

[11] Williams, D.J., 1977, "Waterhammer in non-rigid pipes: precursor waves and mechanical damping", IMechE Journal of Mechanical Engineering Science 19 (6) 237-242.

[12] Wilkinson, D.H., and Curtis, E.M., 1980, "Water hammer in a thin-walled pipe", Proceedings of the BHRA Third International Conference on Pressure Surges, H.S. Stephens, and J.A. Hanson, eds., Canterbury, UK, pp. 221-240.

8 APPENDIX A

ALGORITHMS

Interior point calculation

The subroutine INTERIOR calculates η1 and η2 in the interior points of pipe 1 and 2, respectively, from the boundary and junction values. See Figs 1 and 2. Noting that λ12, λ14, λ22 and

λ24 are negative numbers herein, the pseudo-code reads:

INTERIOR (input: z, t ; output: η1, η2) if (0 < z < z j) then CALL BOUNDARY (0, t – z / λ11 ; η1, η2) η1 := η1₁ CALL BOUNDARY (z j, t – (z− z j) / λ12 ; η1, η2) η2 := η12 CALL BOUNDARY (0, t – z / λ13 ; η1, η2) η3 := η13 CALL BOUNDARY (z j, t – (z− z j) / λ14 ; η1, η2) η4 := η14 1 1 η := η1 2 1 η := η2 3 1 η := η3 4 1 η := η4 if (z j < z < L) then CALL BOUNDARY (z j, t – (z−z j) / λ21 ; η1, η2) η1 := η21 CALL BOUNDARY (L, t – (z−L) / λ22 ; η1, η2) η2 := η2₂ CALL BOUNDARY (z j, t – (z−z j) / λ23 ; η1, η2) η3 := η23 CALL BOUNDARY (L, t – (z−L) / λ24 ; η1, η2) η4 := η2₄ 1 2 η := η1 2 2 η := η2 3 2 η := η3 4 2 η := η4 else

"z is not an interior point" end

Boundary point calculation

The recursive subroutine BOUNDARY calculates η1 and η2 in the boundary and junction points (finally) from constant initial values η10 and η20. See Fig. 2. Noting that λ12, λ14, λ22 and λ24

are negative numbers herein, the pseudo-code reads: BOUNDARY (input: z, t ; output: η1, η2) if (t ≤ 0) then η 1 := η 10 η 2 := η 20 else if (z = 0) then CALL BOUNDARY (z j, t + z j / λ12 ; η1, η2) η2 := η12 CALL BOUNDARY (z j, t + z j / λ14 ; η1, η2) η4 := η1₄ 1 1 η := α12(t) η2 + α14(t) η4 + + β11(t) q1(t) + β13(t) q3(t) 2 1 η := η2 3 1 η := α32(t) η2 + α34(t) η4 + + β31(t) q1(t) + β33(t) q3(t) 4 1 η := η4 if (z = L) then CALL BOUNDARY (z j, t – (L–z j) / λ21 ; η1, η2) η1 := η21 CALL BOUNDARY (z j, t – (L–z j) / λ23 ; η1, η2) η3 := η2₃ 1 2 η := η1 2 2 η := α21(t) η1 + α23(t) η3 + + β22(t) q2(t) + β24(t) q4(t) 3 2 η := η3 4 2 η := α41(t) η1 + α43(t) η3 + + β42(t) q2(t) + β44(t) q4(t) if (z = z j) then CALL BOUNDARY (0, t – z j / λ11 ; η1, η2) η1 := η11 CALL BOUNDARY (L, t – (z j−L) / λ22 ; η1, η2) η2 := η2₂

9 CALL BOUNDARY (0, t – z j / λ13 ; η1, η2) η3 := η1₃ CALL BOUNDARY (L, t – (z j−L) / λ24 ; η1, η2) η4 := η24 1 1 η := η1 2 1 η := γ21(t) η1 + γ22(t) η2 + + γ23(t) η3 + γ24(t) η4 + δ2 3 1 η := η3 4 1 η := γ41(t) η1 + γ42(t) η2 + + γ43(t) η3 + γ44(t) η4 + δ4 1 2 η := γ11(t) η1 + γ12(t) η2 + + γ13(t) η3 + γ14(t) η4 + δ1 2 2 η := η2 3 2 η := γ31(t) η1 + γ32(t) η2 + + γ33(t) η3 + γ34(t) η4 + δ3 4 2 η := η4

if (z ≠ 0 AND z ≠ L AND z ≠ z j) then

"z is not at a boundary or junction" end

Coefficients

The boundary conditions in Table 3 are expressed in terms of general matrix-vector equations in accordance with Ref. [5] by:

0 1 0 0 1 0 1 0 0 0 1 0 0 A2_f 0 A2_s ⎛ ⎞ ⎜ ₋ ⎟ ⎜ ⎟ = ⎜ ⎟ ⎜ ⎟ ⎜ ₋ ⎟ ⎝ ⎠ D and q=0 . (21)

The coefficients α and β determine the unknowns η 1,η 3 and

η 2,η 4 at the left and right boundaries, respectively. The

coefficients γ and δ determine the unknowns η 2,η 4 and η 1,η 3 at

the left and right sides of the junction, respectively. The coefficients α, β, γ and δ are defined in terms of DS (D multiplied by S, see Eq. 16) and ₌ −1

JS JSU JSK (see Eq. 19) in Table 4. Implementation

The algorithms have been implemented in Mathcad 11 and 14 and can be found in Appendix B.

Table 4. Coefficients used in the subroutine BOUNDARY.

boundary z = 0 in pipe 1 boundary z = L in pipe 2 junction at z = z j

11 33 31 13 det13=DS1 DS1 −DS1 DS1 det 24=DS2 DS2_{22 44} −DS2 DS2 _{42 24}

(

12 33 13 32

)

det13 − − α12 = DS1 DS1 DS1 DS1

(

14 33 13 34

)

det13 − − α14 = DS1 DS1 DS1 DS1

(

32 11 31 12

)

det13 − − α32 = DS1 DS1 DS1 DS1

(

34 11 31 14

)

det13 − − α34 = DS1 DS1 DS1 DS1 33 det13 β11 = DS1 13 det13 β13 = −DS1 31 det13 β31 = −DS1 11 det13 β33 =DS1

(

21 44 24 41

)

det 24 − − α21 = DS2 DS2 DS2 DS2

(

23 44 24 43

)

det 24 − − α23 = DS2 DS2 DS2 DS2

(

41 22 42 21

)

det 24 − − α41 = DS2 DS2 DS2 DS2

(

43 22 42 23

)

det 24 − − α43 = DS2 DS2 DS2 DS2 44 det 24 β22 =DS2 det 24 β24 = −DS224 _{misprint in [5]} 42 det 24 β42 = −DS2 22 det 24 β44 =DS2 ij ij γ = JS 1 ( )_i i − δ = JSU r

Mathcad 14 CASA-09-14_Appendix B.xmcd B 1 .

APPENDIX B

Meshless Water Hammer with FSI

"Double" pipe, Poisson and junction coupling, no friction, instantaneous valve closure.

Input data:

modified Delft Hydraulics Benchmark Problem A (

Figure 3, Table 1

)

Liquid

K:= 2.1 10⋅ 9⋅Pa ρ_f 1000 kg m3 ⋅ := Q₀ 0.5 m 3 s ⋅ :=

Pipe 1

L1:= 10 m⋅ R1:= 398.5 mm⋅ ee1:= 16 mm⋅ E1:= 210 10⋅ 9⋅Pa ρ1_t 7900 kg m3 ⋅ := ν1:= 0.30

derived data:

A1_f:= π R1⋅ 2 A1_f= 0.499 m2 A1_t:= π⋅

⎡⎣

(R1+ee1)2−R12

⎤⎦

A1_t=0.041 m2 V1₀ Q₀ A1_f := V1₀= 1.002 m s⋅ −1 V1₀ 1 m s ⋅ := Q₀:= A1_f⋅V1₀ Q₀= 0.499 m3⋅s−1

Pipe 2

L2:= 10 m⋅ R2:= 398.5 mm⋅ ee2:= 8 mm⋅ E2:= 210 10⋅ 9⋅Pa ρ2_t 7900 kg m3 ⋅ := ν2:= 0.30

derived data:

A2_f:= π R2⋅ 2 A2_f= 0.499 m2 A2_t:= π⋅

⎡⎣

(R2+ee2)2−R22

⎤⎦

A2_t=0.02 m2 V2₀ Q₀ A2_f := V2₀= 1 m s⋅ −1 V2₀ 1 m s ⋅ := Q₀:= A2_f⋅V2₀ Q₀= 0.499 m3⋅s−1

Mathcad 14 CASA-09-14_Appendix B.xmcd B 2 .

Longitudinal wave speeds

(

Eqs 9-11, Table 2

)

Pipe 1

c12_f ρ_f 1 K 2 R1⋅ E1 ee1⋅ 1 ν1 2 −

(

)

⋅ +

⎡⎢

⎣

⎤⎥

⎦

⋅

⎡⎢

⎣

⎤⎥

⎦

1 − := c1_f:= c12_f c1_f= 1202.079 m s⋅ −1 c12_t E1 ρ1_t := c1_t:= c12_t c1_t= 5155.8 m s⋅ −1 γ12 c12_f+ c12_t 2 ν1⋅ 2 ρ_f ρ1_t ⋅ R1 ee1 ⋅ ⋅c12_f + := γ12= 2.885×107m2⋅s−2 λ12₁ 1 2 γ12 γ12 2 4 c12⋅ _f⋅c12_t − −

⎛

⎝

⎞

⎠

⋅ := λ1₁:= λ12₁ λ1₁= 1182.973 m s⋅ −1 λ1₂:= −λ1₁ λ1₂= −1182.973m s⋅ −1 λ12₃ 1 2 γ12 γ12 2 4 c12⋅ _f⋅c12_t − +

⎛

⎝

⎞

⎠

⋅ := λ1₃:= λ12₃ λ1₃= 5239.07 m s⋅ −1 λ1₄:= −λ1₃ λ1₄= −5239.07m s⋅ −1

Pipe 2

c22_f ρ_f 1 K 2 R2⋅ E2 ee2⋅ 1 ν2 2 −

(

)

⋅ +

⎡⎢

⎣

⎤⎥

⎦

⋅

⎡⎢

⎣

⎤⎥

⎦

1 − := c2_f:= c22_f c2_f= 1049.497 m s⋅ −1 c22_t E1 ρ1_t := c2_t:= c22_t c2_t= 5155.8 m s⋅ −1 γ22 c22_f+ c22_t 2 ν2⋅ 2 ρ_f ρ2_t ⋅ R2 ee2 ⋅ ⋅c22_f + := γ22= 2.893×107m2⋅s−2 λ22₁ 1 2 γ22 γ22 2 4 c22⋅ _f⋅c22_t − −

⎛

⎝

⎞

⎠

⋅ := λ2₁:= λ22₁ λ2₁= 1024.711 m s⋅ −1 λ2₂:= −λ2₁ λ2₂= −1024.711m s⋅ −1 λ22₃ 1 2 γ22 γ22 2 4 c22⋅ _f⋅c22_t − +

⎛

⎝

⎞

⎠

⋅ := λ2₃:= λ22₃ λ2₃= 5280.511 m s⋅ −1 λ2₄:= −λ2₃ λ2₄= −5280.511m s⋅ −1

Mathcad 14 CASA-09-14_Appendix B.xmcd B 3 .

Matrices of coefficients A and B (for water hammer with FSI) (

Eqs 5-8; Ref 5, Eqs 21-22

)

Note:

The matrices are made dimensionless through multiplication with the appropriate unit.

Scaling P/K and

σ

/E needed?

ORIGIN≡ 1

Pipe 1

A1

_{1 1,} := 1

A1

_{2 2,} 1 ρfc1f 2 ⋅ Pa ⋅ :=

A1

_{3 3,} := 1

A1

_{4 4,} 1 ρ1tc1t 2 ⋅ − ⋅Pa :=

A1

_{4 2,} ν1 E1 R1 ee1 ⋅ ⋅Pa :=

B1

_{2 1,} := 1

B1

_{1 2,} 1 ρf kg m3 ⋅ :=

B1

_{4 3,} := 1

B1

_{3 4,} 1 ρ1t − kg m3 ⋅ :=

B1

_{2 3,} := − ν12⋅

B1

0 1 0 0 0.001 0 0 0 0 0.6 − 0 1 0 0 1.266 − × 10−4 0

⎛

⎜

⎜

⎜

⎜

⎝

⎞

⎟

⎟

⎟

⎟

⎠

=

A1

1 0 0 0 0 6.92× 10−10 0 3.558×10−11 0 0 1 0 0 0 0 4.762 − ×10−12

⎛⎜

⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟

⎟⎠

=

Pipe 2

A2

_{1 1,} := 1

A2

_{2 2,} 1 ρfc2f 2 ⋅ Pa ⋅ :=

A2

_{3 3,} := 1

A2

_{4 4,} 1 ρ2tc2t 2 ⋅ − ⋅Pa :=

A2

_{4 2,} ν2 E2 R2 ee2 ⋅ ⋅Pa :=

B2

_{2 1,} := 1

B2

_{1 2,} 1 ρf kg m3 ⋅ :=

B2

_{4 3,} := 1

B2

_{3 4,} 1 ρ2t − kg m3 ⋅ :=

B2

_{2 3,} := − ν22⋅

B2

0 1 0 0 0.001 0 0 0 0 0.6 − 0 1 0 0 1.266 − × 10−4 0

⎛

⎜

⎜

⎜

⎜

⎝

⎞

⎟

⎟

⎟

⎟

⎠

=

A2

1 0 0 0 0 9.079×10−10 0 7.116×10−11 0 0 1 0 0 0 0 4.762 − ×10−12

⎛⎜

⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟

⎟⎠

=

Mathcad 14 CASA-09-14_Appendix B.xmcd B 4 .

Transformation matrix T (

Ref 5, Eq 27

)

Pipe 1

α1₁ c1t 2 λ1₁2 − λ1₁2 := α1₁= 17.995 t111:= 1 t112 λ11 sec m ⋅ := t113 2 ν1⋅ α1₁ := t114 c1_t2 λ1₁⋅t113 sec m ⋅ := t121:= t111 t122:= −t112 t123:= t113 t124:= −t114 α1₃ c1f 2 λ1₃2 − λ1₃2 := α1₃= −0.947 t131 ρf ν1 E1 ⋅ R1 ee1 ⋅ c1f 2 α1₃ ⋅ := t132 λ13⋅t131 sec m ⋅ := t133 λ1₃2 c1_t2 := t134 λ13 sec m ⋅ := t141:= t131 t142:= −t132 t143:= t133 t144:= −t134

T1

t111 t1₂₁ t1₃₁ t1₄₁ t112 t1₂₂ t1₃₂ t1₄₂ t113 t1₂₃ t1₃₃ t1₄₃ t114 t1₂₄ t1₃₄ t1₄₄

⎛

⎜

⎜

⎜

⎜

⎝

⎞

⎟

⎟

⎟

⎟

⎠

:=

T1

1 1 0.054 − 0.054 − 1182.973 1182.973 − 284.327 − 284.327 0.033 0.033 1.033 1.033 749.227 749.227 − 5239.07 5239.07 −

⎛⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟⎠

=

Pipe 2

α2₁ c2t 2 λ2₁2 − λ2₁2 := α2₁= 24.316 t211:= 1 t212 λ21 sec m ⋅ := t213 2 ν2⋅ α2₁ := t214 c2_t2 λ2₁⋅t213 sec m ⋅ := t221:= t211 t222:= −t212 t223:= t213 t224:= −t214 α2₃ c2f 2 λ2₃2 − λ2₃2 := α2₃= −0.96 t231 ρf ν2 E2 ⋅ R2 ee2 ⋅ c2f 2 α2₃ ⋅ := t232 λ23⋅t231 sec m ⋅ := t233 λ2₃2 c2_t2 := t234 λ23 sec m ⋅ := t241:= t231 t242:= −t232 t243:= t233 t244:= −t234

T2

t211 t2₂₁ t2₃₁ t2₄₁ t212 t2₂₂ t2₃₂ t2₄₂ t213 t2₂₃ t2₃₃ t2₄₃ t214 t2₂₄ t2₃₄ t2₄₄

⎛

⎜

⎜

⎜

⎜

⎝

⎞

⎟

⎟

⎟

⎟

⎠

:=

T2

1 1 0.082 − 0.082 − 1024.711 1024.711 − 430.905 − 430.905 0.025 0.025 1.049 1.049 640.112 640.112 − 5280.511 5280.511 −

⎛⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟⎠

=

Mathcad 14 CASA-09-14_Appendix B.xmcd B 5 .

Transformation matrix S (

see Ref 5

)

Pipe 1

S1

:=

(

T1 A1

⋅

)

−1

S1

0.499 5.905×105 0.026 2.452 − × 105 0.499 5.905 − × 105 0.026 2.452×105 0.016 − 8.444 − × 104 0.483 2.001 − × 107 0.016 − 8.444× 104 0.483 2.001× 107

⎛⎜

⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟

⎟⎠

=

S1

−1 1 1 0.054 − 0.054 − 8.453× 10−7 8.453 − ×10−7 1.036 − ×10−8 1.036× 10−8 0.033 0.033 1.033 1.033 3.568 − × 10−9 3.568× 10−9 2.495 − × 10−8 2.495× 10−8

⎛⎜

⎜

⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟

⎟

⎟⎠

=

T1 A1

⋅ 1 1 0.054 − 0.054 − 8.453×10−7 8.453 − × 10−7 1.036 − × 10−8 1.036×10−8 0.033 0.033 1.033 1.033 3.568 − × 10−9 3.568×10−9 2.495 − × 10−8 2.495×10−8

⎛⎜

⎜

⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟

⎟

⎟⎠

=

Pipe 2

S2

:=

(

T2 A2

⋅

)

−1

S2

0.499 5.114×105 0.039 3.143 − × 105 0.499 5.114 − × 105 0.039 3.143×105 0.012 − 6.199 − × 104 0.476 1.985 − × 107 0.012 − 6.199× 104 0.476 1.985× 107

⎛⎜

⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟

⎟⎠

=

S2

−1 1 1 0.082 − 0.082 − 9.759× 10−7 9.759 − ×10−7 1.545 − ×10−8 1.545× 10−8 0.025 0.025 1.049 1.049 3.048 − × 10−9 3.048× 10−9 2.515 − × 10−8 2.515× 10−8

⎛⎜

⎜

⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟

⎟

⎟⎠

=

T2 A2

⋅ 1 1 0.082 − 0.082 − 9.759×10−7 9.759 − × 10−7 1.545 − × 10−8 1.545×10−8 0.025 0.025 1.049 1.049 3.048 − × 10−9 3.048×10−9 2.515 − × 10−8 2.515×10−8

⎛⎜

⎜

⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟

⎟

⎟⎠

=

Mathcad 14 CASA-09-14_Appendix B.xmcd B 6 .

Check transformation matrix TA

TA1

1 1, := 1

TA1

1 2, λ1₁ ρ_f⋅c12f 2 ν1⋅ 2⋅R1 ρ1_t⋅ee1 λ1₁ c12t− λ121 ⋅ +

⎛

⎜

⎝

⎞

⎟

⎠

⋅Pa sec m ⋅ :=

TA1

1 3, 2 ν1⋅ λ12₁ c12t−λ121 ⋅ :=

TA1

1 4, 2 ν1⋅ ρ1_t − λ11 c12t−λ121 ⋅ ⋅Pa sec m ⋅ :=

TA1

3 1, ρf ν1 R1⋅ E1 ee1⋅ ⋅ c12f⋅λ123 c12f−λ123 ⋅ :=

TA1

3 2, ν1 R1⋅ E1 ee1⋅ c12f⋅λ13 c12f−λ123

⎛

⎜

⎝

⎞

⎟

⎠

⋅ ⋅Pa sec m ⋅ :=

TA1

3 3, λ12₃ c12t :=

TA1

3 4, 1 ρ1_t − λ13 c12t ⋅ ⋅Pa sec m ⋅ :=

TA1

2 1, :=

TA1

1 1,

TA1

2 2, := −

TA1

1 2,

TA1

2 3, :=

TA1

1 3,

TA1

2 4, := −

TA1

1 4,

TA1

4 1, :=

TA1

3 1,

TA1

4 2, := −

TA1

4 1,

TA1

4 3, :=

TA1

3 3,

TA1

4 4, := −

TA1

3 4,

TA1

1 1 0.054270454309537 − 0.054270454309537 − 8.453278589750525×10−7 8.453278589750525 − × 10−7 1.035879477358778 − × 10−8 0.054270454309537 0.033342330199931 0.033342330199931 1.032562272585722 1.032562272585722 3.567746911537571 − ×10−9 3.567746911537571× 10−9 2.494795393246406 − ×10−8 2.494795393246406× 10−8

⎛⎜

⎜

⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟

⎟

⎟⎠

=

T1 A1

⋅ 1 1 0.054270454309537 − 0.054270454309537 − 8.453278589750523× 10−7 8.453278589750523 − ×10−7 1.035879477358777 − ×10−8 1.035879477358777× 10−8 0.033342330199931 0.033342330199931 1.032562272585722 1.032562272585722 3.56774691153757 − ×10−9 3.56774691153757× 10−9 2.494795393246405 − × 10−8 2.494795393246405×10−8

⎛⎜

⎜

⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟

⎟

⎟⎠

=

Mathcad 14 CASA-09-14_Appendix B.xmcd B 7 .

Algebraic transformation matrix S1 (

Eqs 12-15

)

denom1 1 2 ν1⋅ 2 ρf⋅R1 ρ1t⋅ee1 ⋅ c12f c12f− λ123 ⋅ λ121 c12t− λ121 ⋅ − := denom1=1.00175243997939 denom2 1 2 ν1⋅ 2 ρf⋅R1 ρ1t⋅ee1 ⋅ c12f c12f− λ123 ⋅ λ123 c12t− λ121 ⋅ − := denom2=1.034371775993411

SS1

_{1 1,} 1 2 1 denom1 ⋅ :=

SS1

_{1 3,} −ν1 c12t c12t−λ121 ⋅ λ121 λ123 ⋅ 1 denom1 ⋅ :=

SS1

_{2 1,} ρf⋅c12f 2 λ1⋅ 1 1 denom2 ⋅ 1 Pa ⋅ m s ⋅ :=

SS1

_{2 3,} −ν1 c12t c12t−λ121 ⋅ ρf⋅c12f λ13 ⋅ 1 denom2 ⋅ 1 Pa ⋅ m s ⋅ :=

SS1

_{3 1,} −ν1 ρf⋅R1 2 ρ1⋅ t⋅ee1 ⋅ c12f c12f−λ123 ⋅ 1 denom1 ⋅ :=

SS1

_{3 3,} c12t 2 λ12⋅ 3 1 denom1 ⋅ :=

SS1

_{4 1,} ν1 R1 ee1 ⋅ c12f c12f− λ123 ⋅ ρf⋅c12f 2 λ1⋅ 1 ⋅ 1 denom2 ⋅ 1 Pa ⋅ m s ⋅ :=

SS1

_{4 3,} 1 2 ν1⋅ 2 ρf⋅R1 ρ1t⋅ee1 ⋅ c12f c12t− λ121 ⋅ +

⎛

⎜

⎝

⎞

⎟

⎠

− ρ1t⋅c12t 2 λ1⋅ 3 ⋅ 1 denom2 ⋅ 1 Pa ⋅ m s ⋅ :=

SS1

_{1 4,} :=

SS1

_{1 3,}

SS1

_{1 2,} :=

SS1

_{1 1,}

SS1

_{2 4,} := −

SS1

_{2 3,}

SS1

_{2 2,} := −

SS1

_{2 1,}

SS1

_{3 4,} :=

SS1

_{3 3,}

SS1

_{3 2,} :=

SS1

_{3 1,}

SS1

_{4 4,} := −

SS1

_{4 3,}

SS1

_{4 2,} := −

SS1

_{4 1,}

SS1

0.499125312847041 5.90451749043529× 105 0.026233534000585 2.451651349287133 − × 105 0.499125312847041 5.90451749043529 − × 105 0.026233534000585 2.451651349287133×105 0.016117188700315 − 8.443908505541904 − ×104 0.483385192446691 2.000666323972741 − ×107 0.016117188700315 − 8.443908505541904×104 0.483385192446691 2.000666323972741×107

⎛⎜

⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟

⎟⎠

=

S1

0.499125312847041 5.904517490435292× 105 0.026233534000585 2.45165134928713 − × 105 0.499125312847041 5.904517490435292 − ×105 0.026233534000585 2.451651349287129× 105 0.016117188700315 − 8.443908505541906 − ×104 0.483385192446691 2.000666323972741 − ×107 0.016117188700315 − 8.443908505541906× 104 0.483385192446691 2.000666323972741× 107

⎛⎜

⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟

⎟⎠

=

SS1

−

S1

0 1.164 − × 10−10 0 3.201 − × 10−10 0 1.164×10−10 0 3.783×10−10 0 1.455×10−11 0 3.725× 10−9 0 1.455 − ×10−11 0 3.725 − × 10−9

⎛⎜

⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟

⎟⎠

=

Mathcad 14 CASA-09-14_Appendix B.xmcd B 8 .

Algebraic transformation matrix S2 (

Eqs 12-15

)

den1 1 2 ν2⋅ 2 ρ_f⋅R2 ρ2_t⋅ee2 ⋅ c22f c22_f− λ22₃ ⋅ λ221 c22_t− λ22₁ ⋅ − := den1=1.00191960332544 den2 1 2 ν2⋅ 2 ρ_f⋅R2 ρ2_t⋅ee2 ⋅ c22f c22_f− λ22₃ ⋅ λ223 c22_t− λ22₁ ⋅ − := den2=1.050975381119511

SS2

_{1 1,} 1 2 1 den1 ⋅ :=

SS2

_{1 3,} −ν2 c22t c22_t−λ22₁ ⋅ λ221 λ22₃ ⋅ 1 den1 ⋅ :=

SS2

_{2 1,} ρf⋅c22f 2 λ2⋅ ₁ 1 den2 ⋅ 1 Pa ⋅ m s ⋅ :=

SS2

_{2 3,} −ν2 c22t c22_t−λ22₁ ⋅ ρf⋅c22f λ2₃ ⋅ 1 den2 ⋅ 1 Pa ⋅ m s ⋅ :=

SS2

_{3 1,} −ν2 ρf⋅R2 2 ρ2⋅ _t⋅ee2 ⋅ c22f c22_f−λ22₃ ⋅ 1 den1 ⋅ :=

SS2

_{3 3,} c22t 2 λ22⋅ ₃ 1 den1 ⋅ :=

SS2

_{4 1,} ν2 R2 ee2 ⋅ c22f c22_f− λ22₃ ⋅ ρf⋅c22f 2 λ2⋅ ₁ ⋅ 1 den2 ⋅ 1 Pa ⋅ m s ⋅ :=

SS2

_{4 3,} 1 2 ν2⋅ 2 ρ_f⋅R2 ρ2_t⋅ee2 ⋅ c22f c22_t− λ22₁ ⋅ +

⎛

⎜

⎝

⎞

⎟

⎠

− ρ2t⋅c22t 2 λ2⋅ ₃ ⋅ 1 den2 ⋅ 1 Pa ⋅ m s ⋅ :=

SS2

_{1 4,} :=

SS2

_{1 3,}

SS2

_{1 2,} :=

SS2

_{1 1,}

SS2

_{2 4,} := −

SS2

_{2 3,}

SS2

_{2 2,} := −

SS2

_{2 1,}

SS2

_{3 4,} :=

SS2

_{3 3,}

SS2

_{3 2,} :=

SS2

_{3 1,}

SS2

_{4 4,} := −

SS2

_{4 3,}

SS2

_{4 2,} := −

SS2

_{4 1,}

SS2

0.499042037245769 5.113739168000351×105 0.038822500394747 3.142765931134168 − × 105 0.499042037245769 5.113739168000351 − ×105 0.038822500394747 3.142765931134168× 105 0.011739310302031 − 6.198955497245001 − ×104 0.47574853700892 1.984634281749191 − ×107 0.011739310302031 − 6.198955497245001× 104 0.47574853700892 1.984634281749191× 107

⎛⎜

⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟

⎟⎠

=

S2

0.499042037245769 5.11373916800035×105 0.038822500394747 3.14276593113417 − × 105 0.499042037245768 5.11373916800035 − × 105 0.038822500394747 3.14276593113417 ×105 0.011739310302031 − 6.198955497244997 − ×104 0.47574853700892 1.984634281749191 − ×107 0.011739310302031 − 6.198955497244997×104 0.47574853700892 1.984634281749191×107

⎛⎜

⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟

⎟⎠

=

SS2

−

S2

0 5.821×10−11 0 2.328×10−10 0 5.821 − ×10−11 0 2.328 − ×10−10 0 3.638 − ×10−11 0 0 0 3.638×10−11 0 0

⎛⎜

⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟

⎟⎠

=

Mathcad 14 CASA-09-14_Appendix B.xmcd B 9 .

Constant initial conditions

Note:

Vectors are made dimensionless through multiplication with the appropriate unit.

Steady-state pressure loss over valve located at z = L = L1 + L2:

ξ₀:= 0 V2₀= 1 m s⋅ −1 ud2₀ 0 m sec ⋅ := V2_r0:= V2₀−ud2₀ P2₀ ξ₀ 1 2 ⋅ ρ⋅ V2_f⋅ _r0⋅ V2_r0 := σ2₀ A2_f A2_t⋅P20 :=

Pipe 2

ϕ2

_IC V20 sec m ⋅ P20 Pa ud2₀ sec m ⋅ σ2₀ Pa

⎛

⎜

⎜

⎜

⎜

⎜

⎜

⎜

⎜

⎜

⎝

⎞

⎟

⎟

⎟

⎟

⎟

⎟

⎟

⎟

⎟

⎠

:=

η2

_IC:=

S2

−1⋅

ϕ2

_IC

η2

_IC 1 1 0.082 − 0.082 −

⎛⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟⎠

=

ϕ2

_IC 1 0 0 0

⎛⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟⎠

=

Pipe 1

P1₀:= P2₀ σ1₀ A2_t⋅σ2₀+A1_f⋅P1₀− A2_f⋅P2₀ A1_t := ud1₀:= ud2₀

ϕ1

_IC V10 sec m ⋅ P10 Pa ud1₀ sec m ⋅ σ1₀ Pa

⎛

⎜

⎜

⎜

⎜

⎜

⎜

⎜

⎜

⎜

⎝

⎞

⎟

⎟

⎟

⎟

⎟

⎟

⎟

⎟

⎟

⎠

:=

η1

_IC:=

S1

−1⋅

ϕ1

_IC

η1

_IC 1 1 0.054 − 0.054 −

⎛⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟⎠

=

ϕ1

_IC 1 0 0 0

⎛⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟⎠

=

Mathcad 14 CASA-09-14_Appendix B.xmcd B 10 .

Boundary-condition matrices

(

Table 3, Eq 21

)

Note:

Matrices are made dimensionless through multiplication with the appropriate unit.

D

( )t 0 1 0 0 1 0 0 A2_f m2 0 1 − 1 0 0 0 0 A2_t − m2

⎛⎜

⎜

⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟

⎟

⎟⎠

:=

Right-hand-side (excitation or equilibrium) vector

(

Eq 21

)

Note:

Vectors are made dimensionless through multiplication with the appropriate unit.

q

( )t 0 0 0 0

⎛⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟⎠

:=

qc

:=

q

(1 sec⋅ )

qc

0 0 0 0

⎛⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟⎠

=

Matrix DS

Pipe 1

DS1

( )t :=

D

( )t ⋅

S1

DS1

(1.s) 5.905×105 0.473 0.026 2.995×105 5.905 − ×105 0.473 0.026 2.995 − ×105 8.444 − ×104 0.5 − 0.483 3.626× 105 8.444× 104 0.5 − 0.483 3.626 − ×105

⎛⎜

⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟

⎟⎠

=

Pipe 2

DS2

( )t :=

D

( )t ⋅

S2

DS2

(1.s) 5.114×105 0.46 0.039 2.615×105 5.114 − ×105 0.46 0.039 2.615 − ×105 6.199 − ×104 0.487 − 0.476 3.706× 105 6.199× 104 0.487 − 0.476 3.706 − ×105

⎛⎜

⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟

⎟⎠

=

Mathcad 14 CASA-09-14_Appendix B.xmcd B 11 .

Coefficients

α and β

Table 4

Pipe 1, z = 0

det13 t( ):=

DS1

( )t _{1 1,} ⋅

DS1

( )t _{3 3,} −

DS1

( )t _{3 1,} ⋅

DS1

( )t _{1 3,} det13 1 sec( ⋅ )=2.876×105

α12 t( ) −

(

DS1

( )t 1 2, ⋅

DS1

( )t 3 3, −

DS1

( )t 1 3, ⋅

DS1

( )t 3 2,

)

det13 t( ) := α12 1 sec( ⋅ )=0.985 α14 t( ) −

(

DS1

( )t 1 4, ⋅

DS1

( )t 3 3, −

DS1

( )t 1 3, ⋅

DS1

( )t 3 4,

)

det13 t( ) := α14 1 sec( ⋅ )=−0.284 α32 t( ) −

(

DS1

( )t 3 2, ⋅

DS1

( )t 1 1, −

DS1

( )t 3 1, ⋅

DS1

( )t 1 2,

)

det13 t( ) := α32 1 sec( ⋅ )=−0.108 α34 t( ) −

(

DS1

( )t 3 4, ⋅

DS1

( )t 1 1, −

DS1

( )t 3 1, ⋅

DS1

( )t 1 4,

)

det13 t( ) := α34 1 sec( ⋅ )=−0.985 β11 t( )

DS1

( )t 3 3, det13 t( ) := β11 1 sec( ⋅ )=1.681×10−6 β13 t( ) −

DS1

( )t 1 3, det13 t( ) := β13 1 sec( ⋅ )=0.294 β31 t( ) −

DS1

( )t 3 1, det13 t( ) := β31 1 sec( ⋅ )=−9.121× 10−8 β33 t( )

DS1

( )t 1 1, det13 t( ) := β33 1 sec( ⋅ )=2.053

Pipe 2, z = L

det24 t( ):=

DS2

( )t _{2 2,} ⋅

DS2

( )t _{4 4,} −

DS2

( )t _{4 2,} ⋅

DS2

( )t _{2 4,} det24 1 sec( ⋅ )=−2.98×105

α21 t( ) −

(

DS2

( )t 2 1, ⋅

DS2

( )t 4 4, −

DS2

( )t 2 4, ⋅

DS2

( )t 4 1,

)

det24 t( ) := α21 1 sec( ⋅ )=−0.145 α23 t( ) −

(

DS2

( )t 2 3, ⋅

DS2

( )t 4 4, −

DS2

( )t 2 4, ⋅

DS2

( )t 4 3,

)

det24 t( ) := α23 1 sec( ⋅ )=1.212 α41 t( ) −

(

DS2

( )t 4 1, ⋅

DS2

( )t 2 2, −

DS2

( )t 4 2, ⋅

DS2

( )t 2 1,

)

det24 t( ) := α41 1 sec( ⋅ )=0.808 α43 t( ) −

(

DS2

( )t 4 3, ⋅

DS2

( )t 2 2, −

DS2

( )t 4 2, ⋅

DS2

( )t 2 3,

)

det24 t( ) := α43 1 sec( ⋅ )=0.145 β22 t( )

DS2

( )t 4 4, det24 t( ) := β22 1 sec( ⋅ )=1.244 β24 t( ) −

DS2

( )t 2 4, det24 t( ) := β24 1 sec( ⋅ )=−1.636× 10−6 β42 t( ) −

DS2

( )t 4 2, det24 t( ) := β42 1 sec( ⋅ )=−0.877 β44 t( )

DS2

( )t 2 2, det24 t( ) := β44 1 sec( ⋅ )=−1.544× 10−6

Mathcad 14 CASA-09-14_Appendix B.xmcd B 12 .

Constant coefficients

α and β to speed up the calculation

Pipe 1, z = 0

α12c:= α12 1 sec( ⋅ ) α12c=0.985 β11c:= β11 1 sec( ⋅ ) β11c=1.681× 10−6 α14c:= α14 1 sec( ⋅ ) α14c=−0.284 β13c:= β13 1 sec( ⋅ ) β13c=0.294 α32c:= α32 1 sec( ⋅ ) α32c=−0.108 β31c:= β31 1 sec( ⋅ ) β31c=−9.121× 10−8 α34c:= α34 1 sec( ⋅ ) α34c=−0.985 β33c:= β33 1 sec( ⋅ ) β33c=2.053

Pipe 2, z = L

α21c:= α21 1 sec( ⋅ ) α21c=−0.145 β22c:= β22 1 sec( ⋅ ) β22c=1.244 α23c:= α23 1 sec( ⋅ ) α23c=1.212 β24c:= β24 1 sec( ⋅ ) β24c=−1.636× 10−6 α41c:= α41 1 sec( ⋅ ) α41c=0.808 β42c:= β42 1 sec( ⋅ ) β42c=−0.877 α43c:= α43 1 sec( ⋅ ) α43c=0.145 β44c:= β44 1 sec( ⋅ ) β44c=−1.544× 10−6

Mathcad 14 CASA-09-14_Appendix B.xmcd B 13 .

Junction-condition matrices and vector

(

Eq 18

)

Note:

Matrices and vector are made dimensionless through multiplication with the appropriate unit.

J1

( )t A1_f m2 0 0 0 0 1 0 A1_f m2 A1_f − m2 0 1 0 0 0 0 A1_t − m2

⎛⎜

⎜

⎜

⎜

⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟

⎟

⎟

⎟

⎟⎠

:=

J1

(1 s⋅ ) 0.499 0 0 0 0 1 0 0.499 0.499 − 0 1 0 0 0 0 0.041 −

⎛⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟⎠

=

J2

( )t A2_f m2 0 0 0 0 1 0 A2_f m2 A2_f − m2 0 1 0 0 0 0 A2_t − m2

⎛⎜

⎜

⎜

⎜

⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟

⎟

⎟

⎟

⎟⎠

:=

J2

(1 s⋅ ) 0.499 0 0 0 0 1 0 0.499 0.499 − 0 1 0 0 0 0 0.02 −

⎛⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟⎠

=

r

( )t 0 0 0 0

⎛⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟⎠

:=

Matrix JS

Junction side 1

JS1

( )t :=

J1

( )t ⋅

S1

JS1

(1.s) 0.236 5.905× 105 0.026 3.046× 105 0.236 5.905 − × 105 0.026 3.046 − × 105 0.249 − 8.444 − × 104 0.483 7.755×105 0.249 − 8.444×104 0.483 7.755 − × 105

⎛⎜

⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟

⎟⎠

=

Junction side 2

JS2

( )t :=

J2

( )t ⋅

S2

JS2

(1.s) 0.23 5.114× 105 0.039 2.615× 105 0.23 5.114 − × 105 0.039 2.615 − × 105 0.243 − 6.199 − × 104 0.476 3.706×105 0.243 − 6.199×104 0.476 3.706 − × 105

⎛⎜

⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟

⎟⎠

=

Mathcad 14 CASA-09-14_Appendix B.xmcd B 14 .

Matrices and vectors defining 4-equation junction system

(

Eq 20

)

JSU

( )t

JS2

( )t _{1 1,} −

JS2

( )t _{2 1,} −

JS2

( )t _{3 1,} −

JS2

( )t _{4 1,} −

JS1

( )t _{1 2,}

JS1

( )t _{2 2,}

JS1

( )t _{3 2,}

JS1

( )t _{4 2,}

JS2

( )t _{1 3,} −

JS2

( )t _{2 3,} −

JS2

( )t _{3 3,} −

JS2

( )t _{4 3,} −

JS1

( )t _{1 4,}

JS1

( )t _{2 4,}

JS1

( )t _{3 4,}

JS1

( )t _{4 4,}

⎛

⎜

⎜

⎜

⎜

⎜

⎝

⎞

⎟

⎟

⎟

⎟

⎟

⎠

:=

JSU

(1 s⋅ ) 0.23 − 5.114 − ×105 0.039 − 2.615 − ×105 0.236 5.905 − ×105 0.026 3.046 − ×105 0.243 6.199× 104 0.476 − 3.706 − ×105 0.249 − 8.444× 104 0.483 7.755 − ×105

⎛⎜

⎜

⎜

⎜

⎜⎝

=

JSK

( )t

JS1

( )t _{1 1,} −

JS1

( )t _{2 1,} −

JS1

( )t _{3 1,} −

JS1

( )t _{4 1,} −

JS2

( )t _{1 2,}

JS2

( )t _{2 2,}

JS2

( )t _{3 2,}

JS2

( )t _{4 2,}

JS1

( )t _{1 3,} −

JS1

( )t _{2 3,} −

JS1

( )t _{3 3,} −

JS1

( )t _{4 3,} −

JS2

( )t _{1 4,}

JS2

( )t _{2 4,}

JS2

( )t _{3 4,}

JS2

( )t _{4 4,}

⎛

⎜

⎜

⎜

⎜

⎜

⎝

⎞

⎟

⎟

⎟

⎟

⎟

⎠

:=

JSK

(1 s⋅ ) 0.236 − 5.905 − × 105 0.026 − 3.046 − × 105 0.23 5.114 − × 105 0.039 2.615 − × 105 0.249 8.444× 104 0.483 − 7.755 − × 105 0.243 − 6.199×104 0.476 3.706 − × 10

⎛⎜

⎜

⎜

⎜

⎜⎝

=

JS

( )t :=

JSU

( )t −1⋅

JSK

( )t

JS

(1 s⋅ ) 1.071 0.072 0.017 − 0.011 0.073 − 0.928 0.018 − 0.006 0.019 0.03 1.355 0.334 0.028 − 0.022 − 0.333 − 0.655

⎛⎜

⎜

⎜

⎜⎝

⎞⎟

⎟

⎟

⎟⎠

=

JS

(1 s⋅ ) = 1

Constant coefficients

γ and δ to speed up the calculation

Junction side 1, z =z

_j-

_{Junction side 2, z =z}

j+ γ21c:=

JS

(1 s⋅ )_{2 1,} γ21c=0.072 γ11c:=

JS

(1 s⋅ )_{1 1,} γ11c=1.071 γ22c:=

JS

(1 s⋅ )_{2 2,} γ22c=0.928 γ12c:=

JS

(1 s⋅ )_{1 2,} γ12c=−0.073 γ23c:=

JS

(1 s⋅ )_{2 3,} γ23c=0.03 γ13c:=

JS

(1 s⋅ )_{1 3,} γ13c=0.019 γ24c:=

JS

(1 s⋅ )_{2 4,} γ24c=−0.022 γ14c:=

JS

(1 s⋅ )_{1 4,} γ14c=−0.028 δ2

⎝

⎛

JSU

(1 s⋅ )−1⋅

r

(1 s⋅ )

⎞

⎠

2 := δ2= 0 δ1

⎛

⎝

JSU

(1 s⋅ )−1⋅

r

(1 s⋅ )

⎞

⎠

1 := δ1= 0 γ41c:=

JS

(1 s⋅ )_{4 1,} γ41c=0.011 γ31c:=

JS

(1 s⋅ )_{3 1,} γ31c=−0.017 γ42c:=

JS

(1 s⋅ )_{4 2,} γ42c=0.006 γ32c:=

JS

(1 s⋅ )_{3 2,} γ32c=−0.018 γ43c:=

JS

(1 s⋅ )_{4 3,} γ43c=0.334 γ33c:=

JS

(1 s⋅ )_{3 3,} γ33c=1.355 γ44c:=

JS

(1 s⋅ )_{4 4,} γ44c=0.655 γ34c:=

JS

(1 s⋅ )_{3 4,} γ34c=−0.333 δ4

⎝

⎛

JSU

(1 s⋅ )−1⋅

r

(1 s⋅ )

⎞

⎠

4 := δ4= 0 δ3

⎛

⎝

JSU

(1 s⋅ )−1⋅

r

(1 s⋅ )

⎞

⎠

3 := δ3= 0

Mathcad 14 CASA-09-14_Appendix B.xmcd B 15 .

Lengths

z_j:= L1 L:= L1 +L2 z_j=10 m L=20 m

Wave travel times, assuming that λ

₁

= -λ

₂

(

Table 2

)

Pipe 1

Δt1₁ L1 λ11 := Δt1₁=0.008453 s Δt1₂:= Δt1₁ Δt1₃ L1 λ13 := Δt1₃=0.001909 s Δt1₄:= Δt1₃

Pipe 2

Δt2₁ L2 λ21 := Δt2₁=0.009759 s Δt2₂:= Δt2₁ Δt2₃ L2 λ23 := Δt2₃=0.001894 s Δt2₄:= Δt2₃

Recursion "coast to coast"

App. A in CASA-09-14

ε:= 10−15⋅s εε:= 10−15

η

_{BOUNDARY z t}( , )

η

_{1 1,}

η1

_IC 1 ←

η

_{1 2,}

η1

_IC 2 ←

η

_{1 3,}

η1

_IC 3 ←

η

1 4,

η1

IC 4 ←

η

_{2 1,}

η2

_IC 1 ←

η

_{2 2,}

η2

_IC 2 ←

η

_{2 3,}

η2

_IC 3 ←

η

2 4,

η2

IC 4 ← t<ε if

η

←

η

_{BOUNDARY z}

(

_j, t−Δt1₄

)

η4←

η

_{1 4,}

η

←

η

_{BOUNDARY z}

(

j, t−Δt12

)

η2←

η

_{1 2,}

η

1 1, ←α12c η2⋅ +α14c η4⋅ +β11c⋅

q

( )t 1+ β13c⋅

q

( )t 3

η

_{1 2,} ←η2

η

1 3, ←α32c η2⋅ +α34c η4⋅ +β31c⋅

q

( )t 1+ β33c⋅

q

( )t 3

η

_{1 4,} ←η4 z=0 m⋅ if z=z_j if t≥ε if :=