fsi_pseudo_solid_collapsible_channel.cc
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 
31 // Generic oomph-lib includes
32 #include "generic.h"
33 #include "navier_stokes.h"
34 #include "beam.h"
35 #include "solid.h"
36 #include "constitutive.h"
37 
38 // The wall mesh
40 
41 //Include the fluid mesh
43 
44 using namespace std;
45 
46 using namespace oomph;
47 
48 
49 //======================================================================
50 /// Upgrade mesh to solid mesh
51 //======================================================================
52 template <class ELEMENT>
54  public virtual CollapsibleChannelMesh<ELEMENT>,
55  public virtual SolidMesh
56 {
57 
58 
59 public:
60 
61  /// \short Constructor: Build mesh and copy Eulerian coords to Lagrangian
62  /// ones so that the initial configuration is the stress-free one.
64  const unsigned& ncollapsible,
65  const unsigned& ndown,
66  const unsigned& ny,
67  const double& lup,
68  const double& lcollapsible,
69  const double& ldown,
70  const double& ly,
71  GeomObject* wall_pt,
72  TimeStepper* time_stepper_pt=
73  &Mesh::Default_TimeStepper) :
74  CollapsibleChannelMesh<ELEMENT>(nup, ncollapsible, ndown, ny,
75  lup, lcollapsible, ldown, ly,
76  wall_pt,
77  time_stepper_pt)
78  {
79  /// Make the current configuration the undeformed one by
80  /// setting the nodal Lagrangian coordinates to their current
81  /// Eulerian ones
82  set_lagrangian_nodal_coordinates();
83  }
84 
85 };
86 
87 
88 
89 
90 //==========start_of_BL_Squash =========================================
91 /// Namespace to define the mapping [0,1] -> [0,1] that re-distributes
92 /// nodal points across the channel width.
93 //======================================================================
94 namespace BL_Squash
95 {
96 
97  /// Boundary layer width
98  double Delta=0.1;
99 
100  /// Fraction of points in boundary layer
101  double Fract_in_BL=0.5;
102 
103  /// \short Mapping [0,1] -> [0,1] that re-distributes
104  /// nodal points across the channel width
105  double squash_fct(const double& s)
106  {
107  // Default return
108  double y=s;
109  if (s<0.5*Fract_in_BL)
110  {
111  y=Delta*2.0*s/Fract_in_BL;
112  }
113  else if (s>1.0-0.5*Fract_in_BL)
114  {
115  y=2.0*Delta/Fract_in_BL*s+1.0-2.0*Delta/Fract_in_BL;
116  }
117  else
118  {
119  y=(1.0-2.0*Delta)/(1.0-Fract_in_BL)*s+
120  (Delta-0.5*Fract_in_BL)/(1.0-Fract_in_BL);
121  }
122 
123  return y;
124  }
125 }// end of BL_Squash
126 
127 
128 
129 //////////////////////////////////////////////////////////////////////////
130 //////////////////////////////////////////////////////////////////////////
131 //////////////////////////////////////////////////////////////////////////
132 
133 
134 //====start_of_underformed_wall============================================
135 /// Undeformed wall is a steady, straight 1D line in 2D space
136 /// \f[ x = X_0 + \zeta \f]
137 /// \f[ y = H \f]
138 //=========================================================================
139 class UndeformedWall : public GeomObject
140 {
141 
142 public:
143 
144  /// \short Constructor: arguments are the starting point and the height
145  /// above y=0.
146  UndeformedWall(const double& x0, const double& h): GeomObject(1,2)
147  {
148  X0=x0;
149  H=h;
150  }
151 
152 
153  /// \short Position vector at Lagrangian coordinate zeta
154  void position(const Vector<double>& zeta, Vector<double>& r) const
155  {
156  // Position Vector
157  r[0] = zeta[0]+X0;
158  r[1] = H;
159  }
160 
161 
162  /// \short Parametrised position on object: r(zeta). Evaluated at
163  /// previous timestep. t=0: current time; t>0: previous
164  /// timestep. Calls steady version.
165  void position(const unsigned& t, const Vector<double>& zeta,
166  Vector<double>& r) const
167  {
168  // Use the steady version
169  position(zeta,r);
170 
171  } // end of position
172 
173 
174  /// \short Posn vector and its 1st & 2nd derivatives
175  /// w.r.t. to coordinates:
176  /// \f$ \frac{dR_i}{d \zeta_\alpha}\f$ = drdzeta(alpha,i).
177  /// \f$ \frac{d^2R_i}{d \zeta_\alpha d \zeta_\beta}\f$ =
178  /// ddrdzeta(alpha,beta,i). Evaluated at current time.
179  virtual void d2position(const Vector<double>& zeta,
180  Vector<double>& r,
181  DenseMatrix<double> &drdzeta,
182  RankThreeTensor<double> &ddrdzeta) const
183  {
184  // Position vector
185  r[0] = zeta[0]+X0;
186  r[1] = H;
187 
188  // Tangent vector
189  drdzeta(0,0)=1.0;
190  drdzeta(0,1)=0.0;
191 
192  // Derivative of tangent vector
193  ddrdzeta(0,0,0)=0.0;
194  ddrdzeta(0,0,1)=0.0;
195 
196  } // end of d2position
197 
198  private :
199 
200  /// x position of the undeformed beam's left end.
201  double X0;
202 
203  /// Height of the undeformed wall above y=0.
204  double H;
205 
206 }; //end_of_undeformed_wall
207 
208 
209 
210 
211 
212 
213 
214 //====start_of_physical_parameters=====================
215 /// Namespace for phyical parameters
216 //======================================================
217 namespace Global_Physical_Variables
218 {
219  /// Reynolds number
220  double Re=50.0;
221 
222  /// Womersley = Reynolds times Strouhal
223  double ReSt=50.0;
224 
225  /// Default pressure on the left boundary
226  double P_up=0.0;
227 
228  /// Pseudo-solid mass density
229  double Lambda_sq=0.0;
230 
231  /// Pseudo-solid Poisson ratio
232  double Nu=0.1;
233 
234  /// Traction applied on the fluid at the left (inflow) boundary
235  void prescribed_traction(const double& t,
236  const Vector<double>& x,
237  const Vector<double>& n,
238  Vector<double>& traction)
239  {
240  traction.resize(2);
241  traction[0]=P_up;
242  traction[1]=0.0;
243 
244  } //end traction
245 
246  /// Non-dimensional wall thickness. As in Jensen & Heil (2003) paper.
247  double H=1.0e-2;
248 
249  /// 2nd Piola Kirchhoff pre-stress. As in Jensen & Heil (2003) paper.
250  double Sigma0=1.0e3;
251 
252  /// External pressure
253  double P_ext=0.0;
254 
255  /// \short Load function: Apply a constant external pressure to the wall.
256  /// Note: This is the load without the fluid contribution!
257  /// Fluid load gets added on by FSIWallElement.
258  void load(const Vector<double>& xi, const Vector<double>& x,
259  const Vector<double>& N, Vector<double>& load)
260  {
261  for(unsigned i=0;i<2;i++)
262  {
263  load[i] = -P_ext*N[i];
264  }
265  } //end of load
266 
267 
268  /// \short Fluid structure interaction parameter: Ratio of stresses used for
269  /// non-dimensionalisation of fluid to solid stresses.
270  double Q=1.0e-5;
271 
272 
273 } // end of namespace
274 
275 
276 
277 
278 //====start_of_problem_class==========================================
279 ///Problem class
280 //====================================================================
281 template <class ELEMENT>
282 class FSICollapsibleChannelProblem : public Problem
283 {
284 
285  public :
286 
287 /// \short Constructor: The arguments are the number of elements and
288 /// the lengths of the domain.
289  FSICollapsibleChannelProblem(const unsigned& nup,
290  const unsigned& ncollapsible,
291  const unsigned& ndown,
292  const unsigned& ny,
293  const double& lup,
294  const double& lcollapsible,
295  const double& ldown,
296  const double& ly);
297 
298  /// Destructor (empty)
300 
301  /// Access function for the specific bulk (fluid) mesh
303  {
304  // Upcast from pointer to the Mesh base class to the specific
305  // element type that we're using here.
306  return dynamic_cast<
308  (Bulk_mesh_pt);
309  }
310 
311 
312  /// Access function for the wall mesh
314  {
315  return Wall_mesh_pt;
316 
317  } // end of access to wall mesh
318 
319 
320  /// \short Update the problem specs before solve (empty)
322 
323  /// Update the problem after solve (empty)
325 
326  /// \short Update no slip before Newton convergence check
328  {
329  // Moving wall: No slip; this implies that the velocity needs
330  // to be updated in response to wall motion
331  unsigned ibound=3;
332  unsigned num_nod=bulk_mesh_pt()->nboundary_node(ibound);
333  for (unsigned inod=0;inod<num_nod;inod++)
334  {
335  // Which node are we dealing with?
336  Node* node_pt=bulk_mesh_pt()->boundary_node_pt(ibound,inod);
337 
338  // Apply no slip
339  FSI_functions::apply_no_slip_on_moving_wall(node_pt);
340  }
341  }
342 
343 
344  /// Doc the solution
345  void doc_solution(DocInfo& doc_info,ofstream& trace_file);
346 
347  /// Apply initial conditions
348  void set_initial_condition();
349 
350 private :
351 
352  /// Create the prescribed traction elements on boundary b
353  void create_traction_elements(const unsigned &b,
354  Mesh* const &bulk_mesh_pt,
355  Mesh* const &traction_mesh_pt);
356 
357  /// Delete prescribed traction elements from the surface mesh
358  void delete_traction_elements(Mesh* const &traction_mesh_pt);
359 
360  /// \short Create elements that enforce prescribed boundary motion
361  /// by Lagrange multipliers
362  void create_lagrange_multiplier_elements();
363 
364  /// Delete elements that enforce prescribed boundary motion
365  /// by Lagrange multipliers
366  void delete_lagrange_multiplier_elements();
367 
368  ///Number of elements in the x direction in the upstream part of the channel
369  unsigned Nup;
370 
371  /// \short Number of elements in the x direction in the collapsible part of
372  /// the channel
373  unsigned Ncollapsible;
374 
375  ///Number of elements in the x direction in the downstream part of the channel
376  unsigned Ndown;
377 
378  ///Number of elements across the channel
379  unsigned Ny;
380 
381  ///x-length in the upstream part of the channel
382  double Lup;
383 
384  ///x-length in the collapsible part of the channel
385  double Lcollapsible;
386 
387  ///x-length in the downstream part of the channel
388  double Ldown;
389 
390  ///Transverse length
391  double Ly;
392 
393  /// Pointer to the "bulk" mesh
395 
396  /// \short Pointer to the "surface" mesh that applies the traction at the
397  /// inflow
398  Mesh* Applied_fluid_traction_mesh_pt;
399 
400  /// Pointers to mesh of Lagrange multiplier elements
402 
403  /// Pointer to the "wall" mesh
405 
406  /// Geometric object incarnation of the wall mesh
407  MeshAsGeomObject* Wall_geom_object_pt;
408 
409  ///Pointer to the left control node
410  Node* Left_node_pt;
411 
412  ///Pointer to right control node
413  Node* Right_node_pt;
414 
415  /// Pointer to control node on the wall
416  Node* Wall_node_pt;
417 
418  /// Constitutive law used to determine the mesh deformation
419  ConstitutiveLaw *Constitutive_law_pt;
420 
421 
422 };//end of problem class
423 
424 
425 
426 
427 //=====start_of_constructor======================================
428 /// Constructor for the collapsible channel problem
429 //===============================================================
430 template <class ELEMENT>
432  const unsigned& nup,
433  const unsigned& ncollapsible,
434  const unsigned& ndown,
435  const unsigned& ny,
436  const double& lup,
437  const double& lcollapsible,
438  const double& ldown,
439  const double& ly)
440 {
441  // Store problem parameters
442  Nup=nup;
443  Ncollapsible=ncollapsible;
444  Ndown=ndown;
445  Ny=ny;
446  Lup=lup;
447  Lcollapsible=lcollapsible;
448  Ldown=ldown;
449  Ly=ly;
450 
451 
452  // Overwrite maximum allowed residual to accomodate bad initial guesses
453  Problem::Max_residuals=1000.0;
454 
455  // Allocate the timestepper for the Navier-Stokes equations
456  BDF<2>* fluid_time_stepper_pt=new BDF<2>;
457 
458  // Add the fluid timestepper to the Problem's collection of timesteppers.
459  add_time_stepper_pt(fluid_time_stepper_pt);
460 
461  // Create a dummy Steady timestepper that stores two history values
462  Steady<2>* wall_time_stepper_pt = new Steady<2>;
463 
464  // Add the wall timestepper to the Problem's collection of timesteppers.
465  add_time_stepper_pt(wall_time_stepper_pt);
466 
467  // Geometric object that represents the undeformed wall:
468  // A straight line at height y=ly; starting at x=lup.
469  UndeformedWall* undeformed_wall_pt=new UndeformedWall(lup,ly);
470 
471  //Create the "wall" mesh with FSI Hermite beam elements, passing the
472  //dummy wall timestepper to the constructor
474  (Ncollapsible,Lcollapsible,undeformed_wall_pt,wall_time_stepper_pt);
475 
476  // Build a geometric object (one Lagrangian, two Eulerian coordinates)
477  // from the wall mesh
478  Wall_geom_object_pt=
479  new MeshAsGeomObject(Wall_mesh_pt);
480 
481 
482 
483  //Build bulk (fluid) mesh
484  Bulk_mesh_pt =
486  (nup, ncollapsible, ndown, ny,
487  lup, lcollapsible, ldown, ly,
488  Wall_geom_object_pt,
489  fluid_time_stepper_pt);
490 
491  // Set a non-trivial boundary-layer-squash function...
492  Bulk_mesh_pt->bl_squash_fct_pt() = &BL_Squash::squash_fct;
493 
494  // ... and update the nodal positions accordingly
495  bool update_all_solid_nodes=true;
496  Bulk_mesh_pt->node_update(update_all_solid_nodes);
497  Bulk_mesh_pt->set_lagrangian_nodal_coordinates();
498 
499  // Create "surface mesh" that will contain only the prescribed-traction
500  // elements. The constructor just creates the mesh without
501  // giving it any elements, nodes, etc.
502  Applied_fluid_traction_mesh_pt = new Mesh;
503 
504  // Create prescribed-traction elements from all elements that are
505  // adjacent to boundary 5 (left boundary), but add them to a separate mesh.
506  create_traction_elements(5,Bulk_mesh_pt,Applied_fluid_traction_mesh_pt);
507 
508  // Construct the mesh of elements that enforce prescribed boundary motion
509  // by Lagrange multipliers
510  Lagrange_multiplier_mesh_pt=new SolidMesh;
511  create_lagrange_multiplier_elements();
512 
513  // Add the sub meshes to the problem
514  add_sub_mesh(Bulk_mesh_pt);
515  add_sub_mesh(Applied_fluid_traction_mesh_pt);
516  add_sub_mesh(Wall_mesh_pt);
517  add_sub_mesh(Lagrange_multiplier_mesh_pt);
518 
519  // Combine all submeshes into a single Mesh
520  build_global_mesh();
521 
522  //Set the constitutive law
523  Constitutive_law_pt = new GeneralisedHookean(&Global_Physical_Variables::Nu);
524 
525  // Complete build of fluid mesh
526  //-----------------------------
527 
528  // Loop over the elements to set up element-specific
529  // things that cannot be handled by constructor
530  unsigned n_element=Bulk_mesh_pt->nelement();
531  for(unsigned e=0;e<n_element;e++)
532  {
533 
534  // Upcast from GeneralisedElement to the present element
535  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
536 
537  //Set the Reynolds number
538  el_pt->re_pt() = &Global_Physical_Variables::Re;
539 
540  // Set the Womersley number
541  el_pt->re_st_pt() = &Global_Physical_Variables::ReSt;
542 
543  //Set the constitutive law
544  el_pt->constitutive_law_pt() = Constitutive_law_pt;
545 
546  // Density of pseudo-solid
547  el_pt->lambda_sq_pt()=&Global_Physical_Variables::Lambda_sq;
548 
549  } // end loop over elements
550 
551 
552  // Apply boundary conditions for fluid
553  //------------------------------------
554 
555 
556  //Pin the velocity on the boundaries
557  //x and y-velocities pinned along boundary 0 (bottom boundary) :
558  unsigned ibound=0;
559  unsigned num_nod= bulk_mesh_pt()->nboundary_node(ibound);
560  for (unsigned inod=0;inod<num_nod;inod++)
561  {
562  for(unsigned i=0;i<2;i++)
563  {
564  bulk_mesh_pt()->boundary_node_pt(ibound, inod)->pin(i);
565  }
566  }
567 
568  //x and y-velocities pinned along boundaries 2, 3, 4 (top boundaries) :
569  for(ibound=2;ibound<5;ibound++)
570  {
571  num_nod= bulk_mesh_pt()->nboundary_node(ibound);
572  for (unsigned inod=0;inod<num_nod;inod++)
573  {
574  for(unsigned i=0;i<2;i++)
575  {
576  bulk_mesh_pt()->boundary_node_pt(ibound, inod)->pin(i);
577  }
578  }
579  }
580 
581  //y-velocity pinned along boundary 1 (right boundary):
582  ibound=1;
583  num_nod= bulk_mesh_pt()->nboundary_node(ibound);
584  for (unsigned inod=0;inod<num_nod;inod++)
585  {
586  bulk_mesh_pt()->boundary_node_pt(ibound, inod)->pin(1);
587  }
588 
589 
590  //y-velocity pinned along boundary 5 (left boundary):
591  ibound=5;
592  num_nod= bulk_mesh_pt()->nboundary_node(ibound);
593  for (unsigned inod=0;inod<num_nod;inod++)
594  {
595 
596  bulk_mesh_pt()->boundary_node_pt(ibound, inod)->pin(1);
597 
598  }//end of pin_velocity
599 
600 
601  //Pin the nodal position on all boundaries apart from 3 (the moving wall)
602  for (unsigned ibound=0;ibound<6;ibound++)
603  {
604  if (ibound!=3)
605  {
606  unsigned num_nod= bulk_mesh_pt()->nboundary_node(ibound);
607  for (unsigned inod=0;inod<num_nod;inod++)
608  {
609  for(unsigned i=0;i<2;i++)
610  {
611  dynamic_cast<SolidNode*>(bulk_mesh_pt()->
612  boundary_node_pt(ibound, inod))
613  ->pin_position(i);
614  }
615  }
616  }
617  }
618 
619  // Now pin all nodal positions in the rigid bits
620  unsigned nnod=bulk_mesh_pt()->nnode();
621  for (unsigned j=0;j<nnod;j++)
622  {
623  SolidNode* nod_pt=dynamic_cast<SolidNode*>(bulk_mesh_pt()->node_pt(j));
624  if ((nod_pt->x(0)<=lup)||(nod_pt->x(0)>=(lup+lcollapsible)))
625  {
626  for(unsigned i=0;i<2;i++)
627  {
628  nod_pt->pin_position(i);
629  }
630  }
631  }
632 
633 
634  // Complete build of applied traction elements
635  //--------------------------------------------
636 
637  // Loop over the traction elements to pass pointer to prescribed
638  // traction function
639  unsigned n_el=Applied_fluid_traction_mesh_pt->nelement();
640  for(unsigned e=0;e<n_el;e++)
641  {
642  // Upcast from GeneralisedElement to NavierStokes traction element
643  NavierStokesTractionElement<ELEMENT> *el_pt =
644  dynamic_cast< NavierStokesTractionElement<ELEMENT>*>(
645  Applied_fluid_traction_mesh_pt->element_pt(e));
646 
647  // Set the pointer to the prescribed traction function
648  el_pt->traction_fct_pt() = &Global_Physical_Variables::prescribed_traction;
649 
650  }
651 
652 
653 
654 
655  // Complete build of wall elements
656  //--------------------------------
657 
658  //Loop over the elements to set physical parameters etc.
659  n_element = wall_mesh_pt()->nelement();
660  for(unsigned e=0;e<n_element;e++)
661  {
662  // Upcast to the specific element type
663  FSIHermiteBeamElement *elem_pt =
664  dynamic_cast<FSIHermiteBeamElement*>(wall_mesh_pt()->element_pt(e));
665 
666  // Set physical parameters for each element:
667  elem_pt->sigma0_pt() = &Global_Physical_Variables::Sigma0;
668  elem_pt->h_pt() = &Global_Physical_Variables::H;
669 
670  // Set the load vector for each element
671  elem_pt->load_vector_fct_pt() = &Global_Physical_Variables::load;
672 
673  // Function that specifies the load ratios
674  elem_pt->q_pt() = &Global_Physical_Variables::Q;
675 
676  // Set the undeformed shape for each element
677  elem_pt->undeformed_beam_pt() = undeformed_wall_pt;
678 
679  // The normal on the wall elements as computed by the FSIHermiteElements
680  // points away from the fluid rather than into the fluid (as assumed
681  // by default)
682  elem_pt->set_normal_pointing_out_of_fluid();
683 
684  } // end of loop over elements
685 
686 
687 
688  // Boundary conditions for wall mesh
689  //----------------------------------
690 
691  // Set the boundary conditions: Each end of the beam is fixed in space
692  // Loop over the boundaries (ends of the beam)
693  for(unsigned b=0;b<2;b++)
694  {
695  // Pin displacements in both x and y directions
696  wall_mesh_pt()->boundary_node_pt(b,0)->pin_position(0);
697  wall_mesh_pt()->boundary_node_pt(b,0)->pin_position(1);
698  }
699 
700 
701 
702  //Choose control nodes
703  //---------------------
704 
705  // Left boundary
706  ibound=5;
707  num_nod= bulk_mesh_pt()->nboundary_node(ibound);
708  unsigned control_nod=num_nod/2;
709  Left_node_pt= bulk_mesh_pt()->boundary_node_pt(ibound, control_nod);
710 
711  // Right boundary
712  ibound=1;
713  num_nod= bulk_mesh_pt()->nboundary_node(ibound);
714  control_nod=num_nod/2;
715  Right_node_pt= bulk_mesh_pt()->boundary_node_pt(ibound, control_nod);
716 
717 
718  // Set the pointer to the control node on the wall
719  num_nod= wall_mesh_pt()->nnode();
720  Wall_node_pt=wall_mesh_pt()->node_pt(Ncollapsible/2);
721 
722 
723 
724  // Setup FSI
725  //----------
726 
727  // The velocity of the fluid nodes on the wall (fluid mesh boundary 3)
728  // is set by the wall motion -- hence the no-slip condition needs to be
729  // re-applied whenever a node update is performed for these nodes.
730  // Such tasks may be performed automatically by the auxiliary node update
731  // function specified by a function pointer:
732  ibound=3;
733  num_nod= bulk_mesh_pt()->nboundary_node(ibound);
734  for (unsigned inod=0;inod<num_nod;inod++)
735  {
736  bulk_mesh_pt()->boundary_node_pt(ibound, inod)->
737  set_auxiliary_node_update_fct_pt(
738  FSI_functions::apply_no_slip_on_moving_wall);
739  }
740 
741  // Work out which fluid dofs affect the residuals of the wall elements:
742  // We pass the boundary between the fluid and solid meshes and
743  // pointers to the meshes. The interaction boundary is boundary 3 of the
744  // 2D fluid mesh.
745  FSI_functions::setup_fluid_load_info_for_solid_elements<ELEMENT,2>
746  (this,3,Bulk_mesh_pt,Wall_mesh_pt);
747 
748  // Setup equation numbering scheme
749  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
750 
751 
752 }//end of constructor
753 
754 
755 
756 
757 //====start_of_doc_solution===================================================
758 /// Doc the solution
759 //============================================================================
760 template <class ELEMENT>
762  ofstream& trace_file)
763 {
764  ofstream some_file;
765  char filename[100];
766 
767  // Number of plot points
768  unsigned npts;
769  npts=5;
770 
771  // Output fluid solution
772  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
773  doc_info.number());
774  some_file.open(filename);
775  bulk_mesh_pt()->output(some_file,npts);
776  some_file.close();
777 
778  // Document the wall shape
779  sprintf(filename,"%s/beam%i.dat",doc_info.directory().c_str(),
780  doc_info.number());
781  some_file.open(filename);
782  wall_mesh_pt()->output(some_file,npts);
783  some_file.close();
784 
785  // Loop over all elements do dump out previous solutions
786  // (get the number of previous timesteps available from the wall
787  // timestepper)
788  unsigned nsteps=time_stepper_pt(1)->nprev_values();
789  for (unsigned t=0;t<=nsteps;t++)
790  {
791  sprintf(filename,"%s/wall%i-%i.dat",doc_info.directory().c_str(),
792  doc_info.number(),t);
793  some_file.open(filename);
794  unsigned n_elem=wall_mesh_pt()->nelement();
795  for (unsigned ielem=0;ielem<n_elem;ielem++)
796  {
797  dynamic_cast<FSIHermiteBeamElement*>(wall_mesh_pt()->element_pt(ielem))->
798  output(t,some_file,npts);
799  }
800  some_file.close();
801  } // end of output of previous solutions
802 
803 
804 
805  // Write trace file
806  trace_file << time_pt()->time() << " "
807  << Wall_node_pt->x(1) << " "
808  << Left_node_pt->value(0) << " "
809  << Right_node_pt->value(0) << " "
811  << std::endl;
812 
813 } // end_of_doc_solution
814 
815 
816 
817 
818 
819 //=====start_of_create_traction_elements======================================
820 /// Create the traction elements
821 //============================================================================
822 template <class ELEMENT>
824  const unsigned &b, Mesh* const &bulk_mesh_pt, Mesh* const &traction_mesh_pt)
825 {
826 
827  // How many bulk elements are adjacent to boundary b?
828  unsigned n_element = bulk_mesh_pt->nboundary_element(b);
829 
830  // Loop over the bulk elements adjacent to boundary b?
831  for(unsigned e=0;e<n_element;e++)
832  {
833  // Get pointer to the bulk element that is adjacent to boundary b
834  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>
835  (bulk_mesh_pt->boundary_element_pt(b,e));
836 
837  //What is the index of the face of element e along boundary
838  int face_index = bulk_mesh_pt->face_index_at_boundary(b,e);
839 
840  // Build the corresponding prescribed-traction element
841  NavierStokesTractionElement<ELEMENT>* flux_element_pt =
842  new NavierStokesTractionElement<ELEMENT>(bulk_elem_pt,face_index);
843 
844  //Add the prescribed-traction element to the surface mesh
845  traction_mesh_pt->add_element_pt(flux_element_pt);
846 
847  } //end of loop over bulk elements adjacent to boundary b
848 
849 } // end of create_traction_elements
850 
851 
852 
853 //============start_of_delete_traction_elements==============================
854 /// Delete traction elements and wipe the surface mesh
855 //=======================================================================
856 template<class ELEMENT>
858 delete_traction_elements(Mesh* const &surface_mesh_pt)
859 {
860  // How many surface elements are in the surface mesh
861  unsigned n_element = surface_mesh_pt->nelement();
862 
863  // Loop over the surface elements
864  for(unsigned e=0;e<n_element;e++)
865  {
866  // Kill surface element
867  delete surface_mesh_pt->element_pt(e);
868  }
869 
870  // Wipe the mesh
871  surface_mesh_pt->flush_element_and_node_storage();
872 
873 } // end of delete_traction_elements
874 
875 
876 //============start_of_create_lagrange_multiplier_elements===============
877 /// Create elements that impose the prescribed boundary displacement
878 //=======================================================================
879 template<class ELEMENT>
882 {
883  // Lagrange multiplier elements are located on boundary 3:
884  unsigned b=3;
885 
886  // How many bulk elements are adjacent to boundary b?
887  unsigned n_element = bulk_mesh_pt()->nboundary_element(b);
888 
889  // Loop over the bulk elements adjacent to boundary b?
890  for(unsigned e=0;e<n_element;e++)
891  {
892  // Get pointer to the bulk element that is adjacent to boundary b
893  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
894  bulk_mesh_pt()->boundary_element_pt(b,e));
895 
896  //Find the index of the face of element e along boundary b
897  int face_index = bulk_mesh_pt()->face_index_at_boundary(b,e);
898 
899  // Create new element and add to mesh
900  Lagrange_multiplier_mesh_pt->add_element_pt(
901  new ImposeDisplacementByLagrangeMultiplierElement<ELEMENT>(
902  bulk_elem_pt,face_index));
903  }
904 
905 
906  // Loop over the elements in the Lagrange multiplier element mesh
907  // for elements on the top boundary (boundary 3)
908  n_element=Lagrange_multiplier_mesh_pt->nelement();
909  for(unsigned i=0;i<n_element;i++)
910  {
911  //Cast to a Lagrange multiplier element
912  ImposeDisplacementByLagrangeMultiplierElement<ELEMENT> *el_pt =
913  dynamic_cast<ImposeDisplacementByLagrangeMultiplierElement<ELEMENT>*>
914  (Lagrange_multiplier_mesh_pt->element_pt(i));
915 
916  // Set the GeomObject that defines the boundary shape and set
917  // which bulk boundary we are attached to(needed to extract
918  // the boundary coordinate from the bulk nodes)
919  el_pt->set_boundary_shape_geom_object_pt(Wall_geom_object_pt,b);
920 
921  // Loop over the nodes
922  unsigned nnod=el_pt->nnode();
923  for (unsigned j=0;j<nnod;j++)
924  {
925  Node* nod_pt = el_pt->node_pt(j);
926 
927  // Is the node also on boundary 2 or 4?
928  if ((nod_pt->is_on_boundary(2))||(nod_pt->is_on_boundary(4)))
929  {
930  // How many nodal values were used by the "bulk" element
931  // that originally created this node?
932  unsigned n_bulk_value=el_pt->nbulk_value(j);
933 
934  // The remaining ones are Lagrange multipliers and we pin them.
935  unsigned nval=nod_pt->nvalue();
936  for (unsigned j=n_bulk_value;j<nval;j++)
937  {
938  nod_pt->pin(j);
939  }
940  }
941  }
942  }
943 
944 } // end of create_lagrange_multiplier_elements
945 
946 
947 
948 
949 //====start_of_delete_lagrange_multiplier_elements=======================
950 /// Delete elements that impose the prescribed boundary displacement
951 /// and wipe the associated mesh
952 //=======================================================================
953 template<class ELEMENT>
956 {
957  // How many surface elements are in the surface mesh
958  unsigned n_element = Lagrange_multiplier_mesh_pt->nelement();
959 
960  // Loop over the surface elements
961  for(unsigned e=0;e<n_element;e++)
962  {
963  // Kill surface element
964  delete Lagrange_multiplier_mesh_pt->element_pt(e);
965  }
966 
967  // Wipe the mesh
968  Lagrange_multiplier_mesh_pt->flush_element_and_node_storage();
969 
970 } // end of delete_lagrange_multiplier_elements
971 
972 
973 
974 //====start_of_apply_initial_condition========================================
975 /// Apply initial conditions
976 //============================================================================
977 template <class ELEMENT>
979 {
980  // Check that timestepper is from the BDF family
981  if (time_stepper_pt()->type()!="BDF")
982  {
983  std::ostringstream error_stream;
984  error_stream << "Timestepper has to be from the BDF family!\n"
985  << "You have specified a timestepper from the "
986  << time_stepper_pt()->type() << " family" << std::endl;
987 
988  throw OomphLibError(error_stream.str(),
989  OOMPH_CURRENT_FUNCTION,
990  OOMPH_EXCEPTION_LOCATION);
991  }
992 
993  // Loop over the nodes to set initial guess everywhere
994  unsigned num_nod = bulk_mesh_pt()->nnode();
995  for (unsigned n=0;n<num_nod;n++)
996  {
997  // Get nodal coordinates
998  Vector<double> x(2);
999  x[0]=bulk_mesh_pt()->node_pt(n)->x(0);
1000  x[1]=bulk_mesh_pt()->node_pt(n)->x(1);
1001 
1002  // Assign initial condition: Steady Poiseuille flow
1003  bulk_mesh_pt()->node_pt(n)->set_value(0,6.0*(x[1]/Ly)*(1.0-(x[1]/Ly)));
1004  bulk_mesh_pt()->node_pt(n)->set_value(1,0.0);
1005  }
1006 
1007  // Assign initial values for an impulsive start
1008  bulk_mesh_pt()->assign_initial_values_impulsive();
1009  wall_mesh_pt()->assign_initial_values_impulsive();
1010 
1011 } // end of set_initial_condition
1012 
1013 
1014 
1015 
1016 //============start_of_main====================================================
1017 /// Driver code for a collapsible channel problem with FSI.
1018 /// Presence of command line arguments indicates validation run with
1019 /// coarse resolution and small number of timesteps.
1020 //=============================================================================
1021 int main(int argc, char* argv[])
1022 {
1023 
1024  // Store command line arguments
1025  CommandLineArgs::setup(argc,argv);
1026 
1027  // Reduction in resolution for validation run?
1028  unsigned coarsening_factor=4;
1029  if (CommandLineArgs::Argc>1)
1030  {
1031  coarsening_factor=4;
1032  }
1033 
1034  // Number of elements in the domain
1035  unsigned nup=20/coarsening_factor;
1036  unsigned ncollapsible=40/coarsening_factor;
1037  unsigned ndown=40/coarsening_factor;
1038  unsigned ny=16/coarsening_factor;
1039 
1040  // Length of the domain
1041  double lup=5.0;
1042  double lcollapsible=10.0;
1043  double ldown=10.0;
1044  double ly=1.0;
1045 
1046  // Set external pressure (on the wall stiffness scale).
1048 
1049  // Pressure on the left boundary: This is consistent with steady
1050  // Poiseuille flow
1051  Global_Physical_Variables::P_up=12.0*(lup+lcollapsible+ldown);
1052 
1053 
1054 #ifdef TAYLOR_HOOD
1055 
1056  // Build the problem with QtaylorHoodElements
1058  <PseudoSolidNodeUpdateElement<
1059  QTaylorHoodElement<2>,
1060  QPVDElement<2,3> > >
1061  problem(nup, ncollapsible, ndown, ny,
1062  lup, lcollapsible, ldown, ly);
1063 
1064 #else
1065 
1066  // Build the problem with QCrouzeixRaviartElements
1068  <PseudoSolidNodeUpdateElement<
1069  QCrouzeixRaviartElement<2>,
1070  QPVDElement<2,3> > >
1071  problem(nup, ncollapsible, ndown, ny,
1072  lup, lcollapsible, ldown, ly);
1073 
1074 #endif
1075 
1076  // Timestep. Note: Preliminary runs indicate that the period of
1077  // the oscillation is about 1 so this gives us 40 steps per period.
1078  double dt=1.0/40.0;
1079 
1080  // Initial time for the simulation
1081  double t_min=0.0;
1082 
1083  // Maximum time for simulation
1084  double t_max=3.5;
1085 
1086  // Initialise timestep
1087  problem.time_pt()->time()=t_min;
1088  problem.initialise_dt(dt);
1089 
1090  // Apply initial condition
1091  problem.set_initial_condition();
1092 
1093  //Set output directory
1094  DocInfo doc_info;
1095  doc_info.set_directory("RESLT");
1096 
1097  // Open a trace file
1098  ofstream trace_file;
1099  char filename[100];
1100  sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
1101  trace_file.open(filename);
1102 
1103  // Output the initial condition
1104  problem.doc_solution(doc_info, trace_file);
1105 
1106  // Increment step number
1107  doc_info.number()++;
1108 
1109  // Find number of timesteps (reduced for validation)
1110  unsigned nstep = unsigned((t_max-t_min)/dt);
1111  if (CommandLineArgs::Argc>1)
1112  {
1113  nstep=3;
1114  }
1115 
1116  // Timestepping loop
1117  for (unsigned istep=0;istep<nstep;istep++)
1118  {
1119  // Solve the problem
1120  problem.unsteady_newton_solve(dt);
1121 
1122  // Outpt the solution
1123  problem.doc_solution(doc_info, trace_file);
1124 
1125  // Step number
1126  doc_info.number()++;
1127  }
1128 
1129  // Close trace file.
1130  trace_file.close();
1131 
1132 }//end of main
1133 
void actions_before_newton_convergence_check()
Update no slip before Newton convergence check.
OneDLagrangianMesh< FSIHermiteBeamElement > * wall_mesh_pt()
Access function for the wall mesh.
int main(int argc, char *argv[])
void actions_before_newton_solve()
Update the problem specs before solve (empty)
double Lambda_sq
Pseudo-solid mass density.
void set_initial_condition()
Apply initial conditions.
SolidMesh * Lagrange_multiplier_mesh_pt
Pointers to mesh of Lagrange multiplier elements.
void create_lagrange_multiplier_elements()
Create elements that enforce prescribed boundary motion by Lagrange multipliers.
void load(const Vector< double > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
Load function: Apply a constant external pressure to the wall. Note: This is the load without the flu...
MeshAsGeomObject * Wall_geom_object_pt
Geometric object incarnation of the wall mesh.
double Sigma0
2nd Piola Kirchhoff pre-stress. As in Jensen & Heil (2003) paper.
double Q
Fluid structure interaction parameter: Ratio of stresses used for non-dimensionalisation of fluid to ...
double P_up
Default pressure on the left boundary.
ElasticCollapsibleChannelMesh< ELEMENT > * bulk_mesh_pt()
Access function for the specific bulk (fluid) mesh.
double ReSt
Womersley = Reynolds times Strouhal.
void doc_solution(DocInfo &doc_info, ofstream &trace_file)
Doc the solution.
FSICollapsibleChannelProblem(const unsigned &nup, const unsigned &ncollapsible, const unsigned &ndown, const unsigned &ny, const double &lup, const double &lcollapsible, const double &ldown, const double &ly)
Constructor: The arguments are the number of elements and the lengths of the domain.
virtual void d2position(const Vector< double > &zeta, Vector< double > &r, DenseMatrix< double > &drdzeta, RankThreeTensor< double > &ddrdzeta) const
Posn vector and its 1st & 2nd derivatives w.r.t. to coordinates: = drdzeta(alpha,i). = ddrdzeta(alpha,beta,i). Evaluated at current time.
double Delta
Boundary layer width.
Basic collapsible channel mesh. The mesh is derived from the SimpleRectangularQuadMesh so it's node a...
void actions_after_newton_solve()
Update the problem after solve (empty)
void position(const Vector< double > &zeta, Vector< double > &r) const
Position vector at Lagrangian coordinate zeta.
double Fract_in_BL
Fraction of points in boundary layer.
void create_traction_elements(const unsigned &b, Mesh *const &bulk_mesh_pt, Mesh *const &traction_mesh_pt)
Create the prescribed traction elements on boundary b.
UndeformedWall(const double &x0, const double &h)
Constructor: arguments are the starting point and the height above y=0.
void prescribed_traction(const double &t, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction)
Traction applied on the fluid at the left (inflow) boundary.
double H
Non-dimensional wall thickness. As in Jensen & Heil (2003) paper.
ElasticCollapsibleChannelMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the "bulk" mesh.
ConstitutiveLaw * Constitutive_law_pt
Constitutive law used to determine the mesh deformation.
virtual CollapsibleChannelDomain::BLSquashFctPt & bl_squash_fct_pt()
Function pointer for function that squashes the mesh near the walls. Default trivial mapping (the ide...
double squash_fct(const double &s)
Mapping [0,1] -> [0,1] that re-distributes nodal points across the channel width. ...
void position(const unsigned &t, const Vector< double > &zeta, Vector< double > &r) const
Parametrised position on object: r(zeta). Evaluated at previous timestep. t=0: current time; t>0: pre...
double P_ext
External pressure.
void delete_traction_elements(Mesh *const &traction_mesh_pt)
Delete prescribed traction elements from the surface mesh.
double Nu
Pseudo-solid Poisson ratio.