Transcript Diapositiva 1 - Politecnico di Milano

Coupled Fluid-Structural Solver
•CFD incompressible flow solver has been coupled with a FEA code
to analyze dynamic fluid-structure coupling phenomena
•CFD solver uses pressure splitting and orthogonal subgrid subscale
stabilization to prevent non-physical behavior
•Structural solver suitable for large deformation analysis, including
wrinkling of structural membranes
•All coding has been done in C++ inside the KRATOS objectoriented development framework for multiphysics analysis
•KRATOS allows for completely modular code design and ease of
incorporation of new capabilities
KRATOS
•An environment for implementing innovative computational methods
•Under continual development at CIMNE specifically to address
coupled problems
•Based on Object Oriented Approach using C++
•Features a Python-Based programmable input
•Featured coupling strategies:
•STRONG COUPLING “SAFE” but often computationally expensive,
requires iterative solving strategy
•LOOSE COUPLING  is often considered “UNSAFE”, computational
efficiency is potentially very HIGH
Incompressible CFD Solver
•ALE formulation
•Orthogonal subgrid subscale stabilization
•Choice of:
•Second-Order Accurate Fractional Step solver
•Monolithic solver
Structural Solver
•Non-linear large displacement/deformation capability
•Features advanced membrane elements including wrinkling
•Total Lagrangian model
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 subjected to a
predicted pressure field (the simplest choice is the pressure at
the end of the step before)
•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
x

n

1

yn1  Ayn +Lnfn +Ln+1f n+1 y n1   v 


 n1 
it is possible to express the solution of the coupled problem “in
the future” 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 the amplification 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  the procedure can be
made stable irrespective of the mass ratio
•
A suitable value can be estimated from the structure of the
stiffness matrix of the fluid problem
Example: Flag Flutter
flag flutter
0.01
tip disp
0.005
0
0
5
10
15
-0.005
-0.01
time
Fvon Karman = 3.7Hz
Fcoupled simulation = 3.05Hz
Fcoupled experiment = 3.10Hz
20
25
Example: 2D & 3D driven cavity with deformable base
PUMI 3D CFD Solver
Capabilities Overview:
•Finite element unstructured compressible flow solver
•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
Algorithm Overview
NS equations in conservative form
Φ Fk G k


0
t xk xk
for k  1...3
Ui
0





 u U  p 
 

U 
i1 
1i
 i 1


 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

u2 
T


U i   ui , e    cvT   , h  e  p, qi  k
, ij  2 eij  vkk ij
2
xi

Weak form of the NS equations
  Φ Fk G k
W ( x )  t  xk  xk

 d  0  W

Finite element discretization
 ~ j
~ 
~j
Φ( x )  N j ( x ) Φ( x )  N j Φ


W ( x)  Ni ( x)
Weak semi-discrete form
~
~
 ~ j F
G k
k

N
N
Φ


 i  j
xk xk

 d  0 for i  1...nnode


The numerical fluxes are now approximated by
~
~ j
~j
Fk  N j Fk ( x )  N j Fk
 ~ j N j ~ j ~ j
 Ni  N jΦ  xk Fk  G k



 d  0 for i  1...nnode

by introducing the mass matrix, the last expression can be solved for the
time derivatives of the nodal variables
~ j
Φ  M 1 r
M   N i N j d

r   Ni

N j
xk

~j ~ j
d Fk  G k

To improve computational efficiency the residual is split two parts
N j
~j
~i
r    Ni N j ,k d Fk   Ni Ni ,k d Fk where N j ,k 
xk
j i 

i
(from this point on, no sum is assumed on i, and all sums on j are carried out for j≠i)
integrating by parts and rearranging
r   N i ,k
i

~ ij
~i
~j
N j d Fk   N i ,k N j d Fk   N i N j nk d Fk 


1
~i
~ ij ~ i ~ j
  N i N i nk d Fk where Fk  Fk  Fk
2
the expression must now be symmetrized to realize the benefits of the edge
data structure
1
~ ij
~ ij
r   N i ,k N j  N i N j ,k d Fk   N i N j nk d Fk 
2

i
  N i ,k

~i
~j 1
~i
N j d Fk   N i N j nk d Fk   N i N i nk d Fk
2

using the shape function property Ni  1   N j and after some manipulation
j i
~ ij
ij ~ ij
i ~i
r  d Fk  bk Fk  ck Fk
i
ij
k
1
d   N i ,k N j  N i N j ,k d
2
ij
k
1
b    N i N j nk d
2
ij
k
cki    N i N i nk d

Please remark that
d kij  d kji
bkij  cki  0 foranyinternalnodei
thus, only one coefficient need be stored for each internal edge (pair of connected
nodes, i.e. nodes belonging to the same element)
The scheme is conservative because for any given edge e connecting two internal
nodes i-j, the total contribution the residual is zero
~
ji ~ ji
ij ~ ij
ij ~ ij
re  re  d F  d k Fk  d k Fk  d k Fk  0
i
j
ij ij
k k
When solving a viscous problem, nodal values the solution gradient are required to
obtain the nodal diffusive fluxes. These can be recovered by means of a smoothing
step. Using the regular FE interpolation for the gradients we set
N N
i

d  k Φ   N i N j ,k d Φ
j
j

 k Φ j  M 1  N i N j ,k d Φ j

j
It is well known that the basic Galerkin discretization is inherently unstable (it is
equivalent to a centered difference scheme). To overcome this limitation the
interface fluxes are modified according to Roe’s upwind scheme
~ ij
~i ~ j 1
ij
ij  ij
Fk  Fk  Fk  Fk  A u ij  u k
2
 ij
 ij  j  i
 ij
l
~ j ~i
ij
u   ij
l  x x   Φ Φ
l
where the matrix A u represents the positive flux jacobian along the direction of
the edge, evaluated at the Roe average state between states i and j
ij
This scheme, while stable, provides only first order space accuracy. The amount of
artificial dissipation must be reduced. Two additional states i+ and j- are
introduced


~i ~ j 1
~ j  ~ i   ij
F  Fk  Fk  A u ij Φ  Φ u k
2
ij
k
from the backward and forward extrapolated differences
~ i  ij
  Φ  l
i

~ j  ij
  Φ  l
j

the new interface states are calculated as




~ i ~ i 1
Φ  Φ  1  k  i  1  k  ij
4
~ j ~ j 1
Φ  Φ  1  k   j  1  k  ij
4
where the parameter k controls the degree of approximation.
Near discontinuities the scheme must revert to first order. This is accomplished by
limiting the degree of extrapolation



 
~ i ~ i si
Φ Φ 
1  ksi i  1  ksi ij
4
0  loworder schem e
i
s 
1  highorder schem e
There are many possible choices for the limiting parameter. As an example, we
show here the van Albada limiter
i
ij

2


 


i

s  max0 ,

i 2
ij 2

      

   
  1
Remark: It is in theory possible to achieve a higher accuracy by calculating the
interface fluxes using the extrapolated values, i.e.
 
 


~ i
~ j 1
~ j  ~ i   ij
F  Fk Φ  Fk Φ  A u ij Φ  Φ u k
2
ij
k
however there is usually little difference in practice, so this enhancement can de
omitted without noticeable loss of accuracy
Time integration is performed using a n-stage Runge-Kutta scheme
~
Ψ0  Φ
ti
Ψ j  Ψ 0   j ti 1  ti  M 1r (Ψ j 1 )
~
Φ  Ψn
ti 1
n  1
this scheme is conditionally stable, the nodal allowable time step is calculated as
2

hi
hi 


ti  minCFL 
, 
u a  



a 
2
p

hi being the nodal size, n the maximum fluid diffusivity and CFL the allowable
Courant number
By means of implicit residual smoothing the allowable time step can be increased

r i  ri    r j  r i

onlyfor j connectedtoi
j
which is solved using Jacobi iterations
rn 
i
r i    rn 1
j
j
1   1
j
Time derivatives (and solution gradients) can be solved for very efficiently using
the following iterative process
Md   ij   Ni N j d
j 
~
M Φ0  r
~
~

d ~
M Φ m  Φ m1  r  MΦ m1
d


Example: Transonic flow over a commercial airliner test model