oomph-lib

In this example we consider the flow in a 2D channel past a cylinder with an attached elastic "flag". This is the FSI benchmark problem proposed by Turek & Hron,

The problem presented here was used as one of the test cases for oomph-lib's FSI preconditioner; see

In this tutorial we concentrate on the problem formulation. The application of the preconditioner is discussed elsewhere -- the required source code is contained in the driver code.

The Problem

We non-dimensionalise all length and coordinates on the diameter of the cylinder, $ d^* $

, the velocities on the mean velocity, $ U^* $

, and the fluid pressure on the viscous scale. To faciliate comparisons with Turek & Hron's dimensional benchmark data (particularly for the period of the self-excited oscillations), we use a timescale of $T^* = 1 \ \mbox{sec}$ to non-dimensionalise time. The fluid flow is then governed by the non-dimensional Navier-Stokes equations

$Re \left( St \frac{\partial u_i}{\partial t} + u_j \frac{\partial u_i}{\partial x_j} \right) = - \frac{\partial p}{\partial x_i} + \frac{\partial }{\partial x_j} \left[ \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) \right],$

We describe the deformation of the elastic flag by the non-dimensional position vector ${\bf R}(\xi^1, \xi^2, t)$ which is determined by the principle of virtual displacements

$\int_{v} \left\{ \sigma^{ij} \ \delta \gamma_{ij} - \left( {\bf f} - \Lambda^2 \frac{\partial^2 {\bf R}}{\partial t^{2}} \right) \cdot \delta {\bf R} \right\} \ dv -\oint_{A_{tract}} {\bf t} \cdot \delta {\bf R} \ dA =0, \ \ \ \ \ \ \ \ (2)$

where all solid stresses and tractions have been non-dimensionalised on Young's modulus, $ E^* $

; see the Solid Mechanics Tutorial for details. The solid mechanics timescale ratio (the ratio of the timescale $ T^* $

chosen to non-dimensionalise time, to the intrinsic timescale of the solid) can be expressed in terms of the Reynolds and Strouhal numbers, the density ratio, and the FSI interaction parameter as

$\Lambda^2 = \left( \frac{d^*}{T^*} \sqrt{\frac{\rho_s}{E^*}} \right)^2 = St^2 \left( \frac{\rho_s}{\rho_f} \right) Re \ Q.$

Parameter values for the benchmark problems

Geometry

Non-dimensional parameters

Results

FSI2

The constantly adapted mesh contains and average of 65,000 degrees of freedom. A relatively large timestep of $\Delta t = 0.01$ -- corresponding to about 50 timesteps per period of the oscillation -- was used in this computation. With this discretisation the system settles into oscillations with a period of $\approx 0.52$ and an amplitude of the tip-displacement of $0.01 \pm 0.83.$

FSI3

This computation was performed with a timestep of $\Delta t = 0.005$ and resulted in oscillations with a period of $\approx 0.19$ and an amplitude of the tip-displacement of $0.01 \pm 0.36.$

The increase in frequency and Reynolds number leads to the development of thinner boundary and shear layers which require a finer spatial resolution, involving an average of 84,000 degrees of freedom.

Overview of the driver code

Parameter values for the benchmark problems

//=====start_of_global_parameters=================================
/// Global variables
//================================================================
namespace Global_Parameters
{

 /// Default case ID
 string Case_ID="FSI1";

 /// Reynolds number (default assignment for FSI1 test case)
 double Re=20.0; 

 /// Strouhal number (default assignment for FSI1 test case)
 double St=0.5;

 /// \short Product of Reynolds and Strouhal numbers (default
 /// assignment for FSI1 test case)
 double ReSt=10.0;

 /// FSI parameter (default assignment for FSI1 test case)
 double Q=1.429e-6;
 
 /// \short Density ratio (solid to fluid; default assignment for FSI1
 /// test case)
 double Density_ratio=1.0; 

 /// Height of flag
 double H=0.2;
 
 /// x position of centre of cylinder
 double Centre_x=2.0;

 /// y position of centre of cylinder
 double Centre_y=2.0;

 /// Radius of cylinder
 double Radius=0.5;

 /// Pointer to constitutive law
 ConstitutiveLaw* Constitutive_law_pt=0;

 /// \short Timescale ratio for solid (dependent parameter
 /// assigned in set_parameters())
 double Lambda_sq=0.0;

 /// Timestep
 double Dt=0.1;

 /// Ignore fluid (default assignment for FSI1 test case)
 bool Ignore_fluid_loading=false;

 /// Elastic modulus
 double E=1.0;

 /// Poisson's ratio
 double Nu=0.4;

We also include a gravitational body force for the solid. (This is only used for the solid mechanics test cases, CSM1 and CSM2, which will not be discussed here.)


 /// Non-dim gravity (default assignment for FSI1 test case)
 double Gravity=0.0; 

 /// Non-dimensional gravity as body force
 void gravity(const double& time, 
              const Vector<double> &xi,
              Vector<double> &b)
 {
  b[0]=0.0;
  b[1]=-Gravity;
 }

The domain geometry and flow field are fairly complex and it is difficult to construct a good initial guess for the Newton iteration. To ensure its convergence at the beginning of the simulation we therefore employ the method suggested by Turek & Hron: We start the flow from rest and ramp up the inflow profile from zero to its maximum value. The parameters for the time-dependent increase in the influx are defined here:

 
 /// Period for ramping up in flux
 double Ramp_period=2.0;
 
 /// Min. flux 
 double Min_flux=0.0; 
 
 /// Max. flux
 double Max_flux=1.0;

 /// \short Flux increases between Min_flux and Max_flux over 
 /// period Ramp_period
 double flux(const double& t)
 {    
  if (t<Ramp_period)
   {
    return Min_flux+(Max_flux-Min_flux)*
     0.5*(1.0-cos(MathematicalConstants::Pi*t/Ramp_period));
   }
  else
   {
    return Max_flux;
   }
 } // end of specification of ramped influx

Finally, we provide a helper function that assigns the parameters for the various test cases, depending on their ID ("FSI1", "FSI2", "FSI3", "CSM1" or "CSM2"). Here is the assigment for the case FSI1:

 

 /// Set parameters for the various test cases
 void set_parameters(const string& case_id)
 {

  // Remember which case we're dealing with
  Case_ID=case_id;

  // Setup independent parameters depending on test case
  if (case_id=="FSI1")
   {
    // Reynolds number based on diameter of cylinder
    Re=20.0;

    // Strouhal number based on timescale of one second
    St=0.5;

    // Womersley number
    ReSt=Re*St;

    // FSI parameter
    Q=1.429e-6;
    
    // Timestep -- aiming for about 40 steps per period
    Dt=0.1;

    // Density ratio
    Density_ratio=1.0;

    // Gravity
    Gravity=0.0;
    
    // Max. flux
    Max_flux=1.0;

    // Ignore fluid
    Ignore_fluid_loading=false;
 
    // Compute dependent parameters
    
    // Timescale ratio for solid 
    Lambda_sq=Re*Q*Density_ratio*St*St;
  }

In the interest of brevity we omit the listings of the assignments for the other cases. Finally, we select the length of the time interval over which the influx is ramped up from zero to its maximum value to be equal to 20 timesteps, create a constitutive equation for the solid, and document the parameter values used in the simulation:

  // Ramp period (20 timesteps)
  Ramp_period=Dt*20.0;

  // "Big G" Linear constitutive equations:
  Constitutive_law_pt = new GeneralisedHookean(&Nu,&E);
  
  // Doc
  oomph_info << std::endl;
  oomph_info << "-------------------------------------------" 
             << std::endl;
  oomph_info << "Case: " << case_id << std::endl;
  oomph_info << "Re            = " << Re << std::endl;
  oomph_info << "St            = " << St << std::endl;
  oomph_info << "ReSt          = " << ReSt << std::endl;
  oomph_info << "Q             = " << Q << std::endl;
  oomph_info << "Dt            = " << Dt << std::endl;
  oomph_info << "Ramp_period   = " << Ramp_period << std::endl;
  oomph_info << "Max_flux      = " << Max_flux << std::endl;
  oomph_info << "Density_ratio = " << Density_ratio << std::endl;
  oomph_info << "Lambda_sq     = " << Lambda_sq << std::endl;
  oomph_info << "Gravity       = " << Gravity << std::endl;
  oomph_info << "Ignore fluid  = " << Ignore_fluid_loading<< std::endl;
  oomph_info << "-------------------------------------------"
             << std::endl << std::endl;
 }

}// end_of_namespace

The driver code

//=======start_of_main==================================================
/// Driver 
//======================================================================
int main(int argc, char* argv[])
{
 // Store command line arguments
 CommandLineArgs::setup(argc,argv);

 // Get case id as string
 string case_id="FSI1";
 if (CommandLineArgs::Argc==1)
  {
   oomph_info << "No command line arguments; running self-test FSI1" 
              << std::endl;
  }
 else if (CommandLineArgs::Argc==2)
  {
   case_id=CommandLineArgs::Argv[1];
  }
 else
  {
   oomph_info << "Wrong number of command line arguments" << std::endl;
   oomph_info << "Enter none (for default) or one (namely the case id" 
              << std::endl;
   oomph_info << "which should be one of: FSI1, FSI2, FSI3, CSM1" 
              << std::endl;
  } 
 std::cout << "Running case " << case_id << std::endl;

We set up the global parameter values, create a DocInfo object and trace file to record the output, and build the problem.

 // Setup parameters for case identified by command line
 // argument
 Global_Parameters::set_parameters(case_id);

 // Prepare output
 DocInfo doc_info;
 ofstream trace_file; 
 doc_info.set_directory("RESLT");
 trace_file.open("RESLT/trace.dat");
 
 // Length and height of domain
 double length=25.0;
 double height=4.1;

 //Set up the problem
 TurekProblem<AlgebraicElement<RefineableQTaylorHoodElement<2> >,
  RefineableQPVDElement<2,3> > problem(length, height);

Next, we choose the number of timesteps (using a smaller number for a validation run, and for the case FSI1 in which the system rapidly approaches a steady state) and initialise the time-stepping for an impulsive start from the zero flow solution.

 // Default number of timesteps
 unsigned nstep=4000;
 if (Global_Parameters::Case_ID=="FSI1")
  {
   std::cout << "Reducing number of steps for FSI1 " << std::endl;
   nstep=400;
  }

 if (CommandLineArgs::Argc==1)
  {
   std::cout << "Reducing number of steps for validation " << std::endl;
   nstep=2;
  }

 //Timestep: 
 double dt=Global_Parameters::Dt;

 // Initialise timestep 
 problem.initialise_dt(dt);

 // Impulsive start
 problem.assign_initial_values_impulsive(dt);

Finally, we document the initial condition and start the time-stepping procedure, setting the first flag to false because we have not specified an analytical expression for the initial conditions that could be re-assigned after the mesh adapatation when computing the first timestep.

 // Doc the initial condition
 problem.doc_solution(doc_info,trace_file);
 doc_info.number()++; 
 
 // Don't re-set the initial conditions when adapting during first
 // timestep
 bool first = false;
 
 // Max number of adaptation for time-stepping
 unsigned max_adapt=1;
 
 for(unsigned i=0;i<nstep;i++)
  { 
  // Solve the problem
  problem.unsteady_newton_solve(dt,max_adapt,first); 
  
  // Output the solution
  problem.doc_solution(doc_info,trace_file);
  
  // Step number
  doc_info.number()++;
 }
 
 trace_file.close(); 

}//end of main

The Problem class

//====start_of_problem_class=========================================== 
/// Problem class
//====================================================================== 
template< class FLUID_ELEMENT,class SOLID_ELEMENT >
class TurekProblem : public Problem
{

public:

 /// \short Constructor: Pass length and height of domain
 TurekProblem(const double &length, const double &height);
 
 /// Access function for the fluid mesh 
 RefineableAlgebraicCylinderWithFlagMesh<FLUID_ELEMENT>* fluid_mesh_pt() 
  { return Fluid_mesh_pt;}

 /// Access function for the solid mesh
 ElasticRefineableRectangularQuadMesh<SOLID_ELEMENT>*& solid_mesh_pt() 
  {return Solid_mesh_pt;} 

 /// Access function for the i-th mesh of FSI traction elements
 SolidMesh*& traction_mesh_pt(const unsigned& i)
  {return Traction_mesh_pt[i];} 
 
 /// Actions after adapt: Re-setup the fsi lookup scheme
 void actions_after_adapt();

 /// Doc the solution
 void doc_solution(DocInfo& doc_info, ofstream& trace_file);

 /// Update function (empty)
 void actions_after_newton_solve() {}

 /// Update function (empty)
 void actions_before_newton_solve(){}

 /// \short Update the (enslaved) fluid node positions following the
 /// update of the solid variables before performing Newton convergence
 /// check
 void actions_before_newton_convergence_check();

 /// Update the time-dependent influx
 void actions_before_implicit_timestep();

private:

 /// Create FSI traction elements
 void create_fsi_traction_elements();

 /// Pointer to solid mesh
 ElasticRefineableRectangularQuadMesh<SOLID_ELEMENT>* Solid_mesh_pt;

 ///Pointer to fluid mesh
 RefineableAlgebraicCylinderWithFlagMesh<FLUID_ELEMENT>* Fluid_mesh_pt;

 /// Vector of pointers to mesh of FSI traction elements
 Vector<SolidMesh*> Traction_mesh_pt;

 /// Combined mesh of traction elements -- only used for documentation
 SolidMesh* Combined_traction_mesh_pt;

 /// Overall height of domain
 double Domain_height;

 /// Overall length of domain
 double Domain_length;

 /// Pointer to solid control node
 Node* Solid_control_node_pt;

 /// Pointer to fluid control node
 Node* Fluid_control_node_pt;
 
};// end_of_problem_class

The problem constructor

//=====start_of_constructor============================================= 
/// Constructor: Pass length and height of domain
//====================================================================== 
template< class FLUID_ELEMENT,class SOLID_ELEMENT >
TurekProblem<FLUID_ELEMENT,SOLID_ELEMENT>::
TurekProblem(const double &length,
             const double &height) :  Domain_height(height),
                                      Domain_length(length)
 
{
 // Increase max. number of iterations in Newton solver to
 // accomodate possible poor initial guesses
 Max_newton_iterations=20;
 Max_residuals=1.0e4;

 // Build solid mesh
 //------------------

 // # of elements in x-direction
 unsigned n_x=20;

 // # of elements in y-direction
 unsigned n_y=2;

 // Domain length in y-direction 
 double l_y=Global_Parameters::H;

 // Create the flag timestepper (consistent with BDF<2> for fluid)
 Newmark<2>* flag_time_stepper_pt=new Newmark<2>;
 add_time_stepper_pt(flag_time_stepper_pt); 

 /// Left point on centreline of flag so that the top and bottom
 /// vertices merge with the cylinder.
 Vector<double> origin(2);
 origin[0]=Global_Parameters::Centre_x+
  Global_Parameters::Radius*
  sqrt(1.0-Global_Parameters::H*Global_Parameters::H/
       (4.0*Global_Parameters::Radius*Global_Parameters::Radius));
 origin[1]=Global_Parameters::Centre_y-0.5*l_y;

 // Set length of flag so that endpoint actually stretches all the
 // way to x=6:
 double l_x=6.0-origin[0];

 //Now create the mesh
 solid_mesh_pt() = new ElasticRefineableRectangularQuadMesh<SOLID_ELEMENT>(
  n_x,n_y,l_x,l_y,origin,flag_time_stepper_pt);

We create an error estimator for the solid mesh and identify a control node at the tip of the flag to track its motion.

 
 // Set error estimator for the solid mesh
 Z2ErrorEstimator* solid_error_estimator_pt=new Z2ErrorEstimator;
 solid_mesh_pt()->spatial_error_estimator_pt()=solid_error_estimator_pt;


 // Element that contains the control point
 FiniteElement* el_pt=solid_mesh_pt()->finite_element_pt(n_x*n_y/2-1);

 // How many nodes does it have?
 unsigned nnod=el_pt->nnode();

 // Get the control node
 Solid_control_node_pt=el_pt->node_pt(nnod-1);

 std::cout << "Coordinates of solid control point " 
           << Solid_control_node_pt->x(0) << " " 
           << Solid_control_node_pt->x(1) << " " << std::endl;

Finally, we perform one uniform mesh refinement and disable any further mesh adaptation.

 // Refine the mesh uniformly
 solid_mesh_pt()->refine_uniformly();

 //Do not allow the solid mesh to be refined again
 solid_mesh_pt()->adapt_flag()=false;

Next, we attach FSISolidTractionElements to the boundaries of the solid mesh that are exposed to the fluid. We complete their build by specifying which boundary of the bulk mesh they are attached to, as this information is required when setting up the fluid-structure interaction; see Further comments and exercises.


 // Build mesh of solid traction elements that apply the fluid
 //------------------------------------------------------------
 // traction to the solid elements
 //-------------------------------

 // Create storage for Meshes of FSI traction elements at the bottom
 // top and left boundaries of the flag
 Traction_mesh_pt.resize(3);
 
 // Now construct the traction element meshes
 Traction_mesh_pt[0]=new SolidMesh;
 Traction_mesh_pt[1]=new SolidMesh;
 Traction_mesh_pt[2]=new SolidMesh;

 // Build the FSI traction elements
 create_fsi_traction_elements();

 // Loop over traction elements, pass the FSI parameter and tell them 
 // the boundary number in the bulk solid mesh -- this is required so 
 // they can get access to the boundary coordinates!
 for (unsigned bound=0;bound<3;bound++)
  {
   unsigned n_face_element = Traction_mesh_pt[bound]->nelement();
   for(unsigned e=0;e<n_face_element;e++)
    {
     //Cast the element pointer and specify boundary number
     FSISolidTractionElement<SOLID_ELEMENT,2>* elem_pt=
     dynamic_cast<FSISolidTractionElement<SOLID_ELEMENT,2>*>
      (Traction_mesh_pt[bound]->element_pt(e));

     // Specify boundary number
     elem_pt->set_boundary_number_in_bulk_mesh(bound);

     // Function that specifies the load ratios
     elem_pt->q_pt() = &Global_Parameters::Q;
    
    }
  } // build of FSISolidTractionElements is complete

Finally, we create GeomObject representations of the three surface meshes of FSISolidTractionElements. We will use these to represent the curvilinear, moving boundaries of the fluid mesh.


 // Turn the three meshes of FSI traction elements into compound
 // geometric objects (one Lagrangian, two Eulerian coordinates)
 // that determine the boundary of the fluid mesh
 MeshAsGeomObject<1,2,FSISolidTractionElement<SOLID_ELEMENT,2> >*
  bottom_flag_pt=
  new MeshAsGeomObject<1,2,FSISolidTractionElement<SOLID_ELEMENT,2> >
  (Traction_mesh_pt[0]);
 
 MeshAsGeomObject<1,2,FSISolidTractionElement<SOLID_ELEMENT,2> >* tip_flag_pt=
  new MeshAsGeomObject<1,2,FSISolidTractionElement<SOLID_ELEMENT,2> >
  (Traction_mesh_pt[1]);
 
 MeshAsGeomObject<1,2,FSISolidTractionElement<SOLID_ELEMENT,2> >* top_flag_pt=
  new MeshAsGeomObject<1,2,FSISolidTractionElement<SOLID_ELEMENT,2> >
  (Traction_mesh_pt[2]);

The final mesh to be built is the fluid mesh whose constructor requires pointers to the four GeomObjects that represent the cylinder and three fluid-loaded faces of the flag, respectively. We represent the cylinder by a Circle object:


 // Build fluid mesh
 //-----------------

 //Create a new Circle object as the central cylinder
 Circle* cylinder_pt = new Circle(Global_Parameters::Centre_x,
                                  Global_Parameters::Centre_y,
                                  Global_Parameters::Radius);

We build the mesh and identify a control node (a node at the upstream face of the cylinder), before creating an error estimator and performing one uniform mesh refinement.

 // Allocate the fluid timestepper
 BDF<2>* fluid_time_stepper_pt=new BDF<2>;
 add_time_stepper_pt(fluid_time_stepper_pt);
 
 // Build fluid mesh
 Fluid_mesh_pt=
  new RefineableAlgebraicCylinderWithFlagMesh<FLUID_ELEMENT>
  (cylinder_pt, 
   top_flag_pt,
   bottom_flag_pt,
   tip_flag_pt,
   length, height, 
   l_x,Global_Parameters::H,
   Global_Parameters::Centre_x,
   Global_Parameters::Centre_y,
   Global_Parameters::Radius,
   fluid_time_stepper_pt);
 

 // I happen to have found out by inspection that
 // node 5 in the hand-coded fluid mesh is at the 
 // upstream tip of the cylinder
 Fluid_control_node_pt=Fluid_mesh_pt->node_pt(5);

 // Set error estimator for the fluid mesh
 Z2ErrorEstimator* fluid_error_estimator_pt=new Z2ErrorEstimator;
 fluid_mesh_pt()->spatial_error_estimator_pt()=fluid_error_estimator_pt;

 // Refine uniformly
 Fluid_mesh_pt->refine_uniformly();

We now add the various meshes to the Problem's collection of sub-meshes and combine them to a global mesh


 // Build combined global mesh
 //---------------------------

 // Add Solid mesh the problem's collection of submeshes
 add_sub_mesh(solid_mesh_pt());

 // Add traction sub-meshes
 for (unsigned i=0;i<3;i++)
  {
   add_sub_mesh(traction_mesh_pt(i));
  }

 // Add fluid mesh
 add_sub_mesh(fluid_mesh_pt());
 
 // Build combined "global" mesh
 build_global_mesh();

The application of boundary conditions for the solid are straightforward: All displacements of the flag's left end (mesh boundary 3) are suppressed; the other faces are free. Strictly speaking, the pinning of the redundant solid pressure nodes is superfluous since the RefineableQPVDElement used for the discretisation of the flag employ a displacement-based formulation, but it is good practise to perform this step anyway to "future-proof" the code for the use of other element types.

 


 // Apply solid boundary conditons
 //-------------------------------
 
 //Solid mesh: Pin the left boundary (boundary 3) in both directions
 unsigned n_side = mesh_pt()->nboundary_node(3);
 
 //Loop over the nodes
 for(unsigned i=0;i<n_side;i++)
  {
   solid_mesh_pt()->boundary_node_pt(3,i)->pin_position(0);
   solid_mesh_pt()->boundary_node_pt(3,i)->pin_position(1);
  }
 
 // Pin the redundant solid pressures (if any)
 PVDEquationsBase<2>::pin_redundant_nodal_solid_pressures(
  solid_mesh_pt()->element_pt());

The fluid has Dirichlet boundary conditions (prescribed velocity) everywhere apart from the outflow where only the horizontal velocity is unknown.

 
 // Apply fluid boundary conditions
 //--------------------------------
 
 //Fluid mesh: Horizontal, traction-free outflow; pinned elsewhere
 unsigned num_bound = fluid_mesh_pt()->nboundary();
 for(unsigned ibound=0;ibound<num_bound;ibound++)
  {
   unsigned num_nod= fluid_mesh_pt()->nboundary_node(ibound);
   for (unsigned inod=0;inod<num_nod;inod++)
    {
     // Parallel, axially traction free outflow at downstream end
     if (ibound != 1)
      {
       fluid_mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
       fluid_mesh_pt()->boundary_node_pt(ibound,inod)->pin(1);
      }
     else
      {
       fluid_mesh_pt()->boundary_node_pt(ibound,inod)->pin(1);
      }
    }
  }//end_of_pin
 
 // Pin redundant pressure dofs in fluid mesh
 RefineableNavierStokesEquations<2>::
  pin_redundant_nodal_pressures(fluid_mesh_pt()->element_pt());

We impose a parabolic inflow profile with the current value of the influx at the inlet (fluid mesh boundary 3).


 // Apply boundary conditions for fluid
 //-------------------------------------

 // Impose parabolic flow along boundary 3
 // Current flow rate
 double t=0.0;
 double ampl=Global_Parameters::flux(t);
 unsigned ibound=3; 
 unsigned num_nod= Fluid_mesh_pt->nboundary_node(ibound);
 for (unsigned inod=0;inod<num_nod;inod++)
  {
   double ycoord = Fluid_mesh_pt->boundary_node_pt(ibound,inod)->x(1); 
   double uy = ampl*6.0*ycoord/Domain_height*(1.0-ycoord/Domain_height);
   Fluid_mesh_pt->boundary_node_pt(ibound,inod)->set_value(0,uy);
   Fluid_mesh_pt->boundary_node_pt(ibound,inod)->set_value(1,0.0);    
  }

We complete the build of the solid elements by passing them the pointer to the constitutive equation, the gravity vector and the timescale ratio:

 

 // Complete build of solid elements
 //---------------------------------

 //Pass problem parameters to solid elements
 unsigned  n_element =solid_mesh_pt()->nelement();
 for(unsigned i=0;i<n_element;i++)
  {
   //Cast to a solid element
   SOLID_ELEMENT *el_pt = dynamic_cast<SOLID_ELEMENT*>(
    solid_mesh_pt()->element_pt(i));
   
   // Set the constitutive law
   el_pt->constitutive_law_pt() =
    Global_Parameters::Constitutive_law_pt;
   
   //Set the body force
   el_pt->body_force_fct_pt() = Global_Parameters::gravity;

   // Timescale ratio for solid
   el_pt->lambda_sq_pt() = &Global_Parameters::Lambda_sq;
  }

The fluid elements require pointers to the Reynolds and Womersley (product of Reynolds and Strouhal) numbers:

 
 

 // Complete build of fluid elements
 //---------------------------------

 // Set physical parameters in the fluid mesh
 unsigned nelem=fluid_mesh_pt()->nelement();
 for (unsigned e=0;e<nelem;e++)
  {
   // Upcast from GeneralisedElement to the present element
   FLUID_ELEMENT* el_pt = dynamic_cast<FLUID_ELEMENT*>
    (fluid_mesh_pt()->element_pt(e));
   
   //Set the Reynolds number
   el_pt->re_pt() = &Global_Parameters::Re;   
   
   //Set the Womersley number
   el_pt->re_st_pt() = &Global_Parameters::ReSt;
   
  }//end_of_loop

Setting up the fluid-structure interaction is done from "both" sides" of the fluid-solid interface: First we ensure that the no-slip condition is automatically applied to all fluid nodes that are located on the three faces of the flag (mesh boundaries 5, 6 and 7). This is done by passing the function pointer to the FSI_functions::apply_no_slip_on_moving_wall() function to the relevant fluid nodes (recall that the auxiliary node update functions are automatically executed whenever the position of a node is updated by the algebraic node update). Since the no-slip condition (1) involves the Strouhal number (which, in the current problem, is not equal to the default value of FSI_functions::Strouhal_for_no_slip=1.0), we overwrite the default assignment with the actual Strouhal number in the problem.



 // Setup FSI
 //----------
 
 // Pass Strouhal number to the helper function that automatically applies
 // the no-slip condition
 FSI_functions::Strouhal_for_no_slip=Global_Parameters::St;

 // The velocity of the fluid nodes on the wall (fluid mesh boundary 5,6,7)
 // is set by the wall motion -- hence the no-slip condition must be
 // re-applied whenever a node update is performed for these nodes. 
 // Such tasks may be performed automatically by the auxiliary node update 
 // function specified by a function pointer:
 if (!Global_Parameters::Ignore_fluid_loading)
  {
   for(unsigned ibound=5;ibound<8;ibound++ )
    { 
     unsigned num_nod= Fluid_mesh_pt->nboundary_node(ibound);
     for (unsigned inod=0;inod<num_nod;inod++)
      {   
       Fluid_mesh_pt->boundary_node_pt(ibound, inod)->
        set_auxiliary_node_update_fct_pt(
         FSI_functions::apply_no_slip_on_moving_wall);
      }
    } // done automatic application of no-slip

Next, we set up the lookup schemes required by the FSISolidTractionElements to establish which fluid elements affect the traction onto the solid:

  
   // Work out which fluid dofs affect the residuals of the wall elements:
   // We pass the boundary between the fluid and solid meshes and 
   // pointers to the meshes. The interaction boundary are boundaries 5,6,7
   // of the 2D fluid mesh.
   FSI_functions::setup_fluid_load_info_for_solid_elements<FLUID_ELEMENT,2>
    (this,5,Fluid_mesh_pt,Traction_mesh_pt[0]);  
   
   FSI_functions::setup_fluid_load_info_for_solid_elements<FLUID_ELEMENT,2>
    (this,6,Fluid_mesh_pt,Traction_mesh_pt[2]);  
   
   FSI_functions::setup_fluid_load_info_for_solid_elements<FLUID_ELEMENT,2>
    (this,7,Fluid_mesh_pt,Traction_mesh_pt[1]); 
  }

All interactions have now been specified and we conclude by assigning the equation numbers

 // Assign equation numbers
 cout << assign_eqn_numbers() << std::endl; 


}//end_of_constructor

Create traction elements

//============start_of_create_traction_elements==============================
/// Create FSI traction elements 
//=======================================================================
template<class FLUID_ELEMENT,class SOLID_ELEMENT >
void TurekProblem<FLUID_ELEMENT,SOLID_ELEMENT>::create_fsi_traction_elements()
{

 // Container to collect all nodes in the traction meshes
 std::set<SolidNode*> all_nodes;

 // Traction elements are located on boundaries 0-2:
 for (unsigned b=0;b<3;b++)
  {
   // How many bulk elements are adjacent to boundary b?
   unsigned n_element = solid_mesh_pt()->nboundary_element(b);
   
   // Loop over the bulk elements adjacent to boundary b?
   for(unsigned e=0;e<n_element;e++)
    {
     // Get pointer to the bulk element that is adjacent to boundary b
     SOLID_ELEMENT* bulk_elem_pt = dynamic_cast<SOLID_ELEMENT*>(
      solid_mesh_pt()->boundary_element_pt(b,e));
     
     //What is the index of the face of the element e along boundary b
     int face_index = solid_mesh_pt()->face_index_at_boundary(b,e);
          
     // Create new element and add to mesh
     Traction_mesh_pt[b]->add_element_pt(
      new FSISolidTractionElement<SOLID_ELEMENT,2>(bulk_elem_pt,face_index));
     
    } //end of loop over bulk elements adjacent to boundary b

   // Identify unique nodes
   unsigned nnod=solid_mesh_pt()->nboundary_node(b);
   for (unsigned j=0;j<nnod;j++)
    {
     all_nodes.insert(solid_mesh_pt()->boundary_node_pt(b,j));
    }
  }

 // Build combined mesh of fsi traction elements
 Combined_traction_mesh_pt=new SolidMesh(Traction_mesh_pt);
 
 // Stick nodes into combined traction mesh
 for (std::set<SolidNode*>::iterator it=all_nodes.begin();
      it!=all_nodes.end();it++)
  {
   Combined_traction_mesh_pt->add_node_pt(*it);
  }

} // end of create_traction_elements

Actions before Newton convergence check

//====start_of_actions_before_newton_convergence_check===================
/// Update the (enslaved) fluid node positions following the
/// update of the solid variables
//=======================================================================
template <class FLUID_ELEMENT,class SOLID_ELEMENT> 
void TurekProblem<FLUID_ELEMENT,SOLID_ELEMENT>
::actions_before_newton_convergence_check ()
{
 fluid_mesh_pt()->node_update();
}

Actions before the timestep

//===== start_of_actions_before_implicit_timestep=========================
/// Actions before implicit timestep: Update inflow profile
//========================================================================
template <class FLUID_ELEMENT,class SOLID_ELEMENT> 
void TurekProblem<FLUID_ELEMENT,SOLID_ELEMENT>::
actions_before_implicit_timestep()
{
 // Current time
 double t=time_pt()->time();
 
 // Amplitude of flow
 double ampl=Global_Parameters::flux(t);
 
 // Update parabolic flow along boundary 3
 unsigned ibound=3; 
 unsigned num_nod= Fluid_mesh_pt->nboundary_node(ibound);
 for (unsigned inod=0;inod<num_nod;inod++)
  {
   double ycoord = Fluid_mesh_pt->boundary_node_pt(ibound,inod)->x(1); 
   double uy = ampl*6.0*ycoord/Domain_height*(1.0-ycoord/Domain_height);
   Fluid_mesh_pt->boundary_node_pt(ibound,inod)->set_value(0,uy);
   Fluid_mesh_pt->boundary_node_pt(ibound,inod)->set_value(1,0.0);    
  }

} //end_of_actions_before_implicit_timestep

Actions after adapt

//=====================start_of_actions_after_adapt=======================
///  Actions after adapt: Re-setup FSI
//========================================================================
template<class FLUID_ELEMENT,class SOLID_ELEMENT >
void TurekProblem<FLUID_ELEMENT,SOLID_ELEMENT>::actions_after_adapt()
{
 // Unpin all pressure dofs
 RefineableNavierStokesEquations<2>::
  unpin_all_pressure_dofs(fluid_mesh_pt()->element_pt());
 
 // Pin redundant pressure dofs
 RefineableNavierStokesEquations<2>::
  pin_redundant_nodal_pressures(fluid_mesh_pt()->element_pt());

 // Unpin all solid pressure dofs
 PVDEquationsBase<2>::
  unpin_all_solid_pressure_dofs(solid_mesh_pt()->element_pt());
 
 // Pin the redundant solid pressures (if any)
 PVDEquationsBase<2>::pin_redundant_nodal_solid_pressures(
  solid_mesh_pt()->element_pt());


 // The velocity of the fluid nodes on the wall (fluid mesh boundary 5,6,7)
 // is set by the wall motion -- hence the no-slip condition must be
 // re-applied whenever a node update is performed for these nodes. 
 // Such tasks may be performed automatically by the auxiliary node update 
 // function specified by a function pointer:
 if (!Global_Parameters::Ignore_fluid_loading)
  {
   for(unsigned ibound=5;ibound<8;ibound++ )
    { 
     unsigned num_nod= Fluid_mesh_pt->nboundary_node(ibound);
     for (unsigned inod=0;inod<num_nod;inod++)
      {   
       Fluid_mesh_pt->boundary_node_pt(ibound, inod)->
        set_auxiliary_node_update_fct_pt(
         FSI_functions::apply_no_slip_on_moving_wall);
      }
    }
 

   // Re-setup the fluid load information for fsi solid traction elements
   FSI_functions::setup_fluid_load_info_for_solid_elements<FLUID_ELEMENT,2>
    (this,5,Fluid_mesh_pt,Traction_mesh_pt[0]); 
   
   FSI_functions::setup_fluid_load_info_for_solid_elements<FLUID_ELEMENT,2>
    (this,6,Fluid_mesh_pt,Traction_mesh_pt[2]); 
   
   FSI_functions::setup_fluid_load_info_for_solid_elements<FLUID_ELEMENT,2>
    (this,7,Fluid_mesh_pt,Traction_mesh_pt[1]); 
  }

 
}// end of actions_after_adapt

Post-processing

//=====start_of_doc_solution========================================
/// Doc the solution
//==================================================================
template<class FLUID_ELEMENT,class SOLID_ELEMENT >
void TurekProblem<FLUID_ELEMENT,SOLID_ELEMENT>::doc_solution(
 DocInfo& doc_info, ofstream& trace_file)
{
 
//  FSI_functions::doc_fsi<AlgebraicNode>(Fluid_mesh_pt,
//                                        Combined_traction_mesh_pt,
//                                        doc_info);

//  pause("done");

 ofstream some_file;
 char filename[100];

 // Number of plot points
 unsigned n_plot = 5; 

 // Output solid solution
 sprintf(filename,"%s/solid_soln%i.dat",doc_info.directory().c_str(),
         doc_info.number());
 some_file.open(filename);
 solid_mesh_pt()->output(some_file,n_plot);
 some_file.close();
 
 // Output fluid solution
 sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
         doc_info.number());
 some_file.open(filename);
 fluid_mesh_pt()->output(some_file,n_plot);
 some_file.close();


//Output the traction
 sprintf(filename,"%s/traction%i.dat",doc_info.directory().c_str(),
         doc_info.number());
 some_file.open(filename);
// Loop over the traction meshes
 for(unsigned i=0;i<3;i++)
  {
   // Loop over the element in traction_mesh_pt
   unsigned n_element = Traction_mesh_pt[i]->nelement();
   for(unsigned e=0;e<n_element;e++)
    {
     FSISolidTractionElement<SOLID_ELEMENT,2>* el_pt = 
      dynamic_cast<FSISolidTractionElement<SOLID_ELEMENT,2>* > (
       Traction_mesh_pt[i]->element_pt(e) );
     
     el_pt->output(some_file,5);
    }
  }
 some_file.close();  


 // Write trace (we're only using Taylor Hood elements so we know that
 // the pressure is the third value at the fluid control node...
 trace_file << time_pt()->time() << " " 
            << Solid_control_node_pt->x(0) << " " 
            << Solid_control_node_pt->x(1) << " " 
            << Fluid_control_node_pt->value(2) << " "
            << Global_Parameters::flux(time_pt()->time()) << " " 
            << std::endl;
 
 cout << "Doced solution for step " 
      << doc_info.number() 
      << std::endl << std::endl << std::endl;

}//end_of_doc_solution

	$Re = U^* d^* \rho_f/\mu$	$Q = \mu U^* / (E^* d^*)$	$\rho_s/\rho_f$	$\Lambda^2 = (d^* /T^* \sqrt{\rho_s/E^*})^2 = St^2 (\rho_s/\rho_f) Re \ Q$
FSI1		$1.429 \times 10^{-6}$		$7.145 \times 10^{-6}$
FSI2		$7.143 \times 10^{-6}$		$7.143 \times 10^{-6}$
FSI3		$3.571 \times 10^{-6}$		$1.786 \times 10^{-6}$

Demo problem: Turek & Hron's FSI benchmark problem

The Problem

Parameter values for the benchmark problems

Geometry

Non-dimensional parameters

Results

FSI2

FSI3

Overview of the driver code

Parameter values for the benchmark problems

The driver code

The Problem class

The problem constructor

Create traction elements

Actions before Newton convergence check

Actions before the timestep

Actions after adapt

Post-processing

Further comments and exercises

Acknowledgements

Source files for this tutorial

PDF file