PARALLEL COMPUTING FOR

FLUID STRUCTURE INTERACTION PROBLEMS

Introduction

The prediction of nonlinear aeroelastic phenomena requires the coupled solution of fluid and solid dynamics problems. While such a phenomena can be predicted by using a strongly coupled approach, where the coupled field equations for fluid and solid may so solved together, such an approach provides great deal of numerical complexity due to vast differences in the scales of spatial and temporal variations of each medium. Strong coupling at the interfaces usually dictates matching grid points, which would result with major computational inefficiencies. Since the fluid medium typically requires very dense meshes on the solid boundaries, the use of such dense meshes in solid media would lead to unnecessary inefficiencies. As an alternative, the use of a loosely coupled approach has been popular for solution of such problems, where the fluid and solid media are solved separately and the information between the two media exchanged at the solid-fluid interface by interpolations. While this approach is more approximate, since typically there is a one-step stagger in time between the transient solutions of fluid and solid media, it provides a practical and accurate enough alternative. In that case, the fluid problem is solved via a computational fluid dynamics (CFD) code and the solid problem is solved via a computational solid dynamics (CSD) code. The codes communicate with each other at the solid-fluid interfaces where the nodes may not match, since typically coarser meshes are used for the solid and denser meshes used for the fluid.

The loosely coupled approach requires coupling of the meshes at the interfaces which can be achieved by interpolation of nodal quantities from the nodes of one mesh to the nodes of the other mesh. In spite of the relative efficiency of this approach, one is faced with the coupling of two codes, which are usually developed independently. The problem becomes, even more difficult if a domain decomposition based parallel computing approach is utilized for large scale computations.

Code Coupling

  • Parallel Coupling Interface
    • A parallel coupling interface which has been developed in - house, is used for interaction of CFD and CSD codes
  • CFD code for flow solution
    • Cell centered finite volume formulation
    • A parallel Arbitrary Lagrangian and Eulerian (ALE) form of the flow equations
    • Implicit time integrations
    • A domain - decomposition based parallel algorithm
    • A dynamically deforming mesh
  • A CSD code for solution of structural response with :
    • Finite element formulation
    • Lagrangian formulation of solid mechanics equations
    • Structured or unstructured meshes with mid-surface (shell) elements for flexible
    • Structures, such as wings, or solid elasticity elements for non-slender structures
    • Modal superposition algorithm for dynamic response of structures

Coupling configuration for the codes is shown in Figure 1. Since, the CSD code is usually coarse, the computational burden added is much less than the computational requirements of CFD meshes. Hence, only the CFD mesh is partitioned in this study for parallel computing. For parallelization, the CFD mesh is subdivided into sub domains called solution blocks with one-cell overlaps at the interfaces as illustrated in Figure 2. The message passing across the interfaces of the solution blocks of the CFD mesh is handled by using coupling interface via the MPI message passing library. Since the structural computations are very fast, the structural code for modal superposition is assigned to one of the processors used for the CFD code. For the algorithms developed here, more than one process may be assigned to each processor. Once the interface meshes are defined, coupling interface manages the information exchange, such as pressures and displacements, between each mesh via bilinear interpolations. For solid models with mid-plane representation of the structure using shell-like finite elements, virtual meshes are also needed to facilitate the information transfer. Schematic in Figure 3 depicts the CFD and CSD meshes and their virtual counterparts. A staggered approach is used where the structural analysis follows the flow analysis as follows:

1- Start with an initial condition (typically a steady state flow) at t = 0

2- Compute pressures on the nodes of the CFD mesh from flow calculations

3- Pass the load information to the mesh points of the CSD domain via the virtual structural mesh

4- Calculate the nodal displacements with the CSD code

5- Feed the structural deformation information to the CFD domain via the virtual flow mesh

6- Deform the CFD mesh

7- Advance time: 8. Repeat steps 2 through 8


The flow code uses an implicit backward Euler solver for time integrations, while the CFD code uses a modal superposition algorithm. The first few natural frequencies and mode shapes of the structural mesh are computed only once before the coupled time-integrations start. Transient structural responses are computed via modal superposition of generalized displacements. Since, generally the structural analysis is a small fraction of the flow code, the both codes advance with the same time step governed by the accuracy and stability of the flow mesh. This paper focuses on the coupling approach used in this research.


Figure 1. CFD - CSD coupling scheme



Figure 2. Partitioning of the CFD Domain in to sub domains



Figure 3. CFD surface mesh, CSD virtual surface mesh and structural mesh



Figure 4. Flow Chart for solid fluid interaction solution

COMPUTATIONAL FLUID DYNAMICS SOLVER

The Ale Formulation Of The Three-Dimensional Time-Dependent Inviscid Fluid-Flow Equations


Here W represents the physical domain with a boundary  The vectors Q and F are given by

,

Here nx , ny , nz , are the Cartesian components of the exterior surface unit outward normal vector n on the boundary V = u i + v j + w k is the fluid velocity, W = xt i + yt j + zt k is the mesh velocity and W = W . n = xt nx + yt ny + zt nz is the face speed in the normal direction. Pressure can be expressed as


Solution Method

The spatial integration employed in the flow solver is the cell-centered finite volume formulation. The volume-averaged values are adopted to represent the flow variables. For time integration Implicit Backward-Euler Time-Differencing is used.

Geometric Conservation

Due to mesh movements, a geometric conservation law has to be solved in addition to the mass, momentum, and energy conservation laws. This law is expressed in integral form as :


This geometric conservation law provides a self-consistent solution for the local cell volumes and is solved using the same scheme used for all other flow conservation equations.

Boundary Conditions

Characteristic boundary conditions are applied on the farfield using Riemann invariants. For moving boundaries, the wall boundary conditions are modified by taking the mesh movements into account. Thus, the flow tangency condition is imposed by calculating the flow velocity on solid walls from:

The wall pressure is calculated from the normal momentum equation as:


where, aw is the acceleration of the body, subscripts c and w denote cell-centered and wall values, respectively. The above boundary conditions reduce to steady flow conditions by setting the mesh velocity W to zero.

COMPUTATIONAL STRUCTURAL DYNAMIC SOLVER

The finite element equations for dynamic response of a structural element can be expressed as :


[M],[C]and[K] are called the element mass, damping and stiffness matrices, respectively; {R} is the element load vector which can be calculated by solving the fluid-flow equations in aeroelasticity problems; {q} = {q(t)}is the nodal displacement vector which is a function of time t in transient problems; In static analysis, since

{ } = { } = 0, Equation can be rewritten as:


The time history solution can be obtained by using direct integration or mode superposition. Direct integration is the most general method for dynamic problems. Since in each time step full coupled system of equations are needed to be solved for equilibrium, the computational time is very high even with just a few hundreds degrees of freedoms. An alternate but more effective approach is the mode superposition method. In this method, a set of orthogonal vectors are evaluated first, and the large global equilibrium equations are simplified to a small set of uncoupled second order differential equations. The computational cost per time step is much lower when solving these equations than direct integration. For this reason, the mode superposition method is used to reduce the computational time. The first step in the mode superposition analysis is to obtain the generalized eigenvalue solution:


where is free vibration frequency and {Ø} is mode shape vector, respectively. One assumption in superposition analysis is that the n lowest few vibration modes are enough to describe the dynamic response. Hence, the displacement vector {q(t)} can be expressed as :


where [A] is the matrix whose columns are eigenvectors, and {X(t)} is the generalized displacements, n is the number of free-vibration modes. Considering all above discussions, the equation for the ith mode can be written as:


where Xi are generalized displacements,is the damping ratio corresponding to the ith mode and Fi* = {Ø}iT {F, is the generalized aerodynamic load vector. By this method, the solution of Equation is simplified to solve n linear uncoupled second order differential equations.

Your browser does not support inline frames or is currently configured not to display inline frames.


Figure 5. Dynamic Aeroelastic Analysis of AGARD wing 445.6 (Structural
Deformation + Mach Number Distribution)