Multidisciplinary investigation by fluid-structure-motion integrated simulation

974 IV International Conference on Computational Methods for Coupled Problems in Science and Engineering COUPLED PROBLEMS 2011 M. Papadrakakis, E. Oñate and B. Schrefler (Eds) MULTIDISCIPLINARY INVESTIGATION BY FLUID-STRUCTUREMOTION INTEGRATED SIMULATION S. TAKAHASHI* AND N. ARAI* * Tokyo University of Agriculture and Technology (TUAT) 2-24-16, Naka-cho, Koganei, Tokyo, Japan e-mail: [email protected], www.tuat.ac.jp/~arailab Key words: Fluid-Structure Interaction, Biomimetic System, Cartesian Grid method, Overset Grid Method. Abstract. This study is devoted to realize the flight system based on biomimetic mechanisms by using multidisciplinary numerical simulation that consists of flow simulation, structural simulation and motion analysis. In this paper, the effective motion is investigated to generate sufficient force to make a flight. The flow simulation is carried out by finite difference and finite volume method in overset grid while the structural analysis is conducted by linear elastic simulation with finite element method. The coupled algorithm is constructed from the weak coupling with inner iteration. As a result, it was revealed that the fluid force is increased by the deformation of the object. 1 INTRODUCTION The motion of the wings is essential to generate lift and thrust together for the flight of living beings. From the viewpoint of the physics, the phenomenon consists of physics of fluid, structure and motion so that one-sided study is not sufficient to investigate it precisely. Therefore, multidisciplinary analysis should be performed to reveal the coupled mechanism and characteristics. The flight mechanism of living beings has been investigated by experimental and numerical studies [1,2,3]. Especially, one of insects acting in low Reynolds number flow may be closer to the swimming in sticky fluid rather than flying in the air. Many researches are conducted to develop a micro aerial vehicle (MAV) based on the insects because of the high agility flight. These studies are attractive to realize the absolute flight envelope including complete vertical takeoff and landing (VTOL). By the way, what is occurred as the scale becomes larger and the flight speed becomes faster? Obviously, the flow acting on the wings becomes turbulent as Reynolds number becomes high. In the case of birds, the frequency of the flapping motion becomes lower and the shape of the wing section becomes closer to the airfoil with increase of Reynolds number [4,5,6]. It is thought that the trend is caused from the limitation of the muscle power and the efficiency of the wings for the forward flight. It is correct that only small birds can achieve VTOL due to the flapping motion with high frequency. The objective of this study is to develop a manned aerial vehicle that can realize VTOL characteristics based on biomimetic flight systems. In this paper, coupled simulation is 975 S. Takahashi and N. Arai. 2 discussed to perform parametric study for investigation of the aerodynamic characteristics. The coupled simulation consists of incompressible Navier-Stokes simulation and linear elastic structural analysis. The coupling system is constructed by weak coupling algorithm with inner iteration. In this paper, low Reynolds number flows are adopted to investigate the validity of the present solver as a primary stage of this work. 2 NUMERICAL METHODOLOGY 2.1 Numerical method for flow simulation The flow simulation is carried out by two-dimensional incompressible Navier-Stokes equations with overset grid composed of boundary-fitted coordinate (BFC) grid and multiblock Cartesian grid. The present Cartesian grid method is called as Building-Cube Method (BCM) that is used for efficient spatial discretization based on equally-spaced Cartesian grid [7,8,9,10]. Although one of the most important points in using a Cartesian grid method is surface approximation of the object, the object is represented by BFC grid to estimate the fluid force correctly in the present study. Furthermore, the numerical dissipation can be suppressed owing to the background Cartesian grid. It can be significant to capture the hysteresis of the present flowfield to estimate the fluid force exactly. Arbitrary Lagrangian Eulerian (ALE) method is utilized in the BFC grid [11], while Eulerian method is adopted in the background BCM grid [12]. Overset grid system in the present simulation is shown in Fig. 1. Blue thick lines are boundaries of cubes including same number of meshes. Fig. 1 Overset grid system of BFC grid and BCM grid The governing equation in BFC grid is expressed in Eq. (1). The grid velocity term is appeared in Eq. (1) because of the ALE method. (){} uuuu u u ∆+ ∇ −=∇⋅−+ ∂ ∂ = ⋅∇ ρ µ ρ p tg 0 (1) The governing equation is divided into three stages for fractional-step method [13] and discretized by finite difference method in regular grid arrangement. The convective flux and 976 S. Takahashi and N. Arai. 3 the diffusive flux are approximated by UTOPIA scheme and second-order central difference scheme, respectively. Temporal velocity field of fractional-step method is solved by CrankNicolson scheme. Poisson equation of the pressure is discretized by second-order central difference scheme. Both the temporal velocity field and the pressure field are solved by SOR method. Finally, real velocity field is solved by the correction of the pressure gradient to the temporal velocity field. These procedures are implemented based on Eq. (2). (){} ** uuuu uu* ∆=∇⋅−+ ∆ − ρ µ n g n n t (2) t pn ∆ ⋅∇ = ∇∗+ u ρ 12 ρ 11 ++ ∇ −= ∆ −nn p t * uu The governing equation in BCM grid is expressed in Eq. (3). In BCM grid, the equation is also divided into three stages as same as BFC grid for fractional-step method. The spatial discretization is implemented by finite volume method in collocated grid arrangement. () uuu u u ∆+ ∇ −=⋅∇+ ∂ ∂ = ⋅ ∇ ρ µ ρ p t 0 (3) The convective flux and diffusive flux in solving the temporal velocity field are approximated by fully-conservative second-order finite difference scheme and second-order central difference scheme, respectively. In adaptive refined Cartesian grid method, the numerical fluxes through the boundary between the different-sized grids may cause critical numerical instability. In the present method, the convective flux is calculated by using firstorder upwind scheme based on the face flux estimated from Rhie-Chow momentum interpolation [14] and the average of the volume flux. The procedure is shown in the schematic as Fig. 2. The cell-averaged volume flux in large grid fL is adopted as fluxes in fringe cells of the small grid. At the same time, the flux of a large fringe cell is calculated as (fS1+fS2)/2. The face fluxes for small and large grids are evaluated as (fL+fS1)/2, (fL+fS1)/2 and {(fS1+fS2)/2+fL}/2 to calculate the convective fluxes. In the case of the face between samesized grids, the convective flux is calculated by fully-conservative second-order central difference scheme. In the second stage of fractional-step method, Poisson equation is solved for the pressure that should be taken care of the treatment since the gradient of the pressure is also used to evaluate the real velocity field. The right hand side of the Poisson equation is calculated by the face fluxes based on the arithmetic average of the volume fluxes in usual collocated grid method. In the present AMR grid system, however, the face fluxes should be continuous through the different-sized grid face. The process is as same as the one for the estimation of the convective fluxes. Though the calculation process is so simple, it has much importance to 977 S. Takahashi and N. Arai. 4 perform stable computation with fully-conservative scheme in AMR grid system. The real velocity field is calculated by the correction using the pressure gradient to the temporal velocity field in the third stage of fractional-step method. Though the contravariant flux is also calculated by Rhie-Chow momentum interpolation at the same time, it should be treated properly at the boundary of different-sized grid. Therefore, the face flux of the large grid is estimated by the same manner as the RHS of the Poisson equation. Fig. 2 Schematic of interpolation of fluxes at different-sized grid boundary The data exchange process is necessary in overset grid system. In the present study, the relationship between a donor cell and a child node is constructed with the motion of the BFC grid repeatedly. Here, the donor cell means data supplier to the child node. The process is performed with weighted averaged interpolation based on the area-coordinate rule as shown in Fig. 3 and Eq. (4). Blue lines and nodes belong to BCM grid, while red lines and nodes belong to BFC grid. A variable of child node VC is interpolated from weighed values by using the areas of triangles of S12, S23, S34 and S41 of the donor cell. Child nodes of the BFC grid are the vertices belonging in outermost cells as shown in the left schematic of Fig. 3, while one of the BCM grid are cell center points surrounded by BFC cells except the outermost nodes as shown in right schematic of Fig. 3. Fig. 3 Relationship between BCM donor cell and BFC child node (left) and BFC donor cell and BCM child node (right) N 1 N 2 N 3 N 4 S 12 S 23 S 34 S 41 N 1 N 2 N 3 N 4 S 12 S 23 S 34 S 41 N 2 N 1 N 3 N 4 S 12 S 23 S 34 S 41 N 2 N 1 N 3 N 4 S 12 S 23 S 34 S 41 fL fL fL fS1 fS2 1 1 2 2 1 1 2fL fL fL fS1 fS2 1 1 2 2 1 1 2 978 S. Takahashi and N. Arai. 5 41342312 42312312412413413423 SSSS VSSVSSVSSVSS VC+++ + + + = (4) 2.2 Numerical method for structural analysis The deformation of the object is solved by two-dimensional finite element method as linear elastic deformation with triangle elements. Shape function is first-order and the relationship between stress and strain assumes plane-stress state. The deformation at each timestep is calculated from the stiffness matrix at initial state without any stress by BiCGstab method. In this procedure, the phase error can be occurred since the fluid force is estimated around deformed object. This problem is solved by implementation of inner iteration for the coupling of the flow analysis and structural analysis as shown in Fig. 4. After the object deformation, the spatial grid for the flow simulation is deformed as the surface deformation of the object by using simple algebraic method. 2.3 Numerical method for motion analysis The motion analysis is carried out to investigate the trajectory and the stability of the flight. However, fluid-structure coupled simulation is discussed in this paper to find out appropriate motion to achieve a vertical flight. Fig. 4 Flowchart of the present coupled simulation (left) and algorithm for solving temporal velocity field (right) input parameters and computational mesh time integration move boundary make overset info solve flowfield end maximum inner iteration of coupling ? solve deformation make overset info next timestep No Yes input parameters and computational mesh time integration move boundary make overset info solve flowfield end maximum inner iteration of coupling ? solve deformation make overset info next timestep No Yes solve temporal velocity u* send u* BCM→BFC send u* BFC→BCM u* data exchange in BCM solve u* in BCM solve u* in BFC converge ? solve pressure field No Yes next inner iteration of temporal velocity u* solve temporal velocity u* send u* BCM→BFC send u* BFC→BCM u* data exchange in BCM solve u* in BCM solve u* in BFC converge ? solve pressure field No Yes next inner iteration of temporal velocity u* 979 S. Takahashi and N. Arai. 6 2.4 Algorithm for coupled simulation As mentioned before, the present coupled simulation is based on fluid-structure simulation to explore the optimal motion to make a flight. The flowcharts in Fig. 4 show the schematics of the present coupled algorithm and flow simulation with overset grid method. The coupled simulation is implemented by weak coupling algorithm with inner iteration. The flow simulation with overset grid method also includes inner iteration in each stage of fractionalstep method to exchange data through the overlapping boundary. Furthermore, the data exchange process is basically needed for BCM grid because of the multi-block structure. 3 RESULTS AND DISCUSSION In this paper, three kinds of flows are discussed. Firstly, validation of the present method is demonstrated by flow around a circular cylinder. Secondly, flow simulation around translating circular cylinder is discussed as an example of a flow around moving object. Finally, flow simulation around translating elastic circular cylinder is investigated by the present coupled simulation. The diameter of the circular cylinder is 0.10[m]. Freestream velocity of the validation and the maximum translational velocity are 5.0×10-2[m/s]. In all cases, viscosity of fluid is assumed as 1.82×10-5[Pa⋅s] so that Reynolds number based on freestream velocity or maximum translational velocity is around 330. Fig. 5 Vorticity distribution at Re# 330 close to cylinder (left) and far from cylinder (right) 3.1 Flow simulation around fixed circular cylinder The validation of the present overset flow solver is conducted by using fixed BFC grid 980 S. Takahashi and N. Arai. 7 with a circular cylinder at Reynolds number 330. Total numbers of grid points of BFC grid and BCM grid are 4,530 (=151×30) and 274,432 (=268[cube]×322[cell]), respectively. Uniform flow is assumed at the inflow boundary while Neumann boundary condition is used for the outflow velocity. The pressure at inflow and outflow boundaries is also determined by Neumann boundary condition. At upper and lower boundaries, slip boundary condition is used for the velocity and Dirichlet boundary condition is used for the pressure. Figure 5 shows fully developed vorticity fields close to the cylinder and far from the cylinder. The connectivity of the vorticity through the overlapping boundary between BFC and BCM grid is sufficiently smooth. The variation of the aerodynamic coefficient is also comparable to general value. The validity of the present flow solver was confirmed from these results. Fig. 6 Aerodynamic coefficient around circular cylinder at Re# 330 3.2 Flow simulation around translating circular cylinder Next simulation is carried out to investigate the aerodynamic characteristics around translating circular cylinder in quiet fluid. Only BFC grid is moved backward and forward along with x-axis by force in background BCM grid. The translational motion is expressed in trigonometric function as Eq. (5). x(t)=Acos( ω t) (5) The value x in the equation (5) represents the center of the circular cylinder. The variable A and ω mean the amplitude of the translation and the angular velocity. The amplitude and the angular velocity are defined as 0.10[m] and 0.50[rad/s] in the present study. The maximum 0.6 0.8 1.0 1.2 1.4 1.6 1.8 0 20 40 60 80 100 120 140 160 time [sec] Cd -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 Cl Cd Cl 981 S. Takahashi and N. Arai. 8 translational velocity A ω is 5.0×10-2[m/s] so that Reynolds number based on the maximum velocity is around 330. Fig. 7 Vorticity field around translating circular cylinder Fig. 8 Aerodynamic force around translating circular cylinder Figure 8 shows aerodynamic forces history in that the instantaneous vorticity fields are 5.0sec 7.5sec 10.5sec 18.0sec 67.4sec 137.2sec -5.E-04 -4.E-04 -3.E-04 -2.E-04 -1.E-04 0.E+00 1.E-04 2.E-04 3.E-04 4.E-04 5.E-04 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 time [sec] force [N] Drag Lift Location -5.E-04 -4.E-04 -3.E-04 -2.E-04 -1.E-04 0.E+00 1.E-04 2.E-04 3.E-04 4.E-04 5.E-04 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 time [sec] force [N] Drag Lift Location 982 S. Takahashi and N. Arai. 9 captured as shown in Fig. 7. The green symbols in the graph show the timings of capturing the flowfield. When the cylinder is accelerated to the right direction at 5.0 sec twin vortices are grown in the wake of the cylinder. Then, the motion is reversed at 7.5 sec to the left direction. The maximum drag force is generated just after the timing since the cylinder moves into the wake made by the previous motion. The twin vortices are left in the wake and it flows to the right direction due to the self-induced velocity at 10.5 sec. After that, vortex streets are occurred in both sides at 67.4 sec. After further, the flowfield shows asymmetrical vortex streets in both sides, gradually. Finally, complete asymmetric vortex streets are produced at 137.2 sec. The phenomenon is corresponded to the lift generation in Fig. 8. Although Karman vortex street is occurred in this Reynolds number in uniform flow as well known, the present asymmetric flow induced by the motion is similar phenomenon with Karman vortex street in self-induced flow. 3.3 Coupled simulation around translating elastic circular cylinder Finally, the fluid force is investigated around an elastic circular cylinder by coupled simulation. Flow conditions are same as previous simulation around rigid circular cylinder. The deformation of the cylinder is solved by 3,000 grid points with 6,000 degrees of freedom with triangular elements. The computational mesh is shown in Fig. 9. Inner circle drawn in red line is treated as fixed points by the assumption that hard material core is included in the structural mesh. Though the shape of element is quadrangle in Fig. 9, it is treated as triangular element in the analysis. Young’s modulus and Poisson ratio are defined as 4×10-2[Pa] and 0.30 to investigate the affect of the deformation, apparently. The width of the material is 1.0[m]. The deformation is solved as linear elastic problem from the initial state at all times. Inner iteration is performed twice at each timestep to suppress the phase error from the weak coupling algorithm. Fig. 9 Computational mesh for structural analysis Mesh for flow simulation Mesh for structural analysis