Fluid-Structure Interaction in a Pipeline Embedded in Concrete During Water Hammer

1 Introduction

Water conveyance pipeline has been widely used in marine engineering, petrochemical engineering, energy and power engineering, spacecraft power system and also daily life. Repeated stimuli inside and outside the pipeline have been proven to excite pulsations in the pressurized fluid. Water hammer is a typical trigger frequently inducing structural vibration and new pressure fluctuations. This fluctuation further causes pipe expansion or contraction, which in turn affects the hydraulic pressure inside the pipe. This phenomenon is the typical fluid-structure interaction (FSI). However, FSI responses closely depend on the operating condition of the pipe system. Hence, FSI analysis should be carried out according to various constraints and boundary conditions (Mahmoodi, et al., 2019; Tijsseling, 2019; Rajbamshi, et al., 2020).

Friction coupling, Poisson coupling and Junction coupling are three main types when pipeline interacts with fluid. The first two emerge throughout the whole pipeline, whereas Junction coupling only happens in local positions including the elbows, branches, valves, boundaries, and variable cross-sections (Zanganeh, et al., 2015). System coupling in these locations becomes stronger (Alaei, et al., 2019). Furthermore, response to friction coupling is the weakest and has a long duration of oscillation (Huang and Alben, 2016). As the most important coupling form, Junction coupling depends on the robustness of the system and usually causes a pressure head larger than the classical water hammer (Karakouzian, et al., 2019). Poisson coupling and friction coupling can greatly influence Junction coupling. And in turn, Junction coupling has a greater impact on the pipeline system compared with Poisson coupling.

Regarding the constraints of pipeline, Li, et al. (2012) and Liu and Li. (2011) analyzed the elastic support and found the complex boundary conditions can be successfully simulated by the six-equation elastic spring model. Then, Li, et al. (2014) expanded the application of this elastic model. The FSI is continuously induced and transmitted in fluids and pipelines. This interaction becomes much stronger in pipes with few or no support (Liu and Li. 2011). Meanwhile, FSI has a crucial effect on fluid pressure and pipe stress when the pipe axial stiffness is greater than the rigidity of the supports (Li, et al., 2012). In a typical case study, Riedelmeier (Riedelmeier, et al., 2014) analyzed FSI phenomenon under four kinds of axial support and determined proper support to mitigate the displacement and stress in pipes.

The preceding achievements mainly focused on the FSI response of water conveyance systems with various pipeline layouts. Notably, the external constraints of pipes were generally discrete and their effects on FSI response were simplified as boundary conditions. For example, a pipe with continuity constraints degenerates into a multi-span pipeline model (Wu and Shih, 2001; Yang, et al., 2004). This simplification omits the global restriction of the constraint as well as its vibration. So the current model available in the literature is yet to be completely suitable. Herein, this paper investigates the constraints and dynamic responses in a water conveyance pipeline surrounded by concrete, focusing on fluid-pipe-concrete interaction. Pipe and concrete are integrated as a composite structure in this research, similar to that in a pre-stressed concrete cylinder pipe (PCCP) (Lee, et al., 2012; Hu, et al., 2019; Sun, et al., 2020). Compared to the typical FSI, FSI arising in such a system contains an extra pipe-concrete interaction. However, this extra effect is yet to be included in the current PCCP model.

A typical reservoir-pipe-valve system is used in this paper to produce water hammer. With the simplified elastic model of concrete, as well as continuity and motion equations of fluid, pipe and concrete, the fluid-pipe-concrete six-equation model is derived. Then an iterative approach is proposed to solve this model and the results get validated by experiment. Dynamic responses concerning the fluid-pipe-concrete interaction of the reservoir-pipe-valve system are subsequently scrutinized, including hydraulic pressure, pipe wall stress and vibration.

2 Numerical Model

When water conveyance system with a straight pipe suffers from water hammer, the axial structural vibration becomes obvious. In the fluid-pipe-concrete model, system characteristics can be described by solid constitutive equations and fluid governing equations. Then the FSI in such a composite structure is mathematically deduced. Two global hypotheses should be followed.

• The axial component of parameters in the control volume is equivalent to the average value in the corresponding section (Ting, et al., 2017).

• Axial stress, strain and displacement are separately assumed to be constant at the radial direction on a certain cross-section.

2.1 Stress at Inner Interface and Outer Interface

As for a piping system embedded in concrete, the pipe concrete coupling is taken into account. A sketch of stress from concrete imposing on pipe wall is shown in Figure 1. Figure 1A illustrates the external forces applied on an element of pipe wall. In this work, section of the concrete segment is considered to be annular. Then the inner interface is defined by the fluid and pipe wall, while outer interface by the pipe wall and ambient concrete, as seen in Figure 1B.

FIGURE 1

FIGURE 1. Sketch of stress acting on pipe wall. (A) Side view. (B) Cross-sectional view.

Since the pipeline is symmetric, the stress of concrete uniformly distributes around the pipe ring and can be assumed as a function of the radial coordinate (Wu, 2017; Mirjavadi, et al., 2020). By extension, Lame solutions of the radial stress of pipe wall σ_r, hoop stress of pipe wall σ_θ, radial stress of concrete σ_rs and hoop stress of concrete σ_θs in this encased system are separately presented as follows (Daoxiang and Weinlin, 2006):

${\sigma }_r=\frac{\frac{{\left(R+e\right)}^2}{r^2}-1}{\frac{{\left(R+e\right)}^2}{R^2}-1}p{|}_{r=R}-\frac{1-\frac{R^2}{r^2}}{1-\frac{R^2}{{\left(R+e\right)}^2}}{P}_b(1)$

${\sigma }_{\theta }=\frac{\frac{{\left(R+e\right)}^2}{r^2}+1}{\frac{{\left(R+e\right)}^2}{R^2}-1}p{|}_{r=R}-\frac{1+\frac{R^2}{r^2}}{1-\frac{R^2}{{\left(R+e\right)}^2}}{P}_b(2)$

${\sigma }_{rs}=\frac{\frac{R_s^2}{r^2}-1}{\frac{R_s^2}{{\left(R+e\right)}^2}-1}{P}_b,{\sigma }_{\theta s}=\frac{\frac{R_s^2}{r^2}+1}{\frac{R_s^2}{{\left(R+e\right)}^2}-1}{P}_b(3)$

where R is the internal diameter of pipe, e is the thickness of pipe wall, P_b and p|_{r = R} are external and internal pressure of pipe wall, separately.

Pipe wall and concrete have the identical radial strain ε_r|_{r = R+e} on their interface. Due to the assumption of linear radial displacement, ε_r|_{r = R+e} = ε_r|_{r = R}. According to Hooke law, ε_r|_{r = R+e} is defined as:

${\epsilon }_r{|}_{r=R+e}=\frac{1}{E}\left[{\sigma }_r{|}_{r=R+e}-\nu \left({\sigma }_{\theta }{|}_{r=R+e}+{\sigma }_z{|}_{r=R+e}\right)\right](4)$

where σ_z is the axial stress of pipe wall, E is the elastic modulus of pipe wall, and ν is the Poisson’s ratio of pipe wall.

According to the reference (Sinha, et al., 2001) in a concrete-encased piping system, the concrete strain is small enough so that the stress-strain relation degenerates into an elastic model, shown as:

$P_b={\sigma }_{rs}{|}_{r=R+e}=E_s{\epsilon }_{rs}{|}_{r=R+e}\mathrm{+}{\nu }_s\left({\sigma }_{\theta s}{|}_{r=R+e}+{\sigma }_{zs}{|}_{r=R+e}\right)(5)$

where σ_zs is the concrete axial stress because of expansion or contraction of the outer pipe wall. P_b is the constant pressure imposed on the outer pipe wall by ambient concrete, thus, P_b = σ_rs|_{r = R+e}. ν_s and E_s are the Poisson’s ratio and elastic modulus of concrete, separately.

Recalling linear displacement at the radial direction and substituting Eqs 1–4 into Eq. 5, yields the radial stress on the outer interface:

$P_b=\frac{E_s}{E}\left[-P_b-\nu \left(\frac{2R^{2}p{|}_{r=R}}{{\left(R+e\right)}^2-R^2}-\frac{{\left(R+e\right)}^2+R^2}{{\left(R+e\right)}^2-R^2}{P}_b+{\sigma }_z\right)\right]+{\nu }_s{\sigma }_{zs}+\frac{R_s^2+{\left(R+e\right)}^2}{R_s^2-{\left(R+e\right)}^2}{\nu }_s{P}_b(6)$

The fluid radial velocity component V_r is assumed to be constant at a certain section. So V_r is equal to the radial velocity of pipe u_r on the inner interface (Lai, et al., 2020). Then pressure relationship with respect to pipe center-line and inner interface is established:

$P-p{|}_{r=R}=\frac{1}{2}{\rho }_{f}R\frac{\partial u_r{|}_{r=R}}{\partial t}(7)$

where P is the pressure at the pipe center line. ρ_f is the fluid density. And radial motion equation of pipeline is defined as:

${\rho }_t\frac{\partial u_r}{\partial t}=\frac{Rp{|}_{r=R}}{\left(R+e/2\right)e}-\frac{\left(R+e\right)P_b}{\left(R+e/2\right)e}-\frac{{\displaystyle {\int }_R^{R+e}{\sigma }_{r}dr}}{\left(R+e/2\right)e}(8)$

where ρ_t is the density of pipe wall material. Substitute Eq. 1 into Eq. 8:

${\rho }_t\frac{\partial u_r}{\partial t}=\frac{R^{2}p{|}_{r=R}}{{\left(R+e/2\right)}^{2}e}-\frac{{\left(R+e\right)}^2{P}_b}{{\left(R+e/2\right)}^{2}e}(9)$

The following can be obtained by solving Eqs 6, 7, Eq. 9.

$P_b=aP+b{\sigma }_z+c{\sigma }_{zs}(10)$

$p{|}_{r=R}=dP+f{\sigma }_z+h{\sigma }_{zs}(11)$

The definition of coefficients a, b, c, d, f and h can be seen in the Appendix.

2.2 Governing Equations for Fluid

2.2.1 Axial Motion Equation of Incompressible Flow

The equation of axial motion of fluid at pipe centerline is:

${\rho }_f{A}_{f}dz\frac{DV}{Dt}+A_f\left(P+\frac{\partial P}{\partial z}-P\right)dz=0(12)$

where t and z are time and axial direction respectively. V is the fluid velocity at pipe centerline. A_f is the area of the flow section. The total acceleration of fluid is composed of the convective acceleration and the local one. So Eq. 12 can be converted to:

${\rho }_f{A}_{f}dz\left(\frac{\partial V}{\partial t}+V\frac{\partial V}{\partial z}\right)+A_f\left(P+\frac{\partial P}{\partial z}-P\right)dz=0(13)$

The one-dimensional axial motion equation of fluid at pipe centerline is obtained.

$\frac{\partial V}{\partial t}+V\frac{\partial V}{\partial z}+\frac{1}{{\rho }_f}\frac{\partial P}{\partial z}=0(14)$

2.2.2 Continuity Equation of Incompressible Flow

Ignoring the hoop velocity of the pipe, continuity equation of the symmetric flow can be expressed in cylindrical coordinate:

$\frac{\partial {\rho }_f}{\partial t}+\frac{1}{r}\frac{\partial }{\partial r}\left({\rho }_{f}r{V}_r\right)+\frac{\partial }{\partial z}\left({\rho }_{f}V\right)=0(15)$

Recalling the definition of fluid bulk modulus K, ∂P/∂ρ_f = K/ρ_f, Eq. 15 is converted to:

$\frac{1}{K}\frac{\partial P}{\partial t}+\frac{1}{r}\frac{\partial \left(r{V}_r\right)}{\partial r}+\frac{\partial V}{\partial z}=0(16)$

Averaging over the cross-section, a continuity equation of one-dimensional flow is derived as (Tijsseling, 2007):

$\frac{1}{\pi R^2}{\displaystyle {\int }_0^{R}2\pi r\left(\frac{1}{K}\frac{\partial P}{\partial t}+\frac{1}{r}\frac{\partial \left(r{V}_r\right)}{\partial r}+\frac{\partial V}{\partial z}\right)dr=0}(17)$

The Resultant Equation Is:

$\frac{1}{K}\frac{\partial P}{\partial t}+\frac{2V_r}{R}+\frac{\partial V}{\partial z}=\text{0}(18)$

$\frac{V_r}{R}=\frac{u_r}{R}=\frac{1}{E}\frac{\partial \left[{\sigma }_{\theta }{\mathrm{|}}_{r=R}-\nu \left({\sigma }_z{\mathrm{|}}_{r=R}+{\sigma }_r{\mathrm{|}}_{r=R}\right)\right]}{\partial t}(19)$

Substituting Eqs 1, 2 into Eq. 19:

$\frac{V_r}{R}=\frac{1}{E}\left[\frac{{\left(R+e\right)}^2+R^2}{{\left(R+e\right)}^2-R^2}\frac{\partial p{|}_{r=R}}{\partial t}-\frac{2{\left(R+e\right)}^2}{{\left(R+e\right)}^2-R^2}\frac{\partial P_b}{\partial t}-\nu \left(\frac{\partial p{|}_{r=R}}{\partial t}+\frac{\partial {\sigma }_z}{\partial t}\right)\right](20)$

Substituting Eqs 10, 11, into Eq. 20:

$\frac{V_r}{R}=\alpha \frac{\partial P}{\partial t}+\beta \frac{\partial {\sigma }_z}{\partial t}+\gamma \frac{\partial {\sigma }_{zs}}{\partial t}(21)$

in which

$\left(\begin{array}{c}\alpha \\ \beta \\ \gamma \end{array}\right)=-\frac{1}{E}\frac{2{\left(R+e\right)}^2}{{\left(R+e\right)}^2-R^2}\left(\begin{array}{c}a\\ b\\ c\end{array}\right)+\frac{1}{E}\left[\frac{{\left(R+e\right)}^2+R^2}{{\left(R+e\right)}^2-R^2}-\nu \right]\left(\begin{array}{c}d\\ f\\ h\end{array}\right)+\frac{\nu }{E}\left(\begin{array}{c}0\\ 1\\ 0\end{array}\right)$

Substituting Eq. 21 Into Eq. 18:

$\frac{\partial V}{\partial z}+\left(\frac{1}{K}+2\alpha \right)\frac{\partial P}{\partial t}+2\beta \frac{\partial {\sigma }_z}{\partial t}+2\gamma \frac{\partial {\sigma }_{zs}}{\partial t}=0(22)$

2.3 Governing Equations for Pipe

2.3.1 Axial Motion Equation of Structures

For pipe walls, the axial external force is the shear force. Considering pipe-concrete coupling, the motion equation of the pipe is:

${\mathit{\rho }}_{\mathit{t}}{A}_{t}dz\frac{\partial {\mathit{u}}_{\mathit{z}}}{\partial \mathit{t}}+A_t\left[\left({\mathit{\sigma }}_{\mathit{z}}+\frac{\partial {\mathit{\sigma }}_{\mathit{z}}}{\partial \mathit{z}}dz\right)-{\mathit{\sigma }}_{\mathit{z}}\right]=-G{\epsilon }_{s}2\pi \left(R+e\right)dz(23)$

where u_z and σ_z are axial velocity and axial stress of pipe, A_t is the cross-sectional area of the pipe wall. G = E/2(1+ν_s) is the shear modulus of elasticity of concrete. The axial strain of concrete ε_s = ∂w_s/∂z. w_s is the axial displacement of concrete.

Eq. 23 is then Simplified:

${\mathit{\rho }}_{\mathit{t}}\frac{\partial {\mathit{u}}_{\mathit{z}}}{\partial \mathit{t}}+\frac{\partial {\mathit{\sigma }}_{\mathit{z}}}{\partial \mathit{z}}=-\frac{R+e}{\left(R+e/2\right)e}\frac{E_s}{2\left(1+{\nu }_s\right)}\frac{\partial W_s}{\partial z}(24)$

Due to the concrete axial velocity u_s = ∂w_s/∂t, Eq. 24 is updated.

${\mathit{\rho }}_{\mathit{t}}\frac{\partial {\mathit{u}}_{\mathit{z}}}{\partial \mathit{t}}+\frac{\partial {\mathit{\sigma }}_{\mathit{z}}}{\partial \mathit{z}}+\frac{R+e}{\left(R+e/2\right)e}\frac{E_s}{2\left(1+{\nu }_s\right)}\frac{\partial \left({\displaystyle {\int }_t{u}_{s}dt}\right)}{\partial \mathit{z}}=0(25)$

2.3.2 Pipeline Continuity Equation

According to the Hooke law, the strain-stress relation applied to pipe axial direction is:

${\epsilon }_z=\frac{\partial W_z}{\partial z}=\frac{1}{E}\left[{\sigma }_z-\nu \left(\overline{{\sigma }_r}+\overline{{\sigma }_{\theta }}\right)\right](26)$

where w_z is the axial displacement of pipe wall and u_z = ∂w_z/∂t. According to Eqs 1, 2, the sum of radial and hoop mean stresses is:

$\overline{{\sigma }_r}+\overline{{\sigma }_{\theta }}=\frac{1}{2\pi \left(R+e/2\right)e}{\displaystyle \underset{R}{\overset{R+e}{\int }}2\pi \left({\sigma }_r+{\sigma }_{\theta }\right)}dr=\frac{R^2{\left(R+e\right)}^2}{2{\left(R+e/2\right)}^2{e}^2}\ln \left(1+\frac{e}{R}\right)p{|}_{r=R}-\frac{R^2}{2\left(R+e/2\right)e}{P}_b(27)$

For convenience, Eq. 27 is written as:

$\overline{{\sigma }_r}+\overline{{\sigma }_{\theta }}=kP+l{\sigma }_z+m{\sigma }_{zs}(28)$

in which

$\left(\begin{array}{c}k\\ l\\ m\end{array}\right)=\frac{R^2}{2\left(R+e/2\right)e}\left(\begin{array}{c}a\\ b\\ c\end{array}\right)+\frac{R^2{(R+e)}^2}{2{(R+e/2)}^2{e}^2}\ln \left(1+\frac{e}{R}\right)\left(\begin{array}{c}d\\ f\\ h\end{array}\right)$

Substituting Eq. 28 into Eq. 26:

$\frac{\partial u_z}{\partial z}-\left(\frac{1}{E}-\frac{\nu l}{E}\right)\frac{\partial {\sigma }_z}{\partial t}+\frac{\nu k}{E}\frac{\partial P}{\partial t}+\frac{\nu m}{E}\frac{\partial {\sigma }_{zs}}{\partial t}=0(29)$

2.4 Governing Equations for Concrete

Considering structural vibration caused by FSI, the axial motion equation of concrete is:

${\rho }_s{A}_{s}dz\frac{\partial u_s}{\partial t}+A_s\frac{\partial {\sigma }_s}{\partial z}dz=\frac{E_s{\epsilon }_s}{2\left(1+{\nu }_s\right)}2\pi \left(R+e\right)dz(30)$

where ρ_s and A_s are density and sectional area of concrete, respectively. Eq. 30 can be simplified.

${\rho }_s\frac{\partial u_s}{\partial t}+\frac{\partial {\sigma }_s}{\partial z}=\frac{\pi \left(R+e\right)E_s{\epsilon }_s}{\left(1+{\nu }_s\right)A_s}(31)$

According to Hooke law, the strain-stress relation of concrete is:

${\epsilon }_s=\frac{\partial w_s}{\partial z}=\frac{1}{E}\left[{\sigma }_s-\nu \left({\overline{\sigma }}_{sr}+{\overline{\sigma }}_{s\theta }\right)\right](32)$

Taking time derivative of Eq. 32:

$\frac{\partial u_s}{\partial z}=\frac{1}{E}\left[\frac{\partial {\sigma }_s}{\partial t}-\nu \frac{\partial \left({\overline{\sigma }}_{sr}+{\overline{\sigma }}_{s\theta }\right)}{\partial t}\right](33)$

Based on Eq. 3, the sum of radial and hoop mean stresses can be expanded as:

$\begin{array}{l}{\overline{\sigma }}_{sr}+{\overline{\sigma }}_{s\theta }=\frac{1}{2\pi \left(R_s+e_s/2\right)e_s}{\displaystyle \underset{R_s}{\overset{R_s+e_s}{\int }}2\pi r\left({\sigma }_{sr}+{\sigma }_{s\theta }\right)}dr\\ =\frac{R_s^2}{\left(R_s+e_s/2\right)e_s}{P}_b=\frac{R_s^2}{\left(R_s+e_s/2\right)e_s}\left(aP+b{\sigma }_z+c{\sigma }_{zs}\right)\end{array}(34)$

Then substituting Eq. 34 into Eq. 33:

$\frac{\partial u_s}{\partial z}=\frac{1}{E}\frac{\partial {\sigma }_s}{\partial t}-\frac{v_s}{E}\frac{R_s^2}{\left(R_s+e_s/2\right)e_s}\left(a\frac{\partial P}{\partial t}+b\frac{\partial {\sigma }_z}{\partial t}+c\frac{\partial {\sigma }_{zs}}{\partial t}\right)(35)$

Till now, a six-equation model describing the fluid-pipe-concrete system is derived. This model consists of continuity and motion equations of fluid and pipe, and constitutive equations of concrete. And these equations are coupled by the boundary conditions on their interfaces, at r = R and r = R + e.

3 Solution Methodology

3.1 Finite Volume Discretization

The finite volume method (FVM) can provide satisfactory predictions to flow transients by solving equations in each control volume (Ferras, et al., 2017) and by presenting each physical term in integral form (Cardiff, et al., 2016). Based on cell-centered FV discretization, an iterative approach is proposed in this paper to calculate the FSI responses. As shown in Figure 2A, the control volume includes two dimensions, namely the time step ΔT and the control volume length ΔZ. The model is solved at the center of each cell. In order to numerically keep higher accuracy, the discretization grids system is shown in Figure 2B, where ΔZ = 2Δz and ΔT = 2Δt.

FIGURE 2

FIGURE 2. Discretization grids system (A) Cell-centered form of FVM. (B) An equivalent grid system of A.

This method discretizes Eq. 14 and Eq. 22, Eq. 25, Eq. 29, Eq. 31 and Eq. 35 using the control volume and integrates differential equations from t to t + Δt in each control volume. The generalized differential equation can be expressed as:

$\mathit{A}\frac{\partial \mathit{Q}}{\partial t}\mathrm{+}\mathit{B}\frac{\partial \mathit{Q}}{\partial z}=\mathit{S}(36)$

The Crank-Nicolson implicit form of time-centered difference provides a second-order accuracy to predict the partial differential terms. And unconditional stability is guaranteed with respect to the solution process (Ramírez, et al., 2018).

${\displaystyle \underset{\mathrm{\Delta }t\mathrm{\Delta }z}{\iint }\left(\mathit{{\rm A}}\frac{\partial \mathit{Q}}{\partial t}+\mathit{B}\frac{\partial \mathit{Q}}{\partial z}-\mathit{S}\right)}dtdz=0(37)$

$\mathit{A}{\displaystyle \underset{\mathrm{\Delta }z}{\int }\partial \mathit{Q}dz+\mathit{B}{\displaystyle \underset{\mathrm{\Delta }t}{\int }\partial \mathit{Q}dt-\mathit{S}}}{\displaystyle \underset{\mathrm{\Delta }t\mathrm{\Delta }z}{\iint }dtdz}=0(38)$

$\mathit{A}\left({\mathit{Q}}_i^{n+1}-{\mathit{Q}}_i^n\right)\mathrm{\Delta }z+\frac{\mathit{B}}{2}\left[\left({\mathit{Q}}_{i+1}^{n+1}-{\mathit{Q}}_i^{n+1}\right)+\left({\mathit{Q}}_{i+1}^n-{\mathit{Q}}_i^n\right)\right]\mathrm{\Delta }t=\mathit{S}\mathrm{\Delta }z\mathrm{\Delta }t(39)$

$\mathit{A}{\mathit{Q}}_i^{n+1}\mathrm{\Delta }z+\frac{\mathit{B}}{2}\left({\mathit{Q}}_{i+1}^{n+1}-{\mathit{Q}}_i^{n+1}\right)\mathrm{\Delta }t=A{\mathit{Q}}_i^n\mathrm{\Delta }z-\frac{\mathit{B}}{2}\left({\mathit{Q}}_{i+1}^n-{\mathit{Q}}_i^n\right)\mathrm{\Delta }t+\mathit{S}\mathrm{\Delta }z\mathrm{\Delta }t(40)$

where Q denotes the state matrix in each control volume at each time step.

$\mathit{Q}\mathrm{=}{\left[\begin{array}{cccccc}V& u_z& u_s& P& {\sigma }_z& {\sigma }_s\end{array}\right]}^{\mathrm{T}}$

A, B and S are shown in the Appendix. n denotes the nth time step. i is the ith calculation nodes. Δt denotes the time step. Δz denotes the axial interval of each node.

According to Eqs 37–40, the iterative matrix between adjacent time steps for the water conveyance system is built as follows:

${\left(\begin{array}{cccccc}{\mathit{A}}_1& {\mathit{D}}_1& & & & \\ {\mathit{B}}_2& {\mathit{A}}_2& {\mathit{D}}_2& & & \\ & {\mathit{B}}_3& {\mathit{A}}_3& {\mathit{D}}_3& & \\ & & \ddots & \ddots & \ddots & \\ & & & {\mathit{B}}_{N-1}& {\mathit{A}}_{N-1}& {\mathit{D}}_{N-1}\\ & & & & {\mathit{B}}_N& {\mathit{A}}_N\end{array}\right)}_{8N\times 8N}{\left(\begin{array}{c}{\mathit{Q}}_1\\ {\mathit{Q}}_2\\ {\mathit{Q}}_3\\ ⋮\\ {\mathit{Q}}_{N-1}\\ {\mathit{Q}}_N\end{array}\right)}_{8N\times 1}={\left(\begin{array}{c}{\mathit{C}}_1\\ {\mathit{C}}_2\\ {\mathit{C}}_3\\ ⋮\\ {\mathit{C}}_{N-1}\\ {\mathit{C}}_N\end{array}\right)}_{8N\times 1}(41)$

in which

${\mathit{B}}_i=\frac{1}{2}\mathit{B}\left(\begin{array}{cc}{\mathit{I}}_{4\times 4}& \\ & 0_{4\times 4}\end{array}\right)$

${\mathit{D}}_i=\frac{1}{2}\mathit{B}\left(\begin{array}{cc}{0}_{4\times 4}& \\ & {\mathit{I}}_{4\times 4}\end{array}\right)$

${\mathit{C}}_i=\mathit{A}{\mathit{Q}}_i^n\mathrm{\Delta }z-\frac{\mathit{B}}{2}\left({\mathit{Q}}_{i+1}^n-{\mathit{Q}}_i^n\right)\mathrm{\Delta }t+\mathit{S}\mathrm{\Delta }z\mathrm{\Delta }t$

Vector C_i and Q_i separately denote the state vector in the ith node at the nth and (n+1)^th time step. i = 1 to N, where 1 and N refer to the inlet and outlet of the pipeline, respectively. When discontinuity happens between adjacent pipelines or fluid, this discontinuous section will be modeled as an interior boundary condition. A_i and C_i are included in the following equations of boundary conditions.

3.2 Boundary Conditions

Given that the cross-sectional area of the reservoir is much larger than that of the pipeline, the water level in the tank can be treated as a constant (Vardy et al., 1996):

$P=P_0(42)$

where P₀ denotes the pressure provided by the water level in the reservoir. The pipe inlet and outlet are rigidly connected to the system respectively, so the corresponding continuity equations of axial direction are as follows (De Santis and Shams, 2019):

$u_z=0,u_s=0(43)$

When the outlet valve is closed, fluid flow velocity, axial velocity of the pipe wall and discharge satisfy:

$A_f\left(V-u_z\right)=\tau Q_0(44)$

where τ = (1—t/T_c)^1.5, τ is the closing coefficient which gradually changes from 1 to 0 when the valve is closed. T_c is the total shutdown time. Subscript 0 denotes the initial value.

4 Method Validation

4.1 Experimental Data Comparison

A more general apparatus proposed by Vardy and Fan is used for benchmark purposes, as seen in Figure 3 (Vardy et al., 1996). This system consists of a single pipe of 4.5 m and a rod of 5.006 m. The pipe is suspended and kept horizontal by long wires. Both ends of it are sealed by steel caps. And the pipe contains pressurized water to prevent cavitation. Then the steel rod axially impacts the stationary pipe to generate waves in the fluid and in the pipe.

FIGURE 3

FIGURE 3. Experimental apparatus (Vardy et al., 1996).

This system during transients is solved by the six-equation model and the iterative methodology. Herein, the pipe is regarded as a thin-wall pipe encased in the same material (ρ_t = ρ_s, ν = ν_s, E = E_s). The corresponding cross-section is similar to the sketch in Figure 1B.

Results from the experiment and numerical iteration are presented in Figure 4, including pressure fluctuations at the impact end (Figure 4A) and at the remote end (Figure 4B) of the pipe. The time step in the numerical calculation is 10⁻⁴ s to ensure this algorithm can capture the high-frequency fluctuation. While in the test, the max sampling frequency is 80 Hz. So the black curve in Figure 4 is relatively smooth compared to the red line. The latter one contains more fluctuating features. In Figure 4A, pressure values at the first peak are 2.82 MPa for experiment and 2.74 MPa for simulation. The relative error is 2.84%. Moreover, this six-equation model mainly takes axial coupling and radial coupling into account, while failing to include torsion coupling. Due to experimental error, the discrepancy between the two curves is acceptable. Generally, numerical results are in good agreement with the experimental data. This reveals this numerical model is suitable to investigate the dynamics and kinematics of the piping system embedded in concrete. And the FV discretization method is feasible to solve this system.

FIGURE 4

FIGURE 4. Time history of pressure (A) Impact end. (B) Remote end.

4.2 Classical Model Comparison

A typical reservoir-pipe-valve system is built up with pipe length L = 10 m, R = 0.1 m, e = 0.01 m. Figure 5 sketches this system encased in concrete.

FIGURE 5

FIGURE 5. Side view of water conveyance system encased in concrete.

In this system, the water hammer induced by closing the valve was separately calculated by the classical two-equation model and the six-equation model established in this paper. As seen in Figure 6, little deviation arises between these two models when E_S = 10¹⁵ Pa. The two-equation model neglects fluid-pipe interaction based on the classical water hammer theory. With respect to the six-equation model, the elastic modulus of concrete is large enough so that the pipe vibration is negligible and the FSI response is weak. In such a case, the six-equation model degenerates to the classical model and their results are almost identical.

FIGURE 6

FIGURE 6. Dynamic pressure with and without FSI.

To further justify the derived six-equation model, comparisons to the two-equation model, and four-equation model are discussed with respect to the water hammer in a typical reservoir-pipe-valve system. The pressure fluctuations upstream the valve by these models are illustrated in Figure 7. The two-equation model omits the FSI, so the pressure curve is trapezoidal waveform in shape. The four-equation model involves fluid-pipe interaction but neglects constraints on the pipe’s outer surface. This is equivalent to an exposed pipe system.

FIGURE 7

FIGURE 7. Hydraulic pressure by different models.

The peak segment of the black line in Figure 7 can be split into three parts. In region I, the red and blue lines are lower than the black line due to pipe expansion at the valve. In region II, red and blue lines climb over the black line. Because in FSI-involved models, contraction of the pipe wall plays pumping effects on fluid, triggering a larger pressure than the classical result. In region III, the axial stress wave bounces back, expands again the pipe and results in the last pressure drop over the pressure surge. Moreover, an exposed pipe (four-equation model) has a smaller pressure wave speed, but stronger FSI responses due to less support. When the pipe is embedded into concrete, constraints weaken the contributions of the pumping effect. Peak pressure is suppressed, shown as the red line lower than the blue one. When the stiffness of the concrete is large enough, peak pressure and pressure wave speed are similar to those in the classical model, as shown in Figure 6.

5 Results and Discussions

This section deals with the kinematic and dynamic characteristics based on the layout in Figure 5 when the fluid-pipe-concrete interaction is taken into account. Hydraulic transients is triggered by closing the valve in 0.01 s. And the closing law depends on the coefficient τ in Eq. 44.

Sensitive analysis of grid size and time step is first discussed, which determines the result precision. In this process, the Courant number is well controlled to ensure a stable iteration and reliable results. Numerical results of peak pressure vary with these two parameters are given in Figure 8. When Δz > 0.1 m in Figure 8A, the peak pressure drops dramatically with the increase of grid interval. On the contrary, if Δz < 0.1 m, the numerical result is not sensitive to grid size. Figure 8B shows a similar manner when the threshold value of the time step is 10⁻⁴ s. Consequently, the node interval and time step are separately set to be 0.1 m and 10⁻⁴ s in the iterative algorithm.

FIGURE 8

FIGURE 8. Validation of independence (A) grid size. (B) time step.

5.1 Influence of Concrete on Hydraulic Pressure

As known from Eq. 22 and Eq. 31, FSI responses to water hammer greatly depend on parameters of concrete, including its Poisson’s ratio, thickness-radius ratio, elastic modulus and density. Hydraulic pressure upstream the valve fluctuates with these parameters are shown in Figures 9A–D, respectively.

FIGURE 9

FIGURE 9. Dynamic pressure varies with the parameters of concrete (A) Poisson’s ratio. (B) Thickness-radius ratio. (C) Elastic modulus. (D) Density.

It is clear that the Poisson’s ratio and concrete density hardly affect the hydraulic pressure as the two curves are almost overlapped, as seen in Figures 9A,D. Firstly, coefficients in Eqs 10, 11 are the function of ν_s from their definitions in the Appendix. But the difference of R⋅p|_r=R - (R + e)⋅P_b can be ignorable. Moreover, ν_s mainly contributes to the axial shear force between concrete and pipe, while little affects the pressure. With respect to the concrete density, ρ_s only dominates the inertia of concrete, then affects the frictional force on the outer interface. However, hydraulic pressure hardly changes when ρ_s becomes ten times larger. This reveals that the axial vibration of concrete plays a negligible role in FSI responses.

In contrast, the thickness-radius ratio and elastic modulus play an important role in the FSI response. By extension, larger E_s alleviates the pumping effect on fluid and further weakens the FSI response. Hydraulic pressure fluctuates more moderately in this case and decays faster, as seen in Figure 9C. Besides, E_s and e_s/R_s both influences the pressure wave speed, then change the fluid inner pressure. Figure 10 shows the wave speed changing with the product of E_s and e_s/R_s.

FIGURE 10

FIGURE 10. Pressure wave speed varies with E_se_s/R_s.

In the graph when E_se_s/R_s < 10⁶ Pa, pressure wave speed is almost constant and the pumping effect can be ignored. So the pipe encased in concrete is approximate to an exposed one, and the four-equation model is suitable in this case. When E_se_s/R_s > 10¹¹ Pa, influences of wave speed and pumping effect are also negligible. In this case, the FSI response is weak, so the two-equation model is suitable. When E_se_s/R_s increases from 10⁶ Pa to 10¹¹ Pa, the wave speed changes dramatically as well as the FSI responses. Notably, the E_se_s/R_s of concrete is within this range. Thus, studying the FSI in a piping system embedded in concrete is of great importance.

5.2 Pipe-Wall Stress

In this study, pressure gradients applied to the interface cause pipe wall expansion or contraction, and further compress or relax the ambient concrete. Thus, P_b is produced and subsequently weakens the FSI responses. The dynamic pressure on the inner and the outer interfaces are separately shown in Figures 11A,B.

FIGURE 11

FIGURE 11. Pressure fluctuates with E_se_s/R_s (A) On the inner interface. (B) On the outer interface.

Stronger pulsations are caused by larger E_se_s/R_s on both inner and outer interfaces. With larger P_b, a weaker pumping effect is obtained. Hoop stress and radial stress of pipe wall become smaller due to restrictions by concrete from Eqs 1, 2. The result is the piping system gets safer. Moreover, the difference between R⋅p|_r=R and (R + e)⋅P_b is smaller in the case of larger E_se_s/R_s, which can further lead to weaker axial stress pulsation in the pipe wall, as shown in Figure 12.

FIGURE 12

FIGURE 12. Axial stress of pipe wall fluctuates with E_se_s/R_s.

The six-equation model derived in this paper aims to investigate FSI subject to water hammer. Due to the deformation of structures, FSI response to some extent weakens the stress pulsation in the pipe wall. As shown in Figure 13, both the hoop and radial stress are larger using the classical model without FSI.

FIGURE 13

FIGURE 13. Stress increases with p|_r=R (A) Hoop stress. (B) Radial stress.

5.3 Pipe-Wall Motion

In this segment, the axial vibration of the system is analyzed. Axial displacement of the pipe wall relates to the fluid and concrete, while the concrete plays the damping role. Signals of axial displacement in the time domain and frequency domain are shown in Figures 14A,B, separately. When FSI is taken into account, the pipe-wall motion is found to relate to the pressure wave in fluid and stress wave in solid. These two waves propagate throughout the pipeline system in the axial direction, contributing to the vibration amplitude in Figure 14A. In radial direction, pipe expansion or contraction is confined by the ambient concrete, so the radial vibration of the pipe wall is weak.

FIGURE 14

FIGURE 14. Dynamic axial displacement (A) Time domain. (B) Frequency domain.

Pipe vibration is mainly caused by hydraulic pressure fluctuation, especially when FSI is included. All vibration modes in this case are system modes. However, pipe deformation mode can be identified by data processing. The first eight modes of pipe with respect to different models are shown in Table 1. Data in the right column correspond to the local peaks in Figure 14B. According to calculation results, the frequency of pressure wave is 35.2 Hz, the frequency of stress wave in pipe wall is 129.3 Hz, and the frequency of stress wave in concrete is 99.6 Hz.

TABLE 1

TABLE 1. First eight modes of pipe with respect to three models.

The two-equation model predicts the highest frequency value at each mode without FSI. Compared with the four-equation model, constraints by concrete enlarges the frequency value in the six-equation model. Wave speed in this model is faster, in accord with the curve tendency in Figure 10. Additionally, Pipe and concrete absorb more energy from water hammers so that FSI responses are alleviated. Similar results are obtained in Figure 7. Comparison of the three models shows that inherent frequency of the piping system is changed when different factors are taken into account, like the concrete. On the other, the pressure surge acts as an unstable source during FSI responses. When parameters of the concrete change, nature frequency of the pressure pulsation is changed accordingly. This situation may lead to undesirable resonance, especially in a more complex pipeline layout.

6 Conclusion

This work studies the FSI responses of a piping system embedded in concrete during water hammers. Based on area average treatment on a cross-section, a six-equation model is derived to describe the fluid-pipe-concrete interaction. This model is subsequently solved by numerical iteration and results are compared with experimental data and classical models, showing good agreements.

In a reservoir-pipe-valve system, hydraulic pressure with respect to concrete parameters is discussed. Numerical results reveal that Poisson’s ratio and density of concrete hardly contribute to FSI responses, while elastic modulus and thickness-radius ratio play crucial roles. More specifically, E_s mainly weakens the pipe pumping effect and e_s/R_s determines pressure wave speed. Larger E_s and smaller e_s/R_s result in more mild pressure fluctuations. When E_se_s/R_s < 10⁶ Pa, this six-equation model is approximate to the four-equation model, namely an exposed pipeline system. When E_se_s/R_s > 10¹¹ Pa, constraints from concrete are large enough so that the classical two-equation model is suitable.

Structural vibration caused by FSI mitigates stress pulsations separately in circumferential, radial and axial directions of the pipe wall. These pulsations are more moderate in the case of weaker FSI responses. Without FSI, the classical water-hammer theory may lead to an unintended large safety margin and uneconomic design. It should be noted that weaker FSI responses cause higher structural frequency in a water conveyance system, which justifies more explorations in the design and operation stage of the piping system.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

YC: conceptualization, computation, visualization, writing−original draft; CZ: data analysis, writing−original draft and editing; QG: methodology, formal analysis, validation, writing−review and editing; JZ: methodology, resources, supervision; YF: supervision, writing−review and editing; KX: funding acquisition, assistance with analysis.

Funding

This work was supported by the National Natural Science Foundation of China (Grant No. 12002150), the Natural Science Foundation of Jiangsu Province (BK20201041), the High−education Natural Science Foundation of Jiangsu Province (21KJB570001), and the Research Foundation for High−level Talents of Nanjing Institute of Technology (YKJ202131).

Conflict of Interest

YC is employed by the Company Haitian Plastics Machinery Group Co., Ltd.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Alaei, D., Kwon, Y. W., and Ramezani, A. (2019). Fluid-structure Interaction on Concentric Composite Cylinders Containing Fluids in the Annulus. Multiscale Multidiscip. Model. Exp. Des. 2 (3), 185–197. doi:10.1007/s41939-019-00044-3

CrossRef Full Text | Google Scholar

Cardiff, P., Tuković, Ž., Jasak, H., and Ivanković, A. (2016). A Block-Coupled Finite Volume Methodology for Linear Elasticity and Unstructured Meshes. Comput. Struct. 175, 100–122. doi:10.1016/j.compstruc.2016.07.004

CrossRef Full Text | Google Scholar

Daoxiang, Z., and Weinlin, Z. (2006). Foundation of Engineering Elasticity. Hefei, China: Hefei University of Technology Press.

Google Scholar

De Santis, D., and Shams, A. (2019). Scaling of Added Mass and Added Damping of Cylindrical Rods by Means of Fsi Simulations. J. Fluids Struct. 88, 241–256. doi:10.1016/j.jfluidstructs.2019.05.011

CrossRef Full Text | Google Scholar

Ferras, D., Manso, P. A., Schleiss, A. J., and Covas, D. I. C. (2017). Fluid-Structure Interaction in Straight Pipelines with Different Anchoring Conditions. J. Sound Vib. 394, 348–365. doi:10.1016/j.jsv.2017.01.047

CrossRef Full Text | Google Scholar

Hu, S., Sun, Y., Xue, X., and Huang, Y. (2019). Calculation Model for Bar-Wrapping during Prestressing of an Embedded Bar-Wrapped Cylinder Concrete Pressure Pipe. Thin-Walled Struct. 139, 39–45. doi:10.1016/j.tws.2019.02.036

CrossRef Full Text | Google Scholar

Huang, W.-X., and Alben, S. (2016). Fluid-structure Interactions with Applications to Biology. Acta Mech. Sin. 32 (6), 977–979. doi:10.1007/s10409-016-0608-9

CrossRef Full Text | Google Scholar

Karakouzian, M., Karami, M., Nazari-Sharabian, M., and Ahmad, S. (2019). Flow-Induced Stresses and Displacements in Jointed Concrete Pipes Installed by Pipe Jacking Method. Fluids 4 (1), 34. doi:10.3390/fluids4010034

CrossRef Full Text | Google Scholar

Lai, M. H., Liang, Y. W., Wang, Q., Ren, F. M., Chen, M. T., and Ho, J. C. M. (2020). A Stress-Path Dependent Stress-Strain Model for Frp-Confined Concrete. Eng. Struct. 203, 109824. doi:10.1016/j.engstruct.2019.109824

CrossRef Full Text | Google Scholar

Lee, C. K., Chiew, S. P., and Jiang, J. (2012). Residual Stress Study of Welded High Strength Steel Thin-Walled Plate-To-Plate Joints, Part 1: Experimental Study. Thin-Walled Struct. 56, 103–112. doi:10.1016/j.tws.2012.03.015

CrossRef Full Text | Google Scholar

Li, S.-J., Liu, G.-M., and Chen, H. (2012). Pressure Wave Propagation Characteristics in Fluid-Filled Pipes with Fluid-Structure Interaction. J. Vib. Shock 31 (24), 177–182. doi:10.13465/j.cnki.jvs.2012.24.029

CrossRef Full Text | Google Scholar

Li, S.-j., Liu, G.-m., and Kong, W.-t. (2014). Vibration Analysis of Pipes Conveying Fluid by Transfer Matrix Method. Nucl. Eng. Des. 266, 78–88. doi:10.1016/j.nucengdes.2013.10.028

CrossRef Full Text | Google Scholar

Liu, G., and Li, Y. (2011). Vibration Analysis of Liquid-Filled Pipelines with Elastic Constraints. J. Sound Vib. 330 (13), 3166–3181. doi:10.1016/j.jsv.2011.01.022

CrossRef Full Text | Google Scholar

Mahmoodi, R., Zolfaghari, A., and Minuchehr, A. (2019). Laplace Transform Finite Volume Modeling of Water Hammer along Fluid-Structure Interaction. Comput. Math. Appl. 77 (10), 2821–2832. doi:10.1016/j.camwa.2019.01.014

CrossRef Full Text | Google Scholar

Mirjavadi, S. S., Forsat, M., and Badnava, S. (2020). Nonlinear Modeling and Dynamic Analysis of Bioengineering Hyper-Elastic Tubes Based on Different Material Models. Biomech. Model. Mechanobiol. 19 (3), 971–983. doi:10.1007/s10237-019-01265-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Rajbamshi, S., Guo, Q., and Zhan, M. (2020). Model Updating of Fluid-Structure Interaction Effects on Piping System. Dyn. Substruct. 4, 133–139. doi:10.1007/978-3-030-12184-6_12

CrossRef Full Text | Google Scholar

Ramírez, L., Nogueira, X., Ouro, P., Navarrina, F., Khelladi, S., and Colominas, I. (2018). A Higher-Order Chimera Method for Finite Volume Schemes. Arch. Comput. Methods Eng. 25 (3), 691–706. doi:10.1007/s11831-017-9213-8

CrossRef Full Text | Google Scholar

Riedelmeier, S., Becker, S., and Schlücker, E. (2014). Measurements of Junction Coupling during Water Hammer in Piping Systems. J. Fluids Struct. 48, 156–168. doi:10.1016/j.jfluidstructs.2014.03.001

CrossRef Full Text | Google Scholar

Sinha, J. K., Singh, S., and Rama Rao, A. (2001). Finite Element Simulation of Dynamic Behaviour of Open-Ended Cantilever Pipe Conveying Fluid. J. Sound Vib. 240 (1), 189–194. doi:10.1006/jsvi.2000.3113

CrossRef Full Text | Google Scholar

Sun, Y., Hu, S., Huang, Y., and Xue, X. (2020). Analytical Stress Model for Embedded Bar-Wrapped Cylinder Concrete Pressure Pipe under Internal Load. Thin-Walled Struct. 149, 106540. doi:10.1016/j.tws.2019.106540

CrossRef Full Text | Google Scholar

Tijsseling, A. S. (2019). An Overview of Fluid-Structure Interaction Experiments in Single-Elbow Pipe Systems. J. Zhejiang Univ. Sci. A 20 (4), 233–242. doi:10.1631/jzus.A1800564

CrossRef Full Text | Google Scholar

Tijsseling, A. S. (2007). Water Hammer with Fluid-Structure Interaction in Thick-Walled Pipes. Comput. Struct. 85 (11-14), 844–851. doi:10.1016/j.compstruc.2007.01.008

CrossRef Full Text | Google Scholar

Ting, Z., Yongzhen, L., Zhiqiang, Y., Chiaming, F., and Tsunghan, L. (2017). Axial Vibration Response of Viscoelastic Pipe Conveying Fluid Induced by Water Hammer. J. Vib. Eng. 30 (2), 241–248. doi:10.16385/j.cnki.issn.1004-4523.2017.02.009

CrossRef Full Text | Google Scholar

Vardy, A. E., Fan, D., and Tijsseling, A. S. (1996). Fluid-Structure Interaction in a T-Piece Pipe. J. Fluids Struct. 10 (7), 763–786. doi:10.1006/jfls.1996.0052

CrossRef Full Text | Google Scholar

Wu, J.-S., and Shih, P.-Y. (2001). The Dynamic Analysis of a Multispan Fluid-Conveying Pipe Subjected to External Load. J. Sound Vib. 239 (2), 201–215. doi:10.1006/jsvi.2000.3119

CrossRef Full Text | Google Scholar

Wu, Z. (20172017). Elasticity. Beijing: Beijing Institute of Technology.

Google Scholar

Yang, K., Li, Q., and Zhang, L. (2004). Longitudinal Vibration Analysis of Multi-Span Liquid-Filled Pipelines with Rigid Constraints. J. Sound. Vib. 273 (1-2), 125–147. doi:10.1016/S0022-460X(03)00422-X

CrossRef Full Text | Google Scholar

Zanganeh, R., Ahmadi, A., and Keramat, A. (2015). Fluid-structure Interaction with Viscoelastic Supports during Waterhammer in a Pipeline. J. Fluids Struct. 54, 215–234. doi:10.1016/j.jfluidstructs.2014.10.016

CrossRef Full Text | Google Scholar

Appendix

$a=\frac{-\frac{E_s}{E}\nu \frac{2R^2}{{\left(R+e\right)}^2-R^2}\frac{2{\rho }_t{\left(R+e/2\right)}^{2}e{\rho }_t}{{\rho }_f{R}^2+2{\rho }_t\left(R+e/2\right)e}}{1+\frac{E_s}{E}+\frac{E_s}{E}\nu \frac{2R^2}{{\left(R+e\right)}^2-R^2}-{\nu }_s\frac{R_s^2+{\left(R+e\right)}^2}{R_s^2-{\left(R+e\right)}^2}-\frac{2{\rho }_t{\left(R+e/2\right)}^{2}e{\rho }_t}{{\rho }_f{R}^2+2{\rho }_t\left(R+e/2\right)e}}$

$b=\frac{-\frac{E_s}{E}\nu }{1+\frac{E_s}{E}+\frac{E_s}{E}\nu \frac{2R^2}{{\left(R+e\right)}^2-R^2}-{\nu }_s\frac{R_s^2+{\left(R+e\right)}^2}{R_s^2-{\left(R+e\right)}^2}-\frac{2{\rho }_t{\left(R+e/2\right)}^{2}e{\rho }_t}{{\rho }_f{R}^2+2{\rho }_t\left(R+e/2\right)e}}$

$c=\frac{{\nu }_s}{1+\frac{E_s}{E}+\frac{E_s}{E}\nu \frac{2R^2}{{\left(R+e\right)}^2-R^2}-{\nu }_s\frac{R_s^2+{\left(R+e\right)}^2}{R_s^2-{\left(R+e\right)}^2}-\frac{2{\rho }_t{\left(R+e/2\right)}^{2}e{\rho }_t}{{\rho }_f{R}^2+2{\rho }_t\left(R+e/2\right)e}}$

$d=\frac{2{\rho }_t{\left(R+e/2\right)}^{2}e+2{\left(R+e\right)}^{2}a}{{\rho }_f{R}^2+2{\rho }_t{\left(R+e/2\right)}^{2}e}$

$f=\frac{2{\left(R+e\right)}^{2}b}{{\rho }_f{R}^2+2{\rho }_t{\left(R+e/2\right)}^{2}e}$

$h=\frac{2{\left(R+e\right)}^{2}c}{{\rho }_f{R}^2+2{\rho }_t{\left(R+e/2\right)}^{2}e}$

$\mathit{A}=\left(\begin{array}{cccccc}{\rho }_f& & & & & \\ & & & 1/K+2\alpha & 2\beta & 2\gamma \\ & {\rho }_t& & & & \\ & & & \nu k& \nu l-1& \nu m\\ & & {\rho }_s& & & \\ & & & \frac{{\nu }_s{R}_s^{2}a}{\left(R_s+e_s/2\right)e_s}& \frac{{\nu }_s{R}_s^{2}b}{\left(R_s+e_s/2\right)e_s}& \frac{{\nu }_s{R}_s^{2}c}{\left(R_s+e_s/2\right)e_s}-1\end{array}\right)$

$\mathit{B}=\left(\begin{array}{cccccc}{\rho }_f\left|V_i^n\right|& & & 1& & \\ 1& & & & & \\ & & \frac{E_s\left(R+e\right)\mathrm{\Delta }t}{2\left(1+{\nu }_s\right)\left(R+e/2\right)e}& & 1& \\ & 1& & & & \\ & & -\frac{E_s\pi \left(R+e\right)\mathrm{\Delta }t}{2\left(1+{\nu }_s\right)A_s}& & & 1\\ & & 1& & & \end{array}\right)$

$\mathit{S}=\left(\begin{array}{c}0\\ 0\\ \frac{E_s\left(R+e\right)\mathrm{\Delta }t}{2\left(1+\nu \right)\left(R+e/2\right)e}\frac{\partial {\displaystyle \sum _{n=0}^{\left(n-1\right)\mathrm{\Delta }t}{\left(u_s\right)}_n}}{\partial z}\\ 0\\ -\frac{E_s\pi \left(R+e\right)\mathrm{\Delta }t}{2\left(1+\nu \right)A_s}\frac{\partial {\displaystyle \sum _{n=0}^{\left(n-1\right)\mathrm{\Delta }t}{\left(u_s\right)}_n}}{\partial z}\\ 0\end{array}\right)$