Introduction

This page provides code blocks of several examples that may be used to perform fluid structure interaction computations.

Internal structures and corresponding initial conditions

Initialization

The lstelt array is allocated.

 
!===============================================================================
! 1.  Initialization
!===============================================================================
 
! Allocate a temporary array for boundary faces selection
allocate(lstelt(nfabor))
 

Example

The following code block provides an example of definition of two internal structures.

Here one fills array idfstr(nfabor). For each boundary face ifac, idfstr(ifac) is the number of the structure the face belongs to (if idfstr(ifac) = 0, the face ifac doesn't belong to any structure. When using internal coupling, structure number necessarily needs to be positive (as shown in following examples).

The number of "internal" structures is automatically defined with the maximum value of idfstr table, meaning that internal structure numbers must be defined sequentially with positive values, beginning with integer value 1.

In the following example, boundary faces with color 4 belong to internal structure 1. Boundary faces with color 2 belong to internal structure 2. The total number of internal structures equals 2. The boundary faces are identified using the getfbr subroutine.

 
!===============================================================================
! 2.  Definition of internal structures
!===============================================================================
 
call getfbr('4', nlelt, lstelt)
 
do ilelt = 1, nlelt
 
  ifac = lstelt(ilelt)
 
  idfstr(ifac) = 1
 
enddo
 
call getfbr('6', nlelt, lstelt)
 
do ilelt = 1, nlelt
 
  ifac = lstelt(ilelt)
 
  idfstr(ifac) = 2
 
enddo

For each internal structure one can here define :

an initial velocity vstr0
an initial displacement xstr0 (i.e. xstr0 is the value of the displacement xstr compared to the initial mesh at time t=0)
a displacement compared to equilibrium \ref xstreq. \ref xstreq is the initial displacement of the internal structure compared to its position at equilibrium; at each time step t and for a displacement xstr(t), the associated internal structure will be subjected to a force due to the spring:
$-k (xstr(t)+xstreq).$

Note that xstr0, \ref xstreq and vstr0 arrays are initialized at the beginning of the calculations to the value of 0.

When starting a calculation using ALE, or re-starting a calculation with ALE basing on a first calculation without ALE, an initial iteration 0 is automatically calculated in order to take initial arrays xstr0, vstr0 and \ref xstreq into account. In another case, set the following option

italin=1

in subroutine usipsu, so that the code can deal with arrays xstr0, vstr0 or \ref xstreq.

In the following example :

internal structure 1 has got an initial displacement xstr0 of 2 meters in the y direction and a displacement compared to equilibrium \ref xstreq of 1 meter in the y direction.
Initial velocity vstr0 in the z direction of structure 2 equals -0.5 m/s.

 
if (ntpabs.le.1) then  ! Avoid resetting in case of restart
  xstr0(2,1)  = 2.d0
  xstreq(2,1) = 1.d0
  vstr0(3,2)  =-0.5d0
endif
 

Here one can modify the values of the prediction coefficients for displacements anf fluid forces in internal FSI coupled algorithm.

The use of these coefficients is reminded here :

predicted displacement
$\vect{X_{n+1}} = \vect{X_n} + aexxst \Delta t \vect{X^{\prime}_n} + bexxst \Delta t (\vect{X^{\prime}_n}-\vect{X^{\prime}_{n-1}})$
with $\vect{X_n}$ standing for the displacement at iteration , $\vect{X^{\prime}_n}$ and $\vect{X^{\prime}_{n-1}}$ representing the internal structure velocity respectively at iteration and .
fluid force sent to structure internal solver
$\vect{F_{n+1}} = cfopre \vect{F_n} + (1-cfopre) \vect{F_{n-1}}$
$\vect{F_n}$ and $\vect{F_{n-1}}$ standing for the fluid force acting on the structure respectively at iteration and .

 
aexxst =  0.5d0
bexxst =  0.0d0
cfopre =  2.d0
 

Activation of structural history output (i.e. displacement, structural velocity, structural acceleration and fluid force) (ihistr=0, disabled ; ihistr=1, enabled) The value of structural history output step is the same as the one for standard variables (nthist).

 
ihistr = 1
 

Structural parameters in case of Fluid Structure internal coupling

Note that the subroutine usstr2 is called at each time step of the calculation.

For each internal structure one defines here :

its mass $\tens{M}$ (xmstru)
its friction coefficient $\tens{C}$ (xcstru)
its stiffness $\tens{K}$ (xkstru)

forstr array gives fluid stresses acting on each internal structure. Moreover it's possible to take external forces (gravity for example) into account, too.

\ref xstr array indicates the displacement of the structure compared to its position in the initial mesh.

xstr0 array gives the displacement of the structures in initial mesh compared to structural equilibrium.

\ref vstr array stands for structural velocity.

\ref xstr, xstr0, and \ref vstr arrays are data tables that can be used to define arrays mass, friction and stiffness. THOSE ARE NOT TO BE MODIFIED.

The 3D structural equation that is solved is the following one :

$\tens{M} \vect{X^{\prime\prime}} + \tens{C} \vect{X^{\prime}} + \tens{K} (\vect{X}+\vect{X_0}) = \vect{F}$

where $\vect{X}$ stands for the structural displacement compared to initial mesh position (\ref xstr) and $\vect{X_0}$ represents the displacement of the structure in initial mesh compared to equilibrium.

Note that $\tens{M}$ , $\tens{C}$ and $\tens{K}$ are 3x3 matrices.

This equation is solved using a Newmark HHT algorithm. Note that the time step used to solve this equation (dtstr) can be different from the one of fluid calculations. user is free to define dtstr array. At the beginning of the calculation dtstr is initialized to the value of \ref dt (fluid time step).

Initialization

The matrices xmstru, xcstru and xkstru are set to zero.

 
!===============================================================================
! Definition of structural parameters: mass, friction, stiffness and stresses.
!===============================================================================
 
! --- Matrices xmstru, xcstru and xkstru are initialized to the value of 0.
do istr = 1, nbstru
 
  do ii = 1, 3
    do jj = 1, 3
      xmstru(ii,jj,istr) = 0.d0
      xcstru(ii,jj,istr) = 0.d0
      xkstru(ii,jj,istr) = 0.d0
    enddo
  enddo
 
enddo
 

Example 1

Two examples of definition of the mass, the friction coefficient and the stiffness of an internal structure are provided hereafter.

 
! --- Example 1): In following example structure '1' is defined as an isotropic
!     system (i.e. matrices M, C and K are diagonal) : mass equals 5 kg,
!     stiffness equals 2 N/m and friction
!     =  =     =
!     coefficient equals 3 kg.s .
 
do ii = 1, 3
  xmstru(ii,ii,1) = 5.d0
  xcstru(ii,ii,1) = 2.d0
  xkstru(ii,ii,1) = 3.d0
enddo
 

Example 2

 
! --- Example 2): In this example structure '2' is subjected to the following
!                 movement :
!               - In plane xOy the movement is locally defined along an axis
!                 (OX). Structural parameters in X direction are called
!                 xm, xc and xk. The angle of inclination between global (Ox)
!                 axis and local (OX) axis is called theta. Movement in local (OY)
!                 direction is imposed to be rigid.
!               - In 'z' direction the movement is modeled to be oscillating and
!                 harmonic (meaning that there is no friction). Mass equals 1. kg
!                 and stiffness equals 1. N/m. Fluid stresses in that direction
!                 are taken into account. Moreover the structure is also
!                 subjected to an external oscillating
!                 force Fz_ext = 3 * cos(4*t).
 
 
!                 This leads to the following local equations :
!                 xm.X'' + xc.X' + xk.X = FX
!                                     Y = 0
!                    Z''         +    Z = FZ + 3.cos(4.t)
 
theta = pi/6.d0
cost = cos(theta)
sint = sin(theta)
 
! FX, FY, and FZ stand for the local fluid forces components.
! They are defined as follows, using gobal components of
! fluid forces Fx, Fy and Fz.
!   fx =  cost*Fx + sint*Fy
!   fy = -sint*Fx + cost*Fy
!   fz = Fz
 
! After changing of basis, the problem can be described as follows,
! using global coordinates:
 
xm = 1.d0
xc = 3.d-1
xk = 2.d0
fx = forstr(1,2)
fy = forstr(2,2)
 
xmstru(1,1,2) = xm*cost**2
xmstru(1,2,2) = xm*cost*sint
xmstru(2,1,2) = xm*cost*sint
xmstru(2,2,2) = xm*sint**2
xmstru(3,3,2) = 1.d0
 
xcstru(1,1,2) = xc*cost**2
xcstru(1,2,2) = xc*cost*sint
xcstru(2,1,2) = xc*cost*sint
xcstru(2,2,2) = xc*sint**2
 
xkstru(1,1,2) = (xk-1.d0)*cost**2 + 1.d0
xkstru(1,2,2) = (xk-1.d0)*cost*sint
xkstru(2,1,2) = (xk-1.d0)*cost*sint
xkstru(2,2,2) = (xk-1.d0)*sint**2 + 1.d0
xkstru(3,3,2) = 1.d0
 
forstr(1,2) = fx*cost**2   + fy*sint*cost
forstr(2,2) = fx*sint*cost + fy*sint**2
forstr(3,2) = forstr(3,2) + 3.d0*cos(4.d0*ttcabs)
 
do istr = 1, nbstru
  dtstr(istr) = dtcel(1)
enddo
 

Fluid Structure external coupling with Code_Aster

Initialization

The lstelt array is allocated.

 
!===============================================================================
! 1.  initialization
!===============================================================================
 
! Allocate a temporary array for boundary faces selection
allocate(lstelt(nfabor))
 

Example

In the following example, two sets of boundary faces that will be coupled with Code_Aster are defined as well as the fluid force components which are given to structural calculations.

Here one fills array idfstr(nfabor) For each boundary face ifac, idfstr(ifac) is the number of the structure the face belongs to (if idfstr(ifac) = 0, the face ifac doesn't belong to any structure.) When using external coupling with Code_Aster, structure number necessarily needs to be negative (as shown in following examples).

The number of "external" structures with Code_Aster is automatically defined with the maximum absolute value of idfstr table, meaning that external structure numbers must be defined sequentially with negative values beginning with integer value -1.

In following example, boundary faces with color 2 and which abscissa $ x < 2 $ belong to external structure -1. Boundary faces with color 2 and which abscissa $ x > 2 $ belong to external structure -2. The total number of external structures coupled with Code_Aster equals 2.

 
!===============================================================================
! 2.  Definition of external structures
!===============================================================================
 
call getfbr('2 and X < 2.0', nlelt, lstelt)
 
do ilelt = 1, nlelt
 
  ifac = lstelt(ilelt)
 
  idfstr(ifac) = -1
 
enddo
 
call getfbr('2 and X > 2.0', nlelt, lstelt)
 
do ilelt = 1, nlelt
 
  ifac = lstelt(ilelt)
 
  idfstr(ifac) = -2
 
enddo
 
! The movement of external structure called '-1' is blocked in Z direction.
 
asddlf(3,1) = 0
 
! The movement of external structure called '-2' is blocked in Z direction.
 
asddlf(3,2) = 0