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

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.


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.

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.
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: 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 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.
The relation (4) is very useful for solution of stationary problems, when velocity is just a function of coordinates (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.
The volume integral is simplified using (10) and (11) in surface integral, which is shown in (12).
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).
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.

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).  Table 1.  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]. • ). 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 Neo-Hookean model [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.

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).

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. Pressure outlet Velocity inlet

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).

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 Navier-Stokes 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. 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.  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. 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.

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  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/17-19444S is gratefully acknowledged for support of this work.