Abstract
A new virtual baffle methodology is implemented to solve contact/detach problem which is often encountered in fluid and structure interaction simulations while using dynamic grids technique. The algorithm is based on tetrahedral unstructured grid, and a zero thickness baffle face is generated between actually contacted two objects. In computation process, this baffle face is divided into two parts representing convective and blocked area, respectively; the area of each part is calculated according to the actual displacement between the two objects. Convective part in a baffle face is treated as inner interface between cells, and on blocked part wall boundary condition is applied; so convective and blocking effect can be achieved on a single baffle face. This methodology can simulate real detaching process starting from contact, that is, zero displacement, while it has no restriction to minimum grid cell size. The methodology is then applied in modeling of a complicated safety valve opening process, involving multidisciplinary fluid and structure interaction and dynamic grids. The results agree well with experimental data, which proves that the virtual baffle method is successful.
1. Introduction
Numerical simulation of fluid and structure interaction problem or transient flow problem involving multibody in relative motion using dynamic grid technique often encounters difficulties of dealing with contact/detach problem. Some engineering examples are known as canopy ejection from a fighter plane, cowl ejection of rocket, door opening of internal weapon bay of a stealth fighter, valve washer opening from shut, and so forth. Common feature of these fluid and structure interaction problems is that two objects contact originally and then depart from each other; that is, the distance between them increases gradually from zero. In the process of detaching from contact, the computational domain changes its topology, which causes dynamic grid operation less robust and much more time consuming. In fact, there is almost no accurate dynamic grid method that can describe such topology change faithfully. So, dealing with contact/detach problem has been becoming a focus point of dynamic grid technology and a bottle-neck to particular fluid and structure interaction simulation. Therefore some competent method should be developed to meet the more and more eagerly demand from engineering analysis field for solving this problem. Recently state of the art of the method for this problem is presetting a small distance between the two objects to avoid computational domain topology change, which means that the actual initial position is altered. Sharov et al. [1] simulated the canopy ejection process of F/A-18 fighter, and the computation started at a position corresponding to 26 milliseconds after the canopy started to move in reality. Method like this actually does not resolve contact/detach problem but just partially simulates the moving procedure from midway. So the smaller the initial distance between the moving objectiss, the closer will the result be to real case. But, very small gap between objects means accordingly small grid cell size in the gap area, which will affect grid points distribution and restrict computational time step and then affect the whole computational efficiency consequently. Furthermore, flow field of fluid and structure interaction is usually highly unsteady and nonlinear, and computations often begin with an initial steady-state flow field. In the canopy ejection case, for example, computation should start from steady-state flow with the canopy tightly closed. Therefore, computing from midway will start from a wrong initial condition, that a flow has already been established inside the cockpit while in reality there is no flow at all. This is not acceptable in some serious application such as flow field in a valve when the washer opens from shut position. To overcome these difficulties, a boundary condition variable, zero thickness, and virtual baffle methodology is implemented in this article. The main idea is that a zero thickness baffle face is generated between actually contacted two objects at grid generation stage, and, to simulate the initial contacting status, wall boundary condition is applied on two sides of the baffle face; then, when objects depart from each other, the baffle face is divided into two parts representing convective and blocked area, respectively, and the area of each part is calculated according to the actual displacement between the two objects. The integration of convective flux through the element of baffle face is also treated with combination of convective and blocking effect. This new method can simulate real detaching process starting from contact, that is, zero displacement, while it has no restriction to minimum grid cell size.
2. Governing Equations and Numerical Discretization
2.1. Governing Equations
The three-dimensional time-dependent Euler equations under ALE (Arbitrary Lagrangian-Eulerian) description can be expressed in the integral form for a bounded domain with a boundary as where denotes unit outward normal vector to the boundary. The conservative variable vector , and inviscid convective flux vector are defined by where , , and e denote the density, pressure, and specific total energy, respectively, denotes velocity of fluid particles relative to the motion of the grids, denotes velocity of grid motion, and , , are components of in x, y, z direction, respectively. and are given by Here , , and are components of grid velocity in x, y, and z direction.
In many application cases, like air flow within a valve, a large range of flow speed exists, from nearly static to supersonic. To make the flow solver capture all the features of the whole flow field, we would like to employ a preconditioned equation system. The derivation of the preconditioned governing equations begins by transforming the original system of equations from the conservative variables to the primitive variables q as follows: where , and the following preconditioning matrix , which is modified from Jacobian , and originally proposed by Choi and Merkle [2] and extended further by Weiss and Smith [3], is used in the present work: where H is total enthalpy per unit mass, . is the modified , given by ; here is a reference velocity. The choice of the reference velocity is crucial for success of the preconditioning method. If the magnitude of the reference velocity is equal to the local speed of sound, the preconditioned system (2.4) recovers to the nonpreconditioned system (2.1). In order to make all eigenvalues of preconditioned system (2.4) have the same order of magnitude, the reference velocity must be of the same order as a local velocity. In this work, the following reference velocity is used: where is the local velocity magnitude, c is the speed of sound, and is a constant. In the present work, is fixed at 0.5.
2.2. Discretization Scheme
Integrating (2.4) within control volume (3D tetrahedral grid cell) yields, where is the area vector of the kth face of a tetrahedron cell, and the summation exerts over all the 4 faces. In order to obtain a second-order approximation of the convective flux on integral face, the primitive variable is reconstructed at the centroid of the face through Taylor expansion from cell centroid. The reconstructed values on a face from two sides yield a Riemann problem, which is solved by Flux-Vector-Splitting method, and Van-Leer splitting scheme is used in the present work: Reconstructed primitive variables are calculated by where is flux limiter, which will function as a smoother across discontinuities in supersonic region. The gradient of variables is obtained through Green’s formula, given by Equation (2.10) is integrated within 3D grid cell: The variable location and relationship with grid cells are illustrated in Figure 1.
The flux limiter used here is following the work of Barth and Jespersen [4], which is given by Here , , and subscript “j” represents the 4 neighbors of cell i. is unlimited reconstructed primitive variables on face k of cell i, and given by The limiter for cell i is set to be the minimum value of limiters on all of its faces, and for an integral face, the final limiter used is the minimum value of limiters of the two cells sharing this face:
The time integral scheme of (2.7) used in the present work is standard four-stage Runge-Kutta scheme, and for the calculation of volume of deforming grid cells, Geometric Conservation Law (GCL) is satisfied numerically.
2.3. Integral of Convective Flux on Virtual Baffle Face
Generally, (2.7) is applied to ordinary grid cells in computational domain; the summation exerts over all the 4 triangular faces. If one face of a cell locates on virtual baffle face, say for example, face of cell i, shown in Figure 2, then (2.7) should be modified when applied to cell i as The right-hand side of the above equation consists of two terms, one is the summation of flux integral on three ordinary faces, just like what appears in (2.7), and the other is the flux integral on the face locating on virtual baffle. Here and denote the area and the convective flux vector at the centroid of the face locating on virtual baffle, respectively. Notice that the value of a variable at the centroid of a triangle is a second order approximation of the average value of the same variable distributed on this triangle; so solution at centroid can be regarded as average value within a second-order discretization frame. Figure 2 is a sketch of virtual baffle face, showing the location of face . It can be seen from Figure 2 that the preset virtual baffle has the width of two layers of cell face, but actually it can contain any number of layers of cell face subject to requirement of certain application case. A gap marked with shadow represents the actual distance between two contacted objects, therefore, the shadow part in triangular face should be treated as convective, and the rest part should be treated as blocked. So can be defined by where and denote the area of convective part and blocked part, respectively. The convective flux on virtual baffle face is calculated by Here is the flux vector at the centroid of face calculated assuming that the face is fully convective, that is, calculated by (2.8), and is the flux vector at the centroid of face calculated assuming that the face is fully wall boundary, and is given by For inviscid flow, where no-flow-through slipping boundary condition is exerted on wall boundaries, we get ; thus (2.18) can be simplified as What should be noticed in addition is that the solution gradient in a cell with a face locating on virtual baffle, should also be computed in two manners regarding that the virtual baffle face is fully convective or fully blocked, respectively, and then the corresponding gradient is used to calculate the flux and .
3. Dynamic Grid Algorithm
For applications with geometric deformation not too large, one of the most robust and widely used dynamic unstructured grid methods is spring analogy and originally proposed by Batina [5]. The unstructured mesh is considered as a system of interconnected springs, and each edge of any tetrahedron cell is represented by a tension spring. Various attempts at determining the optimum relationship for specifying the spring stiffness have been made by many researchers [5, 6]. Following the research of Blom [6], in the present work, the spring stiffness is basically assumed inversely proportional to the length of its edge and then imposed boundary modification and torsional effect modification and can be written as where is boundary enhancement factor, which is a function of the distance between this edge and wall boundaries, and its effect is to stiffen the spring stiffness near object surface. represents torsional effect, whose value is set to the minimum dihedral angle facing this edge in all of the grid cells sharing this edge. The purpose of is to count for the torsional restriction along with the linear tension restriction of the cell edges. is the length of edge connecting nodes i and j. In the present work, set , for the nearest one layer of grid cells, and , for other grid cells. Assuming the total force exerted on a node by springs sharing this node will keep constant, that is, the total force exerted on a node after grids deforming is the same as that in initial status, which yields Here is the coordinate vector of node i, and Ni is the number of nodes connected to node i. Applying (3.2) to all nodes will construct a linear system , and is coefficient matrix formed by stiffness of springs and is diagonal dominant. The coordinates of grid nodes can be updated by solving the above linear system through iteration method. Because the solution on both the fixed boundaries and the moving boundaries are known for the next time step, the boundary condition for the above linear system is Dirichlet type, and acceptable result will be obtained through 3 or 4 Jacobi iterations. Details about spring analogy can be found in [7].
From (3.1), we can see that factors and introduce dependence on the coordinates of nodes to the stiffness matrix, which will result in a slightly nonlinear system. To avoid solution diverging, update of (3.1) should not be very frequent. In the present work, (3.1) is updated every 500 time steps.
4. Application in Fluid and Structure Interaction Simulation of a Safety Valve
4.1. Description about the Valve
Now the whole methodology is applied in the numerical simulation of the opening process of a safety valve. The valve considered here is a device used to keep the air pressure constant within a container. It can be opened or shut off by itself driven by the air flow within it, is widely used in a large variety of industry field, and is playing an important role to adjust the working condition of all kinds of aerodynamic facilities. Figure 3 illustrates the working principle of the valve considered in this article. The whole valve is made up of two main components, namely, main valve and command valve, which have similar internal structure including spring, diaphragm, washer, and push rod. The main valve is divided by the main diaphragm into two compartments, room A and room B. When the air pressure is low enough in the container which is connected with the valve through the inlet, room A is linked with room B via command valve and feedback pipe; so the air pressure in room A and room B, that is, on the two sides of the main diaphragm, is nearly equal; therefore, the main washer will be at the shut position under the pressure of room A and the pressing force from main spring. The shut state is shown in Figure 3(a). With the pressure increasing in container and room A, the pressure in command valve and room B will increase simultaneously, and at a moment, when the pressure becomes high enough to press the command valve diaphragm back against the spring in command valve, the connectivity of room B with room A will be shut off by command valve mechanism, and then, the connectivity of room B with outlet is established, and Room B will deflate. With pressure decreasing in room B, the pressure difference on the two sides of main diaphragm becomes larger and larger, and finally the total force will be large enough to press the main diaphragm and main washer back against the main spring, which makes the main washer open and room A and container deflate. The opened state is shown in Figure 3(b). With the pressure in room A decreasing to a low enough value, a reverse process will occur, and the main washer will finally recover to shut position.
(a)
(b)
Numerical modeling of the opening process of such safety valve must take the challenge of fluid and structure coupling, and the most crucial thing is to achieve a proper treatment of the washer moving from tightly shut position. Recently, a lot of research work about flow field in a valve using numerical analysis method has been done by some researchers [8–14], some of them simulated transient flow properties without coupling the motion of structure, and some used dynamic grid technique and prescribed structure motion, and others only solved steady-state problem or “water-hammer” problem using characteristic method. Almost no similar research work covering contact/detach problem has been found.
4.2. Computation of Solid Structure Movement
There are two similar sets of moving components in main valve and command valve, respectively; both of them can be treated as 1D mass-spring-damp vibration system. The governing equations of the motion of moving components can be derived from Lagrangian equation. Here only the governing equation for the main valve is given for simplicity; equation for command valve is similar: where denotes the mass of main spring, denotes the effective mass of the main diaphragm, and is the mass of the rest part of the whole moving component assembly including main washer, central rod, and pressing plate, is the stiffness of the main spring, is the stiffness of the plastic sealing layer on the outer side of main washer, and is the stiffness of central rod when colliding with the case in maximum range of opening. or is taken into account only when collision occurs and will be ignored at any position except the two ends of opening range. The values of and are determined in advance by finite element analysis. is the precompressed length of the main spring, and is the length of the total opening range. represents elastic force of the main diaphragm when deformed and is a function of the moving displacement x. The relation between and x, regarded as the stiffness of the main diaphragm, should also be determined in advance. Considering that the thickness of the diaphragm is very small, and the displacement might be relatively large, in the calculation of , elastic-plastic model is used and geometric nonlinearity is considered. The function value of is calculated in discretized pattern and is stored. At any given x, a force value can be obtained by interpolation. In the present work, the motion of the main diaphragm is regarded as a combination of translation of central part and deformation of outer rim part. The central part is that covered by pressing plate. Correspondingly, the aerodynamic force exerted on the whole moving component assembly is also divided into two parts, represents the aerodynamic force on the outer rim part of the diaphragm, and represents the aerodynamic force on the rest part of the assembly. is a factor used to take account of the effect of the fixed outer edge to . Equation (4.1) is solved numerically using classic Newmark method. The fluid module and the solid structure module are computed alternately in a loose coupling manner.
4.3. Results and Discussions
The numerical simulation was carried out using tetrahedral unstructured dynamic grids; a view of the surface grids is shown in Figure 4. A zero thickness 1.6 mm wide virtual baffle face was generated between the outer surface of the main washer and the outlet; see Figure 5(a). At the beginning stage of the movement of the main washer, when the actual displacement is less than 1.6 mm, the virtual baffle face does not participate in grid deforming, and the numerical algorithm described in Section 2.3 is taken into operation. When the actual displacement becomes larger than 1.6 mm, the virtual baffle face would be treated as fully convective ordinary internal cell face and participate in grid deforming as well. Two other virtual baffle faces were preset as well in command valve to achieve its adjusting function; see Figure 5(b).
(a)
(b)
The whole computational geometry includes a 300-litre tank (container) which is not shown in Figure 4. The unsteady computation started at a point of uniform 0.305 Mpa pressure everywhere in the tank and valve, with the tank being inflated at a flux rate of 0.4 kg/s. The computation was terminated after the main washer opened automatically to a balance position. The deformation of the grids on symmetric plane is displayed in Figure 6, from which the motion of main diaphragm can be clearly seen.
(a)
(b)
The calculated result of main washer displacement versus time is shown in Figure 7 and Figure 8 is the displacement of command valve moving assembly versus time. In both Figures 7 and 8, the experimental data is appended for comparison. We can say, from an engineering point of view, that the results of the simulation agree well with the experimental data. The pressure variation history in room A and room B is displayed in Figures 9 and 10, respectively. It is shown in the above figures that the pressure in room A drops down after the main washer opens, and at the moment of opening, a pressure oscillation occurs due to a sudden expansion phenomenon which is attributed to the function of virtual baffle. Comparing Figures 8, 9, and 10, we can see that the connectivity of room A and room B is shut off by the command valve stopper when the displacement of the stopper reaches 0.5 mm, and from then on, the pressure in room A keeps rising while the pressure in room B remains nearly constant. After a short period of waiting at 0.5 mm position, with the pressure in room A rising to a high enough value, the command valve stopper keeps on moving and opens the connectivity of room B with outlet. As a result, the pressure in room B decreases rapidly and finally results in the opening of the main washer. Figure 11 gives a set of sequential images showing the contours of Mach number in the process of main washer opening, from which water-tight sealing effect at the beginning and deflating flux gradually increasing are shown clearly.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
5. Conclusion
A new virtual baffle methodology is implemented for the case where contact/detach problem is involved in fluid and structure interaction simulations. Mathematical expressions in discretization level for this method are established based on tetrahedral unstructured dynamic grids and preconditioned Euler equations. The effectiveness of this new technique is validated by its application for simulation of a complicated safety valve opening process, involving multidisciplinary fluid and structure interaction and dynamic grids. Detailed results of flow field properties and solid structure movement are obtained. The displacements of the moving components in the valve agree well with experimental data quantitatively, and the real working process of the valve is reproduced numerically qualitatively. From the application case, it is proved that the virtual baffle method is successful, it can simulate real detaching process starting from contact, while it has no restriction to minimum grid cell size.