SIMULATION OF FLUID-SOLID INTERACTION USING THE MESH-ADAPTIVE SUBGRID SCALED TURBULENCE MODEL

Simulation of Fluid-Solid Interaction Using the Mesh-Adaptive Subgrid Scaled Turbulence Model

Ulises F. Gonzalez
John W. Holtz
ALGOR, Inc., 150 Beta Drive, Pittsburgh, PA 15238, USA www.algor.com jhalapch@algor.com

ABSTRACT

The primary aim of the study presented in this paper is to apply the principles of physics-based virtual reality to simulate the interaction of fluid flows with solid structures. We utilize a Computational Fluid Dynamics (CFD) method that employs Finite Element Analysis (FEA) and an algebraic eddy viscosity model to account for turbulence [Peyret and Taylor (1983), Deardorff (1970)]. The method accounts for turbulence of the length-scale of the local mesh size. There is no other input parameter directly associated with the modeling of turbulence. The user simply uses a detailed mesh wherever in the model fine-scale turbulence is expected. Such regions are generally associated with obstructions, and are thus easy to identify. Coarser meshes can be used in portions of the model where little or no turbulence occurs. This points out another advantage of the method presented in this paper: it can be used to easily simulate flows in which several turbulent length scales coexist.

The results of the CFD analysis are integrated with Algor's Mechanical Event Simulation (MES) software to calculate the effect that the flow has on solid structures. Of course, the results of the MES can then be input into the CFD analysis. This process is continued until the solution is indeed a virtual simulation of the event. This interaction between fluid flow and solid stress analysis has been generally regarded as the realm of advanced users. Generally, only such users were adept at avoiding the instabilities commonly present in either analysis. Both the CFD method and MES have built in intelligence to avoid these instabilities. The elimination of these instabilities combined with Algor's mesh-based turbulence model allows typical designers, not just experts, to simulate the event associated with the fluid-solid interaction as if it were a physical experiment.

For instance, the user can observe the motion of the solid as it interacts with the fluid. Furthermore, the mesh size in the fluid domain is analogous to the level of detail of a flow visualization technique.

NOMENCLATURE

	D₁ -	Diameter of duct away from orifice
	D₂ -	Diameter of duct at the orifice
	E -	Young's Modulus
	L -	Measure of the local mesh size
	L₁ -	Length of flow domain upstream of orifice
	L₂ -	Length of flow domain downstream of orifice
	Q -	Volumetric flow-rate
	Re -	Reynolds number
	S_ij -	Rate of strain tensor
	t -	Thickness of orifice
	U -	Velocity vector distribution
	X -	Spatial vector coordinate
	t -	Time-step size
	-	Dynamic viscosity
	-	Mass density

INTRODUCTION

In this paper, we consider the multi-physics problem of fluid-solid interaction. We are particularly interested in how turbulent flows interact with solids. Specifically, how do the stresses caused by the turbulent flow deform the solid, and how does such deformation in turn affect the flow. Algor's goal is to incorporate physics-based virtual reality to simulate the interaction of fluid flows with solid structures, including the key element of providing a feedback loop between the flow and the deformation.

The CFD method that we utilize uses FEA. The finite element mesh does not have to be structured. Some CFD schemes rely on the use of structured meshes that may be difficult to obtain for many of today's CAD-generated complex geometries. The characteristics of the mesh are highly intertwined with the means used to model turbulence. We utilize an algebraic eddy viscosity turbulence model in contrast to the more commonplace closure models [Moin and Kim (1982), Smagorinsky (1963)]. Closure models vary in complexity from the so-called one-constant models to the two-equation models that include the widely used k-epsilon model [Peyret and Taylor (1983)]. All of these models assume that the turbulent fluctuations about a mean flow are represented by a Reynolds stress field [Peyret and Taylor (1983)]. This additional stress results in additional terms to the Navier-Stokes equations (the basis of most CFD codes). These additional terms are calculated using differential equations particular to the chosen closure model. These equations include empirical constants that are only known for certain flow configurations. Boundary layers are defined in order to adapt the closure model to a given flow configuration. Significant experience is usually required to establish the shape of these boundary layers and to "tune" the empirical constants.

Even though closure models make gross approximations and rely on empirical data, they have yielded acceptable results for coarsely meshed models [Peyret and Taylor (1983)]. In the past, when computing resources were at a premium, these models provided the most practical means to account for turbulence in CFD codes. Advances in computing power have made it possible to consider much finer meshes. Of course, the eventual goal has to be to use a mesh fine enough so that no turbulence model is required. It is well recognized that the Navier-Stokes equations can be used to describe the smallest scales of turbulence as long as the mesh size is small enough [Peyret and Taylor (1983)]. The current reality is that computing power is still well away from the eventual goal of no turbulence modeling. But, computing power is currently well within the requirements of algebraic eddy viscosity models. In contrast to closure models, algebraic eddy viscosity models do not utilize differential equations to yield the Reynolds stress field - the additional stress caused by turbulence.

The CFD method presented in this paper uses an algebraic eddy viscosity model. When using this type of turbulence model, CFD methods utilize the Reynolds-averaged Navier-Stokes equations. These equations are identical to their non-averaged counterparts, but recognize the high cost of computing small-scale turbulence. Specifically, a non-constant, so called turbulent or eddy, viscosity is used for flow scales below a selected cutoff length. This cutoff length is related to local mesh size. Thus, the smaller the mesh size, the smaller the detail of the turbulence that can be captured. As the mesh size becomes finer the method approaches the solution of the classical Navier-Stokes equations with no special handling of turbulence. It should be noted that the method yields the best possible solution for a given mesh size. In contrast to closure models, there is no need to establish boundary layers or to estimate constants.

Fine meshes should be used in regions where small-scale turbulence is expected to occur. Generally, such turbulence occurs near easily identifiable obstructions to the flow. The user simply reduces the mesh size near each obstruction, and thus allows the method to "capture" turbulent flow scales above this local mesh size. This makes it possible to easily model flows in which several turbulent length scales coexist by simply using different mesh sizes throughout a model. In a sense this method provides the user with a "What You See is What You Get" model of turbulence.

The calculation of the turbulent flow field is only part of the feedback loop required to solve the fluid-solid interaction problem. Care must be taken to prevent instabilities incurred when solving (dynamic) fluid flow and solid stress analysis. Further care has to be taken to ensure that the coupling between the analysis (the feedback loop) is also stable and accurately describes the physics of the event. The work presented in this paper addresses all of these stability and accuracy concerns. Thus, this work provides the foundation for physics-based virtual reality for the problem of fluid-solid interaction.

The solid mechanics aspect of the fluid-solid interaction problem is accomplished using Mechanical Event Simulation (MES). The MES method is geared towards the simulation of mechanical parts that experience motion. Such is the case for many problems involving the interaction of fluids with solids. MES is based on nonlinear FEA, which has generally been regarded as both highly unstable as well as computationally intensive. The instabilities can usually be eliminated using a smaller time-step size, t. Of course, the smaller t, the larger the computing time. This trade-off between stability and efficiency can be significantly diminished by incorporating an "intelligent" (automatic) algorithm to change t throughout the analysis. MES utilizes such an algorithm; t is automatically reduced during critical times and again increased if/when conditions no longer warrant such temporal detail. Thus, this "intelligent" algorithm can actually improve computational efficiency.

METHOD

In this paper we focus on three aspects of Algor's approach to solving problems involving fluid-solid interaction: (1) addressing stability and accuracy concerns, (2) the algebraic viscosity turbulence model, and (3) meshing techniques.

When attempting to solve coupled analysis not only must care be taken to guarantee that each analysis is stable and accurate, but special treatment is usually necessary to ensure that the link between the (individual) analysis is physically valid and that it not cause unnecessary instabilities. In this paper, we are primarily concerned with the dynamic coupling of fluid flow and solid stress analysis. Each respective analysis contains an automatic time-stepping scheme. In either case, the time-step size is automatically reduced at critical points in the solution. Specifically, the time-step size is reduced whenever it is not small enough to adequately capture pertinent underlying physical details. Of course, the time-step size is increased if the critical phase of the event were to end. This ability to automatically adjust time-step sizes provides for highly stable, computationally efficient solution methods. The stability and efficiency of each separate analysis significantly contributes to the overall stability and efficiency of the coupled system. Nevertheless, the link between analyses must be highly scrutinized for accuracy and stability losses. In the case of fluid-solid interaction, special attention must be given to the changing shape of the fluid domain. As in all coupled systems, the possibility that each analysis may require a different time-step size must be considered.

It should be mentioned that the interaction between the CFD solver and MES is explicit. That is, each analysis is run using the results of other analysis as its input. It is not necessary for each analysis to have equal time-step sizes, but the analyses must interchange data at equivalent times. Specifically, M number of MES time-steps and N number of CFD time-steps can occur between the times when data is exchanged. It is possible to increase the accuracy of the coupled system by utilizing a "predictor-corrector" scheme during each dataexchange cycle. A cycle is comprised of the following steps: (1) Resultant fluid pressures at the surface of the solid are used by the MES method to calculate how the solid deforms; (2) This new deformed shape is used as the boundary of a new (generated) mesh for the fluid domain. Care is taken to ensure that nodal locations match at the boundary and (3) The CFD method uses this new domain to obtain fluid pressures. Note how there is no need to remesh the solid domain. Finally, the accuracy of the coupled system should be considered by examining whether solutions that utilize different time-step sizes and mesh densities yield similar and convergent results.

As we mentioned in the prior section, we utilized an algebraic eddy viscosity model to describe the consequence of turbulence. This model is incorporated into Algor's CFD code and is termed as the "mesh adaptive subgrid scaled turbulence model." Our model, just like other algebraic eddy viscosity models [Moin and Kim (1982), Schumman (1975), Smagorinsky (1963)], uses a turbulent or eddy viscosity proportional to the quantity ,

where X is the spatial coordinate, U is the velocity distribution, L is a measure of the local mesh size, and

where S_ij is the rate of strain tensor,

It should be noted that f(U) is a mathematically concise expression for the shear (or distortion) rate of the flow. This time and space dependent eddy viscosity is an integral part of the Reynolds-averaged Navier-Stokes equations. These equations are solved for the velocity and pressure distributions. This solution is capable of accurately describing turbulent scales (i.e., eddy sizes) comparable to the local mesh size without having to directly model the filler details that govern the dynamics at these length scales.

Sophisticated meshing techniques are used in order to obtain the variable mesh sizes required to capture different turbulent length scales throughout a given model. It would not be computationally efficient to use a fine mesh in regions where (relatively) small-scale turbulence does not occur. The meshing scheme should allow the user to easily specify distinct mesh sizes throughout the model. Meshes for 2-D models can easily be constructed to such user-specifications using our planar meshing algorithm. This algorithm generates a mesh composed solely of quadrilaterals - avoiding any triangles. Triangular elements generally result in inaccuracies when solving CFD problems. For 3-D situations, we consider an automatic solid meshing scheme that utilizes the local surface mesh size to generate the solid elements in the nearby region. Thus, the user need only manipulate the mesh density of the surface mesh to obtain the desired inhomogeneous mesh. In order to avoid the inaccuracies encountered when solving CFD problems using degenerate elements, Algor's automatic solid meshing scheme is capable of generating a mesh composed solely of 8-node brick elements. It should be noted that Algor's CFD FEA method does not rely on the use of structured meshes. Thus, neither the 2-D nor the 3-D meshers need to generate structured meshes. Other CFD schemes rely on the use of structured meshes that may be difficult to obtain for many of today's CAD generated complex geometries.

RESULTS

In this section, we simulate the flow of air through a circular duct with an orifice (see Fig. 1). Such a geometric configuration can be used to demonstrate the power of Algor's method for solving turbulent fluid-solid interaction problems. For sufficiently large Reynolds numbers, the flow will be turbulent with small-scale turbulence caused by the orifice. The flow can induce deformations of the confining solid. In this paper we consider the case of a rigid duct wall, but allow for the obstruction that creates the orifice to deform. Of course, once the obstruction deforms, the flow itself is affected. Thus, we have a situation in which there is a feedback loop between turbulent flow and a solid structure.

Figure 1: Sketch of circular duct with orifice.

We first concentrate on the CFD aspect of the problem. The diameter of the flow domain, D₁, is 0.25m. The diameter at the vena contracta (the orifice), D₂, is 0.165 m. The computational domain extends a distance, L₁, of 2.5 m upstream of the orifice to allow effects caused by the inlet to dissipate by the time the flow reaches the obstruction. The domain also extends a distance, L₂, of 3.0 m downstream of the orifice to ensure that the simulation captures all significant flow disturbances caused by the obstruction. The thickness of the orifice, t, is 0.005 m. Air is taken to have a density,

, of 1.20 kg/m³, and a viscosity,

, of 1.80X 10^-5 N-s/m². The flow-rate, Q, at the inlet to a value of 1.0 m³/s. This flow configuration is turbulent with a Reynolds number

In order to properly simulate the disturbance caused by the obstruction, it is necessary to use a fine mesh near the orifice. In regions away from the orifice it is appropriate to use a much coarser mesh. We utilized our 3-D mesher to generate such an inhomogeneous mesh. In this phase of the project we did not consider the interaction of the flow with the solid. That is, we took the mesh of the fluid domain as constant. Figure 2 shows contours of the steady-state velocity field along a length-wise cross-section of the duct near the orifice. It should be noted that for this case, the velocity field is symmetric about the axis of the duct. Algor's 3-D method could have simulated a non-symmetric field had secondary flows occurred. The accuracy of our solution was verified by comparing the head loss generated by the orifice under state-state conditions. Figure 3 shows the pressure contours near the orifice. The pressure drop is measured between cross-sections of negligible radial gradient in pressure. The head loss of 1676 N/m² predicted by our solution compares favorably with the 1650 N/m² predicted by the analytical solution of Fox and McDonald (1985). This analytical solution accounts for turbulence effects by using a correlating equation.

Figure 2: Level contours of the steady-state velocity (m/s) field along a length-wise cross-section of the duct near the orifice. The other half of the duct is not shown for the sake of visualization. The maximum velocity (not shown) was 73.5 m/s. These results are from a simulation that considered the orifice to be rigid.

Figure 3: Level contours of the steady-state pressure (N/m²) field along a length-wise (half) cross-section of the duct near the orifice. The top half of the duct is not shown for the sake of visualization. These results are from a simulation that considered the orifice to be rigid.

The next step in the project presented in this paper was to allow for the flow and the solid to interact. We assumed that the duct's walls were rigid, but we did allow for the constriction to deform. The solid stress analysis was performed using Algor's Mechanical Event Simulation (MES) with Nonlinear Material Models software. This software is capable of simulating geometric large deformations. We considered two possible materials for this constriction, steel and rubber. The steel was taken as having a Young's Modulus, E, of 2.0 X 10¹¹ N/m² while the rubber had an E of 3.4 X 10⁶ N/m².

Figure 4: Side view of steel orifice (shaded rectangle). Primary flow direction is from left to right. Velocity field is depicted using small point vectors.

Figures 4 and 5 show the combined fluid flow and solid stress steady-state solutions for the steel and rubber orifices respectively. It is apparent from Fig. 4 that the steel orifice deflected an insignificant amount (7.5 X 10^-7 m) whereas, Fig. 5 demonstrates that the rubber orifice deflected significantly (1.6 X 10^-2 m). Because of the relatively large deformation of the rubber orifice, the solutions for its corresponding fluid domain were obtained using a different mesh at each time-step.

Figure 5: Side view of rubber orifice (shaded deformed rectangle). Primary flow direction is from left to right. Velocity field is depicted using small point vectors.

From Figures 4 and 5 one can observe how the flow is significantly affected by the deflection, or lack thereof, of the orifice. Both orifices generate similar turbulent eddies downstream of the orifice. Upstream of the orifice, the basically rigid steel orifice yields a noticeable eddy; whereas, the highly deformed rubber orifice results in a relatively smooth flow. The smoother flow in the case of the deflected orifice is expected. Generally, in situations when the flow and the solid significantly affect each other, the combined system usually attains an equilibrium that minimizes instabilities (e.g., turbulent eddies).

It should be pointed out that we were not guaranteed that a steady-state solution could be achieved. If secondary flows had been present, then these flows might have resulted in the periodic shedding of small eddies from the orifice. These eddies would most likely be shed in a non-axisymmetric manner. It was for this reason that we considered a 3-D rather than a 2-D axisymmetric model. Such periodic shedding of eddies is capable of producing periodic deflections of the flexible structure, thus preventing the ability to reach a steady-state solution. It is interesting to note that it was eddies from a secondary flow that resulted in the resonant periodic forcing and ultimate destruction of the first Tacoma Narrows Bridge in Washington state (USA) on November 7, 1940.

CONCLUSION

In this paper we presented a method for the simulation of events involving fluid-solid interaction. The solution of the fluid flow accounted for turbulent behavior using an algebraic eddy viscosity model called the mesh adaptive subgrid scaled turbulence model. This model does not require the expertise of closure models. Specifically, there is no need to define boundary layers nor to estimate flow-specific constants. Instead, this model automatically accounts for turbulence below the local mesh size. Thus, to observe a given turbulence scale (eddy size), the user need only use a similar-sized mesh. It is only necessary to use inhomogeneous meshes when attempting to simulate several turbulent length scales. Our method includes the (automatic) generation of such meshes with minimal user effort.

Both the CFD and the solid stress (MES) FEA processors were made highly stable by the use of "intelligent" automatic time-stepping schemes. Stability and accuracy concerns governed the method used to couple these processors. This coupled - system was used to successfully simulate the turbulent flow of air in a duct with a flexible orifice. Such a simulation verifies the capability of this combined system to apply the principles of physics-based virtual reality to the problem of fluid-solid interaction. It should be mentioned that this combined system can also be expanded to account for heat transfer effects. This grander system can be used to simulate flows as complex as those within a typical personal computer.

REFERENCES

Deardorff, J. W., (1970) "A Numerical Study of 3-D Turbulent Channel Flow at Large Reynolds Numbers," J. Fluid Mech. Vol. 41, pp. 453-480.

Fox, W., and McDonald, A., (1985) "Introduction to Fluid Mechanics," John Wiley & Sons, Inc., New York, N.Y.

Moin, P., and Kim, J., (1982) "Large Eddy Simulation of Turbulent Channel Flow," J. Fluid Mech. Vol. 118, pp. 341- 365.

Peyret, R., and Taylor, T. D, (1983) "Turbulent-Flow Models and Calculations," Springer-Verlag, Inc., New York, N.Y.

Schumman, V., (1975) "Subgrid Scale Model for Finite Difference Simulations of Turbulent Flows in Plane Channel and Annuli," J. Comp. Phys., Vol. 18, pp. 376-404.

Smagorinsky, J. S., (1963) "General Circulation Model of the Atmosphere," Mon. Weather Rev., Vol. 91, pp. 99-164.