A new form of equation for force determination based on Navier-Stokes equations

* Corresponding author: [email protected]utbr.cz A new form of equation for force determination based on NavierStokes equations Simona Fialová1,*, František Pochylý1, and Dominik Šedivý1 1Brno University of Technology, Viktor Kaplan Department of Fluid Engineering, 61669 Brno, Czech Republic Abstract. This work is focused on calculating the force effects of an incompressible homogeneous liquid on a surface of a rigid or a flexible tube. An unsteady flow induced by differential pressure at the beginning and at the end of the tube is assumed. The pressure difference for the unsteady flow is determined experimentally. The mathematical model is based on modified Navier-Stokes equations. The unsteady term is modified in order to be able to use the Gauss-Ostrogradsky theorem to calculate the force. This method of solution will allow the calculation of the force by integration of the Navier-Stokes equations, which will help to refine and simplify the calculations. In the article, both methods of force calculation will be presented and compared both through the ANSYS FEA and CFD ANSYS Fluent solvers and by the integration of the Navier-Stokes equation. The calculation will not only respect the compliance of the tube but also its movement status. 1 Introduction The interaction between a flowing fluid and a structural phase has been extensively studied in the last decade. As a result of this interaction, eigenfrequencies of turbines blades, bridge pillars or pipelines can be impacted by the fluid. Outputs of these solutions can also be directional deformations, stress analyses, which are connected with force effects evaluations. The interaction of fluid with a solid phase is well-known as a fluid-structure interaction (FSI) [2]. FSI solutions of pliable tubes have most likely the biggest application in biomechanics. Materials of blood vessels, especially of arteries, belong to hyperelastic materials which means that they allow high elastic deformations. Many articles were written about this topic. FSI analysis of blood vessels was solved mostly as a numerical simulation. This approach enables the use of different material models, which cannot be used in experiment so easily. Most represented are models such as Money-Rivlin [1], Fung [3] and Neo-Hookean [3]. These constitutive models are used for description of hyperelastic, isotropic, homogeneous and incompressible materials. Another advantage of a numerical solution is the possibility to determine material stress from which force effects can be calculated. However, any numerical solution should be verified by experiments. Outputs from experiments in this area of research are usually velocity profiles, which are measured by particle image velocimetry (PIV) and radial deformation of pipe wall [4, 5]. Stress analysis can be done by tensiometers in the experimental part, but this measurement method is very complicated. This work presents a new form of Navier-Stokes equation, which can be used for force effect evaluations of unsteady flow in a flexible tube. A mathematical theory is supported by results of numerical solutions. 2 Mathematical theory The law of momentum conservation is in fluid mechanics of incompressible fluids described by Navier-Stokes equations (1) and mass conservation law is defined with continuity equation, which has the form of (2) for incompressible flow. 𝜌𝜕𝑣𝑖 𝜕𝑡+𝜌 𝜕 𝜕𝑥𝑗(𝑣𝑖𝑣𝑗)−𝜕𝜎𝑖𝑗 𝜕𝑥𝑗=𝜌𝑔𝑖 (1) 𝜕𝑣𝑖 𝜕𝑥𝑖=𝑑𝑖𝑣(𝑣)=0 (2) Where 𝜌 is density, 𝑣𝑖 and 𝑣𝑗 are velocity components, 𝑡 is time, 𝑥𝑖 and 𝑥𝑗 are components of coordinates, 𝜎𝑖𝑗 is stress tensor and 𝑔𝑖 is an external acceleration. A general equation of force on the wall is determined as a function of a stress tensor and it has the following definition: 𝐹𝑖=−∫ 𝜎𝑖𝑗𝑛𝑗𝑑𝛤 𝛤 (3) Where 𝐹𝑖 are force components (i = 1, 2, 3), 𝜎𝑖𝑗 is the stress tensor on the wall surface and 𝑛𝑗 is a unit external normal vector. In the case that the surface 𝛤 is very indented and also moving, numerical errors could come into existence during the solution. The stress tensor does not enable an analysis of the impact of the individual force components (inertial, pressure, viscous and so on). That is the reason © The Authors, published by EDP Sciences. This is an open access article distributed under the terms of the Creative Commons Attribution License 4.0 (http://creativecommons.org/licenses/by/4.0/). EPJ Web of Conferences 213, 02018 (2019) https://doi.org/10.1051/epjconf/201921302018 EFM 2018 why another relation was derived from the Navier-Stokes equations and the continuity equation. The area considered is shown in figure 1. The area is continuous with volume, labeled V, which is bounded by surface Γ. The surface can be rigid or pliable, stationary or moving. The inflow and outflow of fluid in volume V are on surfaces S. Surfaces of the volume are oriented by unit external normal vector. Fig. 1. Considered area Navier-Stokes equations (1) has in the considered area form (4). 𝐹𝑖=−[∫𝐴𝑑𝑉+∫(𝐵−𝐶)𝑛𝑗𝑑𝛺−∫𝜎𝑖𝑗𝑛𝑗𝑑𝑆 𝑠𝛺𝑉 ] (4) Where 𝐹𝑖 is defined by (3), 𝛺=𝑆∪𝛤 and other parts of (4) are defined by equations (5), (6), (7). 𝐴=𝜌𝜕𝑣𝑖 𝜕𝑡 (5) 𝐵=𝜌𝑣𝑖𝑣𝑗 (6) 𝐶=𝜌𝛿𝑖𝑗𝑔𝑘𝑥𝑘 (7) The relation (4) is very useful for solution of stationary problems, when velocity is just a function of coordinates (8). 𝒗=𝒗(𝒙)=> 𝜕𝑣𝑖 𝜕𝑡=0 (8) The volume integral of A has to be solved to determine the transient effects the force. This can be also affected by numerical errors. This problem can be easily solved for an incompressible fluid using the continuity equation (2). Equation (9) is derived from (2) by multiplying by 𝑥𝑗 and then it is integrated over the whole volume V. ∫𝜕𝑣𝑗 𝜕𝑥𝑗𝑥𝑖𝑑𝑉=∫𝑣𝑗𝑥𝑖𝑛𝑗𝑑𝛺 𝛺−∫𝑣𝑗𝜕𝑥𝑖 𝜕𝑥𝑗𝑑𝑉 𝑉= 𝑉0 (9) The volume integral is simplified using (10) and (11) in surface integral, which is shown in (12). 𝜕𝑥𝑖 𝜕𝑥𝑗=𝛿𝑖𝑗 (10) 𝑣𝑗𝛿𝑖𝑗=𝑣𝑖 (11) ∫𝑣𝑖𝑑𝑉 𝑉=∫𝑣𝑗𝑥𝑖𝑛𝑗𝑑𝛺 𝛺 (12) An application of Gauss-Ostrogradsky theorem has to be valid (13), (14). 𝑣𝑖=𝜕 𝜕𝑥𝑗(𝑣𝑗𝑥𝑖) (13) 𝜕𝑣𝑖 𝜕𝑡=𝜕 𝜕𝑥𝑗(𝜕𝑣𝑗 𝜕𝑡𝑥𝑖) (14) Then transient part (5) of equation (4) is modified to form (15). ∫𝐴𝑑𝑉 𝑉=∫𝜌 𝜕 𝜕𝑥𝑗(𝜕𝑣𝑗 𝜕𝑡𝑥𝑖)𝑑𝑉 𝑉= ∫𝜌 𝜕 𝜕𝑥𝑗(𝜕𝑣𝑗 𝜕𝑡𝑥𝑖)𝑑𝑉 𝑉= ∫𝜌𝜕𝑣𝑗 𝜕𝑡𝑥𝑖𝑛𝑗𝑑𝛺 𝑆∪𝛤 (15) It is obvious from (15), that the impact of local accelerations on the force, which is generated by interaction between fluid and wall, has its origin on the boundaries of solved area. This has a significant influence on force analysis. The new form of equation for force determination (16) is derived using (15) in (4). 𝐹𝑖=−[∫(𝐷+𝐵−𝐶)𝑛𝑗𝑑𝛺−∫𝜎𝑖𝑗𝑛𝑗𝑑𝑆 𝑠𝛺 ] (16) Where D is defined by (17). 𝐷=𝜌𝜕𝑣𝑗 𝜕𝑡𝑥𝑖 (17) It is possible from (16) to calculate the force components and execute an analysis of the influence of boundary conditions on force dependence in the solution of fluid-structure interaction. 3 MESH AND BOUNDARY CONDITIONS Two types of numerical solutions were done for comparison of the new type of force evaluation. Domains for numerical solutions were created in commercial software ANSYS Design Modeler. Two types of numerical solutions were solved. The first one simulated simply a flow through the fluid domain with rigid wall, so the mechanical parameters of the pipe were not included (CFD). The dimensions of the tube for CFD were the following: pipe length L = 500 mm; inner pipe radius Ri = 6.35 mm. The Fluid domain for FSI analysis (wall is not rigid but flexible) had the same dimension as in the first solution. The FSI solution included also a structural domain, which had the same length and inner radius, but it also had thickness in radial direction. The thickness was T = 1.6 mm, which means that the outer radius of the tube was Re = 7.95 mm. Both domains (CFD, FSI) were simplified for the sake of the FSI analysis, because it is a very time-consuming process. The cross-section of the pipe was not a circle, but rather a quarter circle (see figure 2). 2 EPJ Web of Conferences 213, 02018 (2019) https://doi.org/10.1051/epjconf/201921302018 EFM 2018 Fig. 2. Details of the meshes for fluid and structural domain The details of meshes are shown in figure 2. Both meshes were created in the commercial software ANSYS Meshing. The mesh element had the same length in axial direction for the fluid and the structural domain which was equal to 1 mm. 20 elements were used in tangential direction in both meshes. The fluid domain had 22 elements in radial direction. 3 elements were defined in radial direction in case of the structural domain. All elements were hexahedral. The mesh of the structural domains was composed of quadratic elements. The parameters of both meshes are stated in Table 1. Table 1. Parameters of meshes Domain Fluid Structure Number of elements 170 000 30 000 Number of nodes 186 873 155 727 Maximal aspect ratio 10.46 2.78 Maximal skewness 0.63 0.03 A decomposition of boundary conditions in computational domains is shown in Figure 3. Fig. 3. Boundary conditions A velocity boundary condition velocity inlet was applied at the inlet to the fluid domain. Pressure boundary condition pressure outlet was used at the outlet of the fluid domain. The velocity and pressure were specified by user defined function (UDF). The velocity was from the start of the solution to t = 0.5 s constant and its value was equal to 0.2 𝑚∙ 𝑠−1. The pressure was during this time linear increasing from 0 Pa to a value of 57 100 Pa. Afterwards, both boundary conditions (velocity and pressure) were defined by goniometric functions (see figure 4), which were obtained from an experiment in work [6]. Fig. 4. Shapes of boundary conditions Water was applied as a fluid (density 𝜌=998.2 𝑘𝑔∙ 𝑚−3 and dynamic viscosity 𝜇=1 𝑚𝑃𝑎∙𝑠). The inner wall of the pipe was prescribed as wall. The remaining two boundaries of the solved fluid volume were defined as symmetry. These fluid boundary conditions were the same for both types of the solution (CFD and FSI), but dynamic mesh was used in FSI analysis. The structural domain, which was solved only in the case of the FSI simulation, had fixed support at the pipe’s inlet and outlet boundary conditions, which allowed 0 degrees of freedom. The boundary conditions frictionless support were determined for symmetry simulation. Tygon was defined as the material of the tube. Its density was 1210 𝑘𝑔∙𝑚−3 and it was considered as a homogenous, isotropic and incompressible material. Tygon belongs to hyperelastic materials and it was specified by the NeoHookean model [3](18). 𝑊=𝐺 2(𝐼1−3) (18) Where 𝑊 [Pa] is strain energy density, 𝐺 [Pa] is an initial shear modulus, 𝐼1 [-] is the first invariant of the right Cauchy-Green deformation tensor. The parameters of the material were found out from a uniaxial tension test, which was done in work [6]. The initial shear modulus had value 1.3 MPa. 4 Numerical solutions Both computational simulations were solved with commercial software from ANSYS. CFD was done in ANSYS Fluent 17.2 and System Coupling was used for the need of the FSI analysis, which connected ANSYS Fluent 17.2 with Transient structural (FEM). 4.1 CFD settings Calculations were done as transient and the size of the time step was Δt = 0.005 s. 30 iterations were set for one time step. Scheme SIMPLE was used for the solution and schemes for space and time discretization were defined at second order precision. The simulation was solved as laminar for relatively low Reynolds numbers. The solution was calculated for 700 time steps. 0,202 0,204 0,206 0,208 0,21 0,212 0,214 0,216 0 20 40 60 80 100 0,5 1,5 2,5 3,5 v [m/s] p [kPa] t [s] Pressure outlet Velocity inlet 3 EPJ Web of Conferences 213, 02018 (2019) https://doi.org/10.1051/epjconf/201921302018 EFM 2018 4.2 FSI settings The coupling of FEM with CFD was realized in the environment of the ANSYS Workbench by software called System Coupling (SC). Simulations were solved as transient and the size of the time step was the same as in the case of CFD (i.e. Δt = 0.005 s). 5 iterations in CFD were specified for one iteration in SC. Scheme SIMPLE was used for the solution and schemes for space and time discretization were defined at second order precision. The simulation was solved as laminar for relatively low Reynolds numbers. The Dynamic Mesh, which is included in ANSYS Fluent, was applied for the sake of the wall deformation. The Dynamic Mesh was controlled with Smoothing which was specified by Diffusion method. The diffusion function was selected as cellvolume and its parameter was equal to 0. All fluid zones (velocity inlet, symmetry etc.) except the wall deforming were chosen as a deforming type. The wall deforming was defined as the system coupling type. A stabilization parameter was used for solution stabilization. Its scale factor was equal to the value of 2.5∙10−4 and the stabilization method was coefficient-based. Maximal values of residual for the convergence criterion were set to be equal to 10-5. The size of the time step in FEM was the same as in CFD. Large deflections were allowed and weak springs were turned off. Minimal count of coupled iteration in SC was equal to 5 and maximal value was 20. The duration of solution was specified by end time, which was set for 3.5 s. Maximal values of residual for convergence criterion were set as in CFD (10-5). 5 Results Three force components had an influence on the inner wall of the pipe. These components were directed in meaning of the coordinate system, which is depicted in Figure 3. The force effects were evaluated in the tube’s section from z = 0.2 m to z = 0.3 m. This decision was done for the sake of an effort to minimalize the influence of boundary conditions at the inlet and the outlet (fixed supports). Force components, which were oriented in the direction of x and y axes, had arisen because of the used boundary conditions (symmetry). Their time dependence for both solutions (CFD and FSI) is shown in Figure 5. Fig. 5. Time dependence of Fy Fx and Fy had the same values for the given solution, therefore were displayed just as the force components in the direction of y axis. These force effects were relatively big, but they would not exist, if it was a solution of a whole tube’s cross-section and not only a solution of a quarter circle. That is the reason, why these forces are not dealt with further in this work. Radial deformations are shown in Figure 6. The shapes of the deformations were similar to the shapes of the force components, which were affected in radial direction (in the case of FSI analysis). Fig. 6. Inner radius changes The determination of force, which was affected in axial direction (Fz), was first of all done for the CFD solution, because this solution had zero deformations of wall and therefore the evaluation was easier. Transient parts of (16) were solved by the central differences method. The comparison between the axial force determined straight from the CFD solution and the axial force evaluated by equation (16) is displayed in figure 7. Fig. 7. CFD axial forces The evaluation of force from the new form of NavierStokes equation is precise in case of the CFD solution. Peaks, which were around time 0.5 s, were caused by a numerical error during the change of boundary conditions progress at the inlet and outlet of the domain. An analogous determination was done also for the FSI analysis. However, in this case calculated with velocities on wall due to deformations of the wall. The comparison between the axial force determined straight from the FSI solution and the axial force evaluated by equation (16) is displayed in figure 8. Transient parts of (16) were again solved by the central differences method. 0 10 20 30 40 50 60 70 0 1 2 3 F [N] t [s] Fy - FSI Fy - CFD 6,2 6,4 6,6 6,8 7 7,2 0 1 2 3 R [mm] t [s] R - FSI R - CFD 0 0,0001 0,0002 0,0003 0,0004 0,0005 0 1 2 3 Fz [N] t [s] CFD 4 EPJ Web of Conferences 213, 02018 (2019) https://doi.org/10.1051/epjconf/201921302018 EFM 2018 Fig. 8. FSI axial forces The evaluation of the force from the new form of the Navier-Stokes equation is not precise as in the case of the CFD solution. Correctness of the solution was achieved by repetition of loading cycles. Higher peaks of the evaluated force were from the beginning of the solution caused by big volume changes (see figure 9). Several incorrectnesses, which were at the end of the solution, were probably caused by numerical errors in postprocessing. Numerical errors at the end of the solution could be eliminated using a smaller time step. Fig. 9. Volume and volume changes Fig. 10. Detail of axial forces An interesting fact was that the values of axial forces were very different between the CFD and FSI solutions (for comparison see figure 10). The differences between the forces were in the order of magnitudes. The axial force in the CFD solution had always a positive value, while the axial force in the FSI analysis acted also in the negative direction of axis z. These differences were probably caused by different velocity fields, which had arisen due to wall deformations. Velocity profiles were determined for both solutions for 10 different time moments in the xy plane for z = 0.25 m. Velocity profiles from the CFD solution are shown in Figure 11. The changes of the velocity profiles in CFD solutions were very small. This was caused by the velocity boundary condition. Velocities at the inlet to the domain were oscillating, but the amplitudes were very small, and therefore the velocities did not change much in the solid tube. Fig. 11. Velocity profiles - CFD Fig. 12. Velocity profiles - FSI Velocity profiles from the FSI analysis changed their shapes rapidly over time (see figure 12). It can be concluded that velocity profiles had definitely an influence on the different axial force in the case of FSI. The differences were also in velocity magnitudes, which affected magnitudes of the shear rate. The change of velocity magnitudes was probably caused by pulsatile motion of the wall of the pipe. A change of diameter in few time periods caused a backward flow in the vicinity of the wall. Reverse axial velocities had impact on the orientation of the axial force. 6 Conclusions The new form of Navier-Stokes equation was applied on force determination for two different types of numerical solution. The evaluation of force effects in case of the CFD solution was very precise. This was caused by simplicity of the solved problem which led to a very good convergence. The determination of the axial force in the -0,015 -0,010 -0,005 0,000 0,005 0,010 0,015 0 1 2 3 Fz [N] t [s] FSI Evaluation 0,0E+00 5,0E-06 1,0E-05 1,5E-05 2,0E-05 2,5E-05 -2,0E-05 -1,5E-05 -1,0E-05 -5,0E-06 0,0E+00 5,0E-06 1,0E-05 1,5E-05 2,0E-05 2,5E-05 0 2 V [m3] dV/dt [m3/s] t [s] dV/dt Volume -0,001 0 0,001 0,002 0,003 2,5 3 3,5 Fz [N] t [s] Fz - FSI Fz - CFD 0 0,05 0,1 0,15 0,2 0,25 0,3 0,35 00,005 vz [m/s] r [m] t = 3.5 s t = 3.375 s t = 3.25 s t = 3.125 s t = 3 s t = 2.875 s t = 2.75 s t = 2.625 s t = 2.5 s t = 2.375 s -0,2 -0,1 0 0,1 0,2 0,3 0,4 0,5 0 0,005 vz [m/s] r [m] t = 3.5 s t = 3.375 s t = 3.25 s t = 3.125 s t = 3 s t = 2.875 s t = 2.75 s t = 2.625 s t = 2.5 s t = 2.375 s 5 EPJ Web of Conferences 213, 02018 (2019) https://doi.org/10.1051/epjconf/201921302018 EFM 2018 case of the FSI analysis was more complicated. Integration of the transient part from the Navier-Stokes equation had to be done. Thanks to its new form, its integration did not have to be done in whole volume and it was done only at boundary conditions of this control volume. But the evaluation was still affected by numerical errors. Generally, the FSI analysis is more complicated than just a CFD solution. Results show that it is necessary to let run the FSI analysis for several loading cycles to get a precise solution. A smaller time step should also help with the correctness of the results. The outputs of both solutions showed the necessity of using FSI for specific problems. If the deformations of the wall were neglected for the simulation of flow in the Tygon tube, then the results would be very different from reality. The velocity profiles differed widely in case of the FSI analysis from the velocity profiles, which were determined from CFD solution. The same was valid also for the axial forces, which corresponded with shear rate near to the wall. General benefits of the new form of equation for force determination based on Navier-Stokes equations:  the new method converts the non-stationary part of the forces, induced by the local acceleration, from the volume integral to the spatial, which allows its control by boundary conditions, for example, in the optimization;  simplifies and refines the method of control volumes, since it converts the volume integral of local acceleration to spatial one;  for the non-stationary force calculation form the force definition it is necessary to determine the velocity gradients in the boundary layer which is considerably inaccurate - the new method is determined by velocity gradients at the inlet and outlet from the area where the gradient calculation is burdened only by a small error;  based on the new method, it is possible to perform a qualitative analysis of the influence of boundary conditions, as the possibility of nonstationary force effects can only be expressed at the border of the area;  this method can be applied also to a vorticity which non-stationary part can be evaluated only at the border of the area  the method can be extended to solution of nonstationary Maxwell equations Grant Agency of Czech Republic, within the project GA101/1719444S is gratefully acknowledged for support of this work. References 1. Paidoussis, M. Fluid-Structure Interactions. Edtion ed. London: Academic Press, 2004. 2. Kim, Y. H., J. E. Kim, Y. Ito, A. M. Shih, et al. Hemodynamic analysis of a compliant femoral artery bifurcation model using a fluid-structure interaction framework. Ann Biomed Eng, Nov 2008, 36(11), 1753-1763. 3. Torii, R., M. Oshima, T. Kobayashi, K. Takagi, et al. Influence of wall thickness on fluid-structure interaction computations of cerebral aneurysms. International Journal for Numerical Methods in Biomedical Engineering, Mar-Apr 2010, 26(3-4), 336-347. 4. Bazilevs, Y., M. C. Hsu, Y. Zhang, W. Wang, et al. A fully-coupled fluid-structure interaction simulation of cerebral aneurysms. Computational Mechanics, Jun 2010, 46(1), 3-16. 5. Pielhop, K., M. Klaas AND W. Schroder Experimental analysis of the fluid-structure interaction in finite-length straight elastic vessels. European Journal of Mechanics B-Fluids, Mar-Apr 2015, 50, 71-88. 6. Pielhop, K., C. Schmidt, S. Zholtovski, M. Klaas, et al. Experimental investigation of the fluid-structure interaction in an elastic 180 degrees curved vessel at laminar oscillating flow. Experiments in Fluids, Oct 2014, 55(10), 13. 7. Šedivý, D., J. Burša, S. Fialová. IAHR 2018 (to be published) 6 EPJ Web of Conferences 213, 02018 (2019) https://doi.org/10.1051/epjconf/201921302018 EFM 2018