Transcript Diapositiva 1 - Politecnico di Milano

Innovative Finite Element Methods For
Aeroelastic Analysis
Compiled by:
Roberto Flores*
Presented by:
Gabriel Bugeda*
(*) CIMNE: International Center for Numerical Methods in Engineering, Barcelona, Spain
MATEO ANTASME Meeting, 21/05/2007
Main Objectives of Task:
• Analysis of thin walled structures with little or no
bending stiffness subject to unsteady aerodynamic
loads
•Development of efficient FE techniques for the nonlinear (large strain & large displacement) analysis of
membrane behavior, including wrinkling effects
•Improvements to FE flow solvers to allow for fast
solution of complex flow patterns
•Robust coupling of structural solver and CFD codes
for aeroelastic analysis
Structural FE Solver
•Non-linear large displacement/deformation capability
•Features advanced membrane elements including wrinkling
•Implicit dynamic solver (allows for large time steps)
•Total Lagrangian formulation
Inflated airbag showing wrinkles
Sail deployment
CFD Solvers (I)
Implicit incompressible solver for low speed flows
•ALE formulation: Allows for mesh deformation
•Orthogonal subgrid subscale stabilization: Technique developed
at CIMNE. Achieves stabilization with minimum numerical
diffusion by using assumed forms for unresolved flow scales
•Choice of:
•Second-Order Accurate Fractional Step (pressure
segregated) solver
•Monolithic solver
CFD Solvers (II)
Explicit compressible solver for high speed flows
•Edge-based data structure for minimum memory footprint
and optimum performance
•Second order space accuracy
•Explicit multistage Runge-Kutta time integration scheme
•Convective stabilization through limited upwinding
•Implicit residual smoothing for convergence acceleration
•Parallel execution on shared memory architectures via
OPEN-MP directives
Edge oriented data structure
NS equations in conservative form
Φ Fk G k


0
t xk xk
for k  1...3
Ui
0





 u U  p 
 

U 
i
1
i
1
1
i
1




 
Φ  U 2  Fi  uiU 2  p i 2  G i     2i 




 
u
U

p



U
i3 
2i
 i 3


 3


qi   ik uk 
 e 
ui h
Approximate solution using FE discretization
 ~ j
~ 
~j
Φ( x )  N j ( x ) Φ( x )  N j Φ
Weak semi-discrete form of the NS equations
~
~
 ~ j F

G
k
k

N
N
Φ


i
j
 
xk xk


 d  0 for i  1...nnode


The same finite element interpolation is used for fluxes:
~
~ j
~
Fk  N j Fk ( x )  N j Fkj
Solving for the nodal unknowns yields:
~
Φ  M 1 r
M   N i N j d
~
~
~
~ ~ ~
r i  d kij Fkij  bkij Fkij  cki Fki (Fkij  Fki  Fkj )

d kij 
1
1
ij
i


N
N

N
N
d

b


N
N
n
d

c
   N i N i nk d
i
,
k
j
i
j
,
k
k
i
j
k
k
2 
2 

The coefficients dij and bij are non-zero only for pairs nodes connected by an edge ( i.e.
nodes belonging to the same element).
The resulting algorithm is equivalent to a finite volume scheme in which the interface
flux is the average of the nodal values of the edge. Furthermore, for any interior node:
d kij  d kji
bkij  cki  0
thus, the scheme is conservative because the total contribution of internal edges to the
residual is zero.
The basic scheme is equivalent to a centered finite difference stencil which is inherently
unstable due to the odd-even decoupling phenomenon.
The interface fluxes are modified according to Roe’s upwind scheme in order to
suppress instabilities:

~i ~ j 1
~ j ~ i
F  Fk  Fk  A u ij Φ  Φ
2
ij
k



 
~ i ~ i si
Φ Φ 
1  ksi i  1  ksi ij
4

~ j ~i
  Φ Φ
ij
~ i  ij
  Φ  l
i

The factor k controls the extrapolation order for the interface fluxes, which can range
from first to third order.
The coefficients si represent the flux limiters which revert the scheme to first order
near discontinuities and sharp gradients. In areas where the flow field is smooth the
high order scheme is used instead. The limiters are calculated from the ratio of the
solution gradients at the ends of the edge.
Coupled Euler+Boundary Layer Solver
•Solution of viscous problems at high Re numbers requires use of turbulence
models and hybrid meshes to resolve the boundary layer
•Preparation of a suitable mesh is a lengthy task which cannot be easily automated
•To reduce computational costs and speed up the preprocessing stage a coupled
Euler+Boundary Layer solved has been developed
•Uses boundary mesh of 3D volume to create a “virtual” hybrid boundary layer
mesh (extruded prisms)
•In order to capture 3D effects no integral solution is sought, 3D boundary layer
equations are solved directly
•Mapping of arbitrary 3D surface to a plane using unstructured surface mesh
considered too involved  Flux balances calculated in global coordinate system
and projected to local curvilinear coordinates at each point.
•Cell-centered finite volume scheme
•Boundary layer solution coupled to external inviscid flow through transpiration
boundary conditions
Finite Volume Discretization
“Virtual” boundary layer cell
Ci
n ij 
C
j
ni
ni  n j
ni  n j
h
lij
Cij  h n ij  l ij
n
j
nij
Outer boundary of Euler 3D mesh
The flow of a conservative variable from cell i to cell j is then calculated as:
~
 F  n d  Fij  Cij
ij
Solution scheme for boundary layer equations
Solve approximate momentum equation in global coordinate system
dU i
~
Vi
  Fij  Cij  U*i t  t 
dt
Remove normal component

U*i *  U*i  ni U*i  ni

n
Correct momentum using continuity equation
U i  U*i *  n i U in
U  
n
i
di
0
U 


U

dn
This integral is calculated
establishing the mass balance for
the cell


Coupling of boundary layer solution with external flow
Determine displacement thickness * and evaluate transpiration velocity
 *
vn  U Euler
s
Remarks:
•As the boundary layer thickness is replaced with a transpiration velocity, the
Euler mesh does not need to be replaced
•The scheme is not self-starting, for cells around a stagnation point a similarity
solution for the flow near a stagnation area is used
•The FV scheme is cell centered whereas the FE algorithm is vertex centered, the
variables can by transferred by means of:
in
V


V
j
j
c
j
Coupled Fluid-Structural Solver
•CIMNE’s Kratos multiphysics development framework enables coupling
of CFD solver with a FEA structural code to analyze dynamic fluidstructure coupling phenomena
•KRATOS has been completely developed in C++ using a modular objectoriented data structure to enable efficient coupling of single field solvers in
a straightforward way
•Features a Python-Based programmable input
•Available coupling strategies:
•STRONG COUPLING: “SAFE” but often computationally expensive,
requires iterative solving strategy
•LOOSE COUPLING: Often considered “UNSAFE”, computational efficiency
is potentially very HIGH
Coupled Fluid-Structure Interaction Problem
Structural Deformation
Change in fluid Boundary conditions
Change in the pressure field
•Boundary conditions for the fluid are not known until the structure
displacement is calculated
BUT
•Loads on the structure cannot be determined until the flow field has been
solved for
Coupled “Fractional Step” Strategy
•Structural Prediction
Prediction is done by SOLVING the structure subject to a
predicted pressure field (the simplest choice is the pressure at
the end of the previous step)
•Mesh movement step
•Fluid Solution
•Structural Correction
It follows the same rationale as the fractional step (pressure
segregation) procedures used for the solution of the NavierStokes equations
Error due to the coupling algorithm
Assuming that the pressure can be described in the form
p  t   M ae x  Cae x  Kae x
 
and that the structural time integrator can be expressed in a form of the type
yn1  Ayn +Lnfn +Ln+1f n+1
xn1 
y n1  v 
n1 




it is possible to express the solution of the coupled problem as





y 
yn+1 
 n  +E  y 
 = Aex 
 n

pn+1 


 pn 



where yn is an error term, for the coupling procedure to be stable this term must not
grow without bounds
The amplification factor of the error term is

M ae
M
Convergence is achieved when this factor is less than one
Remark: The amplification factor does not depend on the particular time
integration scheme selected
The basic scheme:
Mxi1  Cxi1  Kxi1  p  xi , xi , xi ,t 

can be replaced with

i1  Cxi1  Kxi1  p  xi , xi , xi ,t   Mxi
M


M
x






the procedure remains consistent, as there is no change when Δt→0
Inserting the assumed form of the pressure into the modified algorithm we have
   i

i

1
i

1
i

1
Mx  Cx  Kx 
Mx  Cae xi  Kae xi
1 
now the scheme is stable when
  1
1

•
By choosing an appropriate value for
irrespective of the mass ratio
•
A suitable value can be estimated from the structure of the stiffness matrix of
the fluid problem
the procedure can be made stable
Remark: Fluid and structural meshes need not be congruent, therefore loads on
the structure are calculated remapping the flow solution. Loads are transferred
by means of:
S
S
S
F
N
i
N j t j d 

N
i
N j (τ v  pn) d

where NS and NF represent the shape functions for the structural and fluid
meshes respectively. This is a conservative mapping scheme in the sense that
energy conservation is preserved.
Fvon Karman = 3.7Hz
Fcoupled simulation = 3.05Hz
Fcoupled experiment = 3.10Hz
Example: Flag Flutter
Example applications
Main topic of interest is structural membranes (e.g. inflatable
structures & airbags)
Deployment of inflatable structure
Airbag deployment
Contact algorithms have been implemented to analyze problems
involving solids impacting the membranes
Solid contacting inflatable structure
Solid impacting airbag
(blue ball is attached to membrane)
Thank you for your attention