Introduction

This page provides code blocks of several examples that may be used to perform Fluid-Structure Interaction (FSI) computations.

Declaring mobile structures.

In addition to implicit declarations with boundary conditions in the GUI, internal structures may be declared in the cs_user_model or cs_user_parameters function. For example, to declare 3 coupled structures:

cs_mobile_structures_add_n_structures(3);

cs_mobile_structures_add_n_structures

void cs_mobile_structures_add_n_structures(int n_structures)

Add internal mobile structures.

Definition: cs_mobile_structures.c:1046

Coupling with code_aster can be declared in a similar manner, using:

cs_ast_coupling_add();

cs_ast_coupling_add

void cs_ast_coupling_add(void)

Define coupling with code_aster.

Definition: cs_ast_coupling.c:603

Selecting coupled faces.

Whether coupled structures are defined using the GUI or as above function, their association with boundary zones may be completed or modified in the cs_user_fsi_structure_num function from the cs_user_fluid_structure_interaction.c file.

In the following example, 2 different structures are each associated to a boundary group:

  int n_structs = 2;
  const char *name[] = {"wall_1", "wall_2"};
 
  for (int st_id = 0; st_id < n_structs; st_id++) {
 
    const cs_zone_t  *zn = cs_boundary_zone_by_name(name[st_id]);
 
    for (cs_lnum_t e_idx = 0; e_idx < zn->n_elts; e_idx++) {
 
      const cs_lnum_t face_id = zn->elt_ids[e_idx];
      structure_num[face_id] = (st_id + 1);
 
    }
 
  }

This type of association is not necessary using the GUI in simple cases, but could be used to associate multiple boundary zones to a common structure for example, or conversely, associated different structures to different portions of a given zone (for example, if a single zone contains multiple tubes, so as to simplify fluid BC definitions, separate structures could be associated to each tube using sub-selections).

In the next example, the same face groups are coupled with an external (code_aster) structure.

  int n_structs = 2;
  const char *name[] = {"wall_1", "wall_2"};
 
  for (int st_id = 0; st_id < n_structs; st_id++) {
 
    const cs_zone_t  *zn = cs_boundary_zone_by_name(name[st_id]);
 
    for (cs_lnum_t e_idx = 0; e_idx < zn->n_elts; e_idx++) {
 
      const cs_lnum_t face_id = zn->elt_ids[e_idx];
      structure_num[face_id] = - (st_id + 1);
 
    }
 
  }

Note that in this case, only the (negative) sign of the structure number is used, not the actual number, so both groups of faces are coupled to the same code_aster instance.

Internal structure parameters

Various global and individual parameters associated to internal structures can be controlled with functions from cs_user_fluid_structure_interaction.c.

Initial positions and parameters

Initial positions and other parameters associated with internal structure coupling may be defined using the cs_user_fsi_structure_define function.

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 xstreq. 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).$

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 xstreq into account. In another case, set the following option

cs_glob_ale_need_init = 1;

cs_glob_ale_need_init

int cs_glob_ale_need_init

in the cs_user_parameters function, so that the code can deal with arrays xstr0, vstr0 or xstreq.

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

In the following example:

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

  if (! is_restart) { /* Avoid resetting in case of restart */
    xstr0[0][1]  = 2.0;
    xstreq[0][1] = 1.0;
    vstr0[1][2]  =-0.5;
  }

Here one can modify the values of the prediction coefficients for displacements and 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 .

The following code snippet redefines theses advanced coefficients for predicted displacements and predicted force:

  *aexxst = 0.5;
  *bexxst = 0.0;
  *cfopre = 2.0;

The time plotting behavior of internal structures may also be modified in this same function. The following example shows how the default output format and interval may be modified:

  *plot = 2;  /* .csv format */
 
  cs_time_control_init_by_time_step(plot_time_control,
                                    -1,
                                    -1,
                                    5,       /* nt_interval */
                                    true,    /* at_start */
                                    false);  /* at end */

The output interval is based on a standard cs_time_control_t logic, so time-based or more advanced user-defined behavior may de used here. It is recommended not to use a too large interval, so that sampling does not miss important signal characteristics.

Note also that the movement of the internal fluid structure is computed using a Newmark method, whose parameters can also be modified by calling the cs_mobile_structures_set_newmark_coefficients function.

Structural parameters

Some internal structure parameters, although usually fixed, may be time-varying in complex models, so they may be defined and changed using the cs_user_fsi_structure_values function, which is called at each time step when internal mobile structures are present.

For each internal structure one defines here:

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

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

The xstr array indicates the displacement of each structure compared to its position in the initial mesh.

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

The vstr array contains the structure velocities.

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

The 3D structure equation that is solved is:

$\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 (xstr) and $\vect{X_0}$ represents the displacement of the structure in the 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 dt (fluid time step).

Example 1

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

  /* In the following example structure 0 is defined as an isotropic
   * system (i.e. matrices M, C and K are diagonal): mass equals 5 kg,
   * damping coefficient equals 2 N.s/m and stiffness equals 3 N/m. */
 
  for (int i = 0; i < 3; i++) {
    xmstru[i][i][0] = 5.0;
    xcstru[i][i][0] = 2.0;
    xkstru[i][i][0] = 3.0;
  }

Example 2

  /* In this example structure '1' 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 the '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)
   */
 
  double theta = cs_math_pi / 6.0;
  double cost = cos(theta);
  double 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: */
 
  double xm = 1.0;
  double xc = 3.e-1;
  double xk = 2.0;
  double fx = forstr[1][0];
  double fy = forstr[1][1];
 
  xmstru[1][0][0] = xm*cost*cost;
  xmstru[1][1][0] = xm*cost*sint;
  xmstru[1][0][1] = xm*cost*sint;
  xmstru[1][1][1] = xm*sint*sint;
  xmstru[1][2][2] = 1.0;
 
  xcstru[1][0][0] = xc*cost*cost;
  xcstru[1][1][0] = xc*cost*sint;
  xcstru[1][0][1] = xc*cost*sint;
  xcstru[1][1][1] = xc*sint*sint;
 
  xkstru[1][0][0] = (xk-1.0)*cost*cost + 1.0;
  xkstru[1][1][0] = (xk-1.0)*cost*sint;
  xkstru[1][0][1] = (xk-1.0)*cost*sint;
  xkstru[1][1][1] = (xk-1.0)*sint*sint + 1.0;
  xkstru[1][2][2] = 1.0;
 
  forstr[1][0] = fx*cost*cost + fy*sint*cost;
  forstr[1][1] = fx*sint*cost + fy*sint*sint;
  forstr[1][2] += 3.0*cos(4.0*ts->t_cur);
 
  for (int i = 0; i < n_structs; i++) {
    dtstr[i] = ts->dt[0];
  }

Fluid Structure external coupling with code_aster

For code_aster coupling, it is sufficient to declare a coupling as described in the earlier section and possibly select associated faces.

Currently, multiple structures associated to faces associated with code_aster are handled in a grouped manner,so the structure number has no importance, as long as it is negative. In the future, in case of coupling with multiple instances of code_aster, we could assign different structures to different coupling instances.