FSI ANALYSIS OF SUBMARINE OUTFALL

In the scope of this study, main pipe of the diffuser, risers, ports, internal and external environments forming the discharge system which is used in application are modelled by Finite Elements Analysis (FEA) program to obtain discharge and structural behaviour. The last two spans of the system (20 m) and four ports on these spans are investigated. While the diameter and geometry of the risers and ports remain constant, the diffuser pipe is modelled in three different ways. These are constant sectioned (Model 1), contracting with sharp edge entrance sectioned (Model 2) and gradually contracting sectioned (Model 3) respectively. Among them, only Model 1 is treated as Single Degree of Freedom (SDOF) system and it is simulated by FEA to verify FEA solver in the first place. After structural suitability is confirmed, rest of the models are analysed to determine reaction forces and stresses. The discharge is performed as unsteady external flow as well as steady external flow assumption which is widely used in external flow model in the literature. The discharge analyses are performed in two different ways to verify FEA program. Iterative method is accompanying to FEA program. As a result of this study, proper model for structural and discharge behaviour and external flow effects on discharge velocities are obtained.


68
E (ML −1 T −2 ) Young modulus F (MLT −2 ) Force f (1) Friction losses g (LT -2 ) Gravity H (L) Wave height ht (L) Total head I (ML 2 ) Moment of inertia I J (MLT −2 )  Internal force i (1) Port number k (MT −2 ) Stiffness L (L) Diffuser length LW (L) Wave length m (M) Mass N (1) Total number of risers p (ML −1 T −2 ) Pressure Flow rate of total discharge q (L 3 T -1 ) Flow rate of port T (T) Wave period t (T) Time u (LT -1 )  Fluid velocity component at x direction V (LT −1 )  Average fluid velocity v (LT -1 )  Fluid velocity component at y direction w (LT -1 )  Fluid velocity component at z direction X (L) Displacement of diffuser X (LT −1 ) Velocity of diffuser Acceleration of diffuser y (L) Geometrical head α (1) Port parameter 0 Grüneisen ratio Slope of the Us−Up curve Square of natural frequency μ (ML −1 T −1 ) Dynamic viscosity Poisson ratio Natural frequency Mass density Local losses

Introduction
During the last few decades, production of hot salty water (brine) in the world has increased abundantly due to rapid increase in various industrialized and mining processes.Submarine outfall diffusers occupy 41 % of total brine disposal capacity due to high dilution capabilities [1].Besides, submarine outfall systems and diffusers are mostly used in many industrial applications from jet engineers to manifold and air conditioning systems [2][3][4][5][6].
Gradually contacting and constant sections are widely implemented in main pipe design of diffusers.Different geometries can be used in ports together with diffuser pipes.In addition to these section types, contracting with sharp edge entrance sections are generally used in buried diffusers [7][8][9].Different geometries can be used in ports as well as diffuser pipes.
Internal and external hydraulic characteristics must be taken into account while determining the diffuser geometry.However, it is known that internal flow in liquid-filled pipe systems is not steady.It is generally assumed as steady by treating internal flow in submarine outfalls.Discharge velocity varies with ambient pressure and offshore environment flow parameters (velocity, pressure etc.) which are time-varying unlike usage in wide range of diffuser applications.This assumption makes the calculation simpler.According to literature survey, it is stated that turbulence effects are negligible except water hummer situation in discharge systems.It is known that effective and controlled usage of submarine outfalls can prevent systems from water hummer effects.Both uniform discharge distribution along the diffuser and preventing salty water intrusion are the most important conditions for discharge systems.Turbulence effects are said to be not effective on these conditions by references [7][8][9][11][12][13].
Submarine outfall models are composed of fluid and structure domains that interacts each other.Fluid-structure interaction (FSI) techniques are commonly implemented to model dynamic behaviour of both fluid and structure [14][15][16][17].FSI analyses are including Computational Fluid Dynamics (CFD) to realise the fluid flow (CFD).In this day and age, however there are many experimental techniques, an alternative method is provided by (CFD) because experimental setup is extremely expensive and laborious and down-scale models are not accurate enough at times [18].
In this study, firstly three different diffuser models which are widely used in applications have been investigated in terms of structural behaviour to determine the sufficient one.After performing structural analysis, effects of changes in diffuser pipe geometry and ambient surrounding diffusers on discharge velocity are investigated.For this purpose, as it is seen Figure 1, three different diffuser pipe sections are modelled; such as constant (Model 1), contracting with sharp edge entrance (Model 2) and gradually contracting (Model 3) respectively.[19,20], CFD results are confirmed by computing energy and continuity equations given by Eqns.(15,16) [7].Flow domains consist of internal and external flows.Internal flow is modelled to be steady.Linear Wave Theory is used to simulate unsteady external flow.Wave is occurred in intermediate water depth region according to wave parameters.As it is stated in [21], underneath the free surface, wave-induced pressure oscillations reduce with depth below the free surface.Since the diffuser is in the sea bottom, effect of free surface is ignored in the transition regions as mentioned in most of the studies in the literature.External flow is separately modelled as steady besides unsteady flow.In this manner, the steady flow assumption, as performed in previous studies, changes the flow characteristics according to unsteady flow that will be revealed.

Methodology
In the scope of this study, it is aimed to perform computer aided numerical FSI analysis of submarine outfall diffusers under internal and external flows.The fluid part of the analysis results are compared with iterative technique that is based on the equilibrium of pressure at the same point obtaining from different points.The comparison is made for all geometric models for both conditions; steady and unsteady outflow.The structure part is compared with SDOF.The solid models under consideration consist of three different diffuser geometries.These are compared by considering reaction force and stress values.All geometries are transporting internal steady flow interacting with external flow.In contrast to previous studies, external flow is modelled not only steady but also unsteady.Mentioned diffusers are discharging salty cooling water with the flow rate of 0.173 m 3 /s, having the same density with ambient to sea.Linear Wave Theory is charged to model unsteady external flow.Ambient flow parameter effects on effluent flow parameters are determined by modelling external flow in two different types.The material for the models is steel, with Young modulus (E) of 2.1x10 11 N/m 2 , Poisson ratio () of 0.30 and mass density (s) of 7850 kg/m 3 .The models are divided into small elements in finite elements method to perform and analyse the complex models.10-node modified tetrahedron elements (C3D10M), which are compatible with contact problems, are utilized in analyses.Distances between meshes are taken as 0.01 m on ports and riser which is the same value of wall thickness, and 0.05 m on diffuser pipe.Smaller values than these values cause duration and volume problems.Ultimate mesh structures and port numbers of models are presented in Figure 2.
u, v and w are the velocity, gx, gy and gz are gravitational components at the x, y and z directions respectively.ρf is the mass density, μ is the dynamic viscosity and p is the pressure.Turbulence terms are ignored as well as former studies [23][24][25][26][27].
Internal and external flows are created in the same geometry as shown in Fig. 3 Fluid domains Fig. 4 Ports mesh structures and node numbers FC3D4 (4-node modified tetrahedron) typed members which are proper for FSI problems are used in the analyses.Distances between meshes are taken as 0.01 m on ports, which is the same value of wall thickness, and 1 m on the rest of the geometry as presented in Figure 4 for Model 2. In Model 1, 23125 nodes and 105102 elements, in Model 2, 22895 nodes and 103897 elements and in Model 3, 24053 nodes and 109152 elements constitute the domains.The analysis is performed by an applied boundary condition on diffuser pipe as the fluid inlet velocity of 0.956 m/s corresponding to flow rate.Simultaneously, equation of velocity that represents the Linear Wave Theory is applied to outlet domain, is given below.
In this paper, employed parameters are taken into account respectively: water depth (d) is 20 m, wave period (T) is 8 s and wave height (H) is 2.50 m.Wave length (LW=95.72 m) is calculated by considering these parameters.Following this, steady external velocity (u=0.5 m/s) which is introduced to outlet domain to observe external flow effects on discharge velocities.

Modelling of FSI
The first step of FSI problem is tasked to determine the contact surfaces as seen in Figure 2. By determining the contact surfaces, where the forces are transferred from fluid to structure and deformations are transferred from structure to fluid is identified.Structural and fluid equations are solved independently.Finite Elements program employs Eq. (1-3) for fluid solver to obtain pressure forces.Subsequently, the (Eq.5) is used to obtain displacement values by Explicit analysis in which the values are transferred to fluid by FSI technique.

  
In (Eq.5), m NJ is the mass matrix, X is acceleration, t is time, F J symbolizes external applied load vector transferred from CFD, I J is internal force vector which is occurred by stresses in the elements.The equations of motion for the body are integrated due to equations given below.
N X and N X are degree of freedom, (N) of displacement and velocity components in the (Eqns.6-7) respectively.The nodal accelerations are calculated by using (Eq.8).
Velocity and displacement values can be obtained after determining accelerations.Modal analyses are also performed simultaneously to find natural frequencies beside explicit analysis.The finite element of the model is given by the matrices in (Eq.9).

 
Where  is square of natural frequency [28].Lanczos Method is utilized solving matrices [21].In this paper, the analyses are completed by 2e -5 time increment for 8 s and the In (Eq.10), I is moment of inertia and (z) is shape function given below.The right side of (Eq.10) which is given in parenthesis is Morrison Equation that is used to represent externally applied wave forces.Force components include the force coefficients as CD is 2.40 and CM is 0.70.(z), given by (Eq.11), should satisfy the geometric boundary conditions,(0)='(0)=0, (L)='(L)=0.The final equation of motion of 1-DOF model is obtained by determining shape function as given in (Eq.12).
157230 476356 In the equation above, F(t) is resultant applied force that is extracted by computing right side of (Eq.10).Time varying numerical value of F(t) is presented by Figure 5.In this paper, Runge-Kutta method is performed to solve (Eq.10) via [29] under dynamic initial boundary conditions for t=0 X(0)=0, (0) =0

X
. This method evaluates the simple relationships given below at the beginning, middle and end of each overall time step (t) [30].
  After time varying displacement function, X(t), is derived.So, time and location varying displacements can be obtained by multiplying X(t) with (z).Natural frequency value (0) is obtained from (Eq. 14) where k and m can be described as stiffness and mass.
Engin GÜCÜYEN, R. Tuğrul ERDEM, FSI Analysis of Submarine Outfall Ümit GÖKKUŞ 74 2.3.2Verification CFD model using iterative method Hydraulic analysis, which is based on the energy and continuity equations, is carried on simultaneously by FEM program.In this analysis, well known Energy Equation given in (Eq.15) is utilized.p, V, g, , y and ht represent pressure, velocity, gravity, specific weight, geometrical head and total loses respectively.
22 11 1 This equation is implemented from end to beginning of the diffuser starting with the nodes i and i-1 in the same streamline.Designation of the node numbers which are utilized in analysis is seen in Figure 4.In this method, energy equation is applied between the two successive nodes (i and i-1) along a streamline following the diffuser pipe centreline and between nodes where one is on the diffuser pipe centreline (i) other one is in the same streamline on the port.In this way, it is aimed to obtain the same value for pressure (pdi).
Energy equation is applied along a streamline following the diffuser pipe centreline results in (Eq.16).22 22 11 11 11 In (Eq.16),A is cross section area and q is flow rate of port.Diffuser pressure (pd,i) equals the sum of upstream the port/riser branch with the known downstream diffuser pressure (pd,i-1), the known static pressure difference due to the elevation difference, the dynamic pressure difference and the known losses occurring in the main diffuser pipe.The losses are divided into friction losses (f) that are calculated Darcy Weisbach Equations via Moddy Diagram and local losses () like bends and diameter changes or the passage of a branch opening.In this study, the elevation differences equal to zero due to zero bottom slope.Total losses between mentioned nodes are given as fallows.
The pressure value (pdi) can be found by writing the Energy Equation between node (i) and one in the same streamline on the port contacting the ambient.It is given by (Eq.18).c ,i p ,i   a,i  p,i  d ,i  d ,i

 
Pressure in diffuser pipe (pd,i) equals to sum of the upstream diffuser pressure with the ambient pressure (pa,i), the static pressure difference due to the elevation difference between diffuser centrelines and port centreline, dynamic pressure difference between the diffuser and one single port and the losses occurring in all pipe segments between these points.Cc is the jet contraction coefficient given by (Eq.19) for bell mouthed ports.

 
The total loses (ht-dp) between mentioned nodes are given below.

FSI Analysis of Submarine Outfall
Engin GÜCÜYEN, R. Tuğrul ERDEM Ümit GÖKKUŞ 75 2 11 1 2 f i p ,i , j p ,i , j r ,i , j r ,i , j p ,i , j r ,i , j p ,i , j d ,i , j r ,i , j r ,i , j Diffuser pipe, riser and port are parameterized by d, r and p respectively in equations above.The same pressure value must be obtained from (Eq.16) and (Eq.18) separately.Therefore, the two equations are equalled to each other.New equation obtained from this operation is given by (Eq.21).
  αi is a parameter with αi = 1/(number of ports at a riser at position i).In this study, ambient pressure is modelled for both dynamic and static conditions contrary to previous studies.Equation of pressure for Linear Wave Theory is given as fallows.
In this study, (Eq.21) is computed [29] for total head.First estimation is used as a starting value and further iterations lead to the final value.At the first port/riser on the seaward side (i=1) an initial discharge q 1 is estimated, for example q 1 = Q/N with Q=total discharge and N=total number of risers.(Eq.16) then allows to calculate the first internal pressure of the diffuser p d,1 . The further discharges q 2 until q N are calculated by using (Eq.21).
A final application of (Eq.16) allows to calculate p d,N+1 , the necessary pressure at the headwork to drive the system.While the internal steady flow is conveyed to unsteady ambient, the internal flow pattern is disrupted and becomes unsteady.This situation induces velocity dissimilarities.As a result of this study, the quantities of these dissimilarities are determined and time varying discharge velocities are obtained in the end.

Results
The study includes evaluations of ABAQUS-SDOF to determine structural behaviour under unsteady external flow and ABAQUS-MATLAB iterations to obtain discharge velocities of submarine outfall diffusers under unsteady and steady external flows.Both fluid and structure results can be obtained from ABAQUS/FSI analyses.Equation of motion of SDOF system is evaluated by Runge-Kutta method.Thus, maximum displacement value on the main pipe is given Table 1.Similar values obtained from finite elements analysis are also given in the same table.In addition to displacement values, first natural frequency values are comparatively presented in Table 1.After one of the FEM structural results are verified, all FEM models are comparatively investigated in terms of Von-Mises Stresses and reaction forces.Visual presentations of the results are given in following figures.As it is presented in Figure 6, maximum value of Von Mises Stress reaches the limit of 1.9810 7 N/m 2 on Model 1, 2.7710 7 N/m 2 on Model 2 and 3.0410 7 N/m 2 on Model 3.According to the reaction forces, it can be said that maximum value on Model 1 is 1.0310 3 N, on Model 2 is 1.7010 3 N and on Model 3 is 1.5610 3 N. Node varying discharge velocities on the pipe axis at certain nodes (i-1, i, i+1, i+2) are obtained by ABAQUS via FSI as shown in Figure 8 for Model 2. These discharge velocities are converted to average ones to compare the analysis results each other.Same outputs are computed in iterative manner by using [29].In this case, average discharge velocities on the ports are obtained.In Figure 8, lines in sinusoidal form, given by Eq. legends, are derived from iterating the (Eq.21).However, the other ones are derived from FEM program.Node numbers on ports of Model 2 are seen in Figure 4. Velocity values are given in Figure 8 in accordance with these nodes.Having the similar results for velocities shows that uniform discharge between ports are satisfied.This situation shows that the diffusers work sufficiently.Average discharge velocities are presented in Table 2 and Table 3 from FEM program and iterating the (Eq.21) for different geometries where the external flow is unsteady and steady.

Conclusions
In this study structural and discharge behaviour of submarine diffusers are investigated simultaneously by FEM Program via FSI.Three different structural models (Model 1, 2, 3), used in applications are utilized when examining the models in terms of geometry.Von-Mises stresses and reaction forces are studied in terms of structural behaviour to indicate the proper model.At the same time, these models are analysed for discharge velocities under steady and unsteady external flow conditions.Unsteady flow is characterized by Linear Wave Theory.Analyses are performed by FEM Program (ABAQUS).Verification of structural results of FEM is confirmed by SDOF model.Subsequently FEM/CFD results are confirmed by iterating the (Eq.21).
FEM/Explicit solutions of Model 1 differ from the approximate solutions obtained from Runge-Kutta method.In this study, it is observed that displacement and first natural frequency values differences are not exceeding 4.16% and 3.91% respectively.The values would be different if another shape function satisfying the boundary conditions was chosen instead of the one given by (Eq.11).After FEM solver is examined through Model 1, all models are analysed.Minimum values of Von-Mises stresses and reaction forces are obtained from Model 1.On the other hand, Von-Mises stresses reach maximum values on Model 3.While stresses are concentrated at the pipe-riser connections and supports for Model 1 and Model 2, stress propagation is observed for Model 3. Unlike stress results, maximum reaction force occurs on Model 2 where pipe connection is sudden.Reaction force reaches maximum value at sudden contraction zone in Model 2. According to both the stress and reaction force, Model 1 is the safest one among all.Comparing to other models, the most impractical model for FEM solvers is Model 3 due to consisting of more nodes and elements and the most sufficient for fabrication is Model 1 due to having the simplest geometry.
The second purpose of this study is to determine the effects of variation of external flow and structural model on discharge velocities.External flow is modelled under two different conditions as steady and unsteady.Uniform discharge distribution along the diffuser and salty water intrusion conditions are said to be provided by reference to the results.Different ports which effect the discharge velocities were examined in previous studies.In this study, it is aimed to expose whether discharge velocities are effected by diffuser pipe geometries or not.According to structural models, the results are said to be unvaried apparently.Geometric structure of the diffuser pipe has no effects on the discharge velocities.
Compability of CFD and iterative analysis is observed in Tables 2 and 3.In case of only discharge velocities are needed, iterative technique would be sufficient due to taking less time than FEM program.FEM program shall be performed when visual results are need as it is seen Figure 9.
As a conclusion, it is suggested to model diffuser main pipe as constant sectioned according to presented results of this study.In terms of fluid results, it is observed that unsteady external flow has no significant effects on discharge velocities.Although, remarkable effects are not detected, in the sense of water hammer effects the unsteady flow shall be taken into account.Whether unsteady flow causes water hammer effect or not shall be examined in forthcoming studies.Finally, in cases that do not require detailed analysis of external flow shall be studied as steady flow.

Fig. 1
Fig. 1 Diffuser geometries All sections are modelled interacting with flow domains.FSI calculations of diffuser and internal-external flows are performed by using ABAQUS FEA program.Structural modelling of diffuser is performed by ABAQUS/Explicit and internal-external flow modelling is created ABAQUS/CFD simultaneously.While the structural verification of FEM is determined by SDOF model[19,20], CFD results are confirmed by computing energy and continuity equations given by Eqns.(15,16) [7].Flow domains consist of internal and external flows.Internal flow is modelled to be steady.Linear Wave Theory is used to simulate unsteady external flow.Wave is occurred in intermediate water depth region according to wave parameters.As it is stated in[21], underneath the free surface, wave-induced pressure

2. 1
Modelling of diffusers Two-span diffuser with 20 m length (L) is modelled with ABAQUS/CAE [22].Fixed supports with a length of 0.40 m are used.Vertical ports having 1 m length are placed as the distance between them to be 5 m.The diameter and geometry of the risers and ports of the diffuser remain constant.However, the diffuser pipe is modelled in three different ways.While diameter of the two-span main diffuser pipe is 0.50 m in Model 1, diameter of the main pipe (D) is 0.50 m in the first span and 0.40 m in the second span in Model 2 and beginning diameter of the first span in 0.50 m and ending diameter of the second span is 0.40 m in Model 3. Diameters of the risers and ports are taken as 0.13 m in all models.Bell mouthed ports are also implemented in this study.Pipe thicknesses of the models are 0.01 m.

Figure 3 .
Three different fluid domains are created as diffuser models, based on the diameter of the diffuser sections.The dimensions of the domains are 25 m in the direction of diffuser, 1.5 m perpendicular to diffuser direction and 20 m in vertical direction (d) for all solid models.Domain sizes are determined in accordance with diffuser geometries and the conditions are citied from [23-27].The detailed explanations about these conditions are given by mentioned references.EOS material is utilized to model fluid with velocity of sound in salty water, co=1560 m/s and the constants (, 0) are equal to zero. is the slope of the Us−Up curve and 0 is the Grüneisen ratio.The properties of salty water are used to model internal and external flows at temperature of 20 C with mass density (ρf) of 1025 kg/m 3 and dynamic viscosity (μ) of 0.0015 Ns/m 2 .

Fig. 8
Fig. 8 Discharge velocities of Model 2 for unsteady external flow

Fig. 9
Fig. 9 Internal flow vectors of Model 2 for unsteady external flow Verification structural model using SDOF modelIn this section, Model 1 is selected to check on FSI results.For this purpose, mentioned model is modelled as SDOF model.Maximum displacements and modal behaviour of FSI model is compared with SDOF model.Equation of motion of SDOF model is given by (Eq.10).

Table 1
Structural results of FEM and SDOF for Model 1

Table 2
Average discharge velocities for unsteady external flow

Table 3
Average discharge velocities for steady external flow