fsi_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// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 //LIC//
8 //LIC// This library is free software; you can redistribute it and/or
9 //LIC// modify it under the terms of the GNU Lesser General Public
10 //LIC// License as published by the Free Software Foundation; either
11 //LIC// version 2.1 of the License, or (at your option) any later version.
12 //LIC//
13 //LIC// This library is distributed in the hope that it will be useful,
14 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 //LIC// Lesser General Public License for more details.
17 //LIC//
18 //LIC// You should have received a copy of the GNU Lesser General Public
19 //LIC// License along with this library; if not, write to the Free Software
20 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 //LIC// 02110-1301 USA.
22 //LIC//
23 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 //LIC//
25 //LIC//====================================================================
26 #include <iostream>
27 
28 // Generic oomph-lib includes
29 #include "generic.h"
30 #include "navier_stokes.h"
31 #include "beam.h"
32 
33 // The wall mesh
34 #include "meshes/one_d_lagrangian_mesh.h"
35 
36 //Include the fluid mesh
37 #include "meshes/collapsible_channel_mesh.h"
38 
39 using namespace std;
40 
41 using namespace oomph;
42 
43 
44 //==========start_of_BL_Squash =========================================
45 /// Namespace to define the mapping [0,1] -> [0,1] that re-distributes
46 /// nodal points across the channel width.
47 //======================================================================
48 namespace BL_Squash
49 {
50 
51  /// Boundary layer width
52  double Delta=0.1;
53 
54  /// Fraction of points in boundary layer
55  double Fract_in_BL=0.5;
56 
57  /// Mapping [0,1] -> [0,1] that re-distributes
58  /// nodal points across the channel width
59  double squash_fct(const double& s)
60  {
61  // Default return
62  double y=s;
63  if (s<0.5*Fract_in_BL)
64  {
65  y=Delta*2.0*s/Fract_in_BL;
66  }
67  else if (s>1.0-0.5*Fract_in_BL)
68  {
69  y=2.0*Delta/Fract_in_BL*s+1.0-2.0*Delta/Fract_in_BL;
70  }
71  else
72  {
73  y=(1.0-2.0*Delta)/(1.0-Fract_in_BL)*s+
74  (Delta-0.5*Fract_in_BL)/(1.0-Fract_in_BL);
75  }
76 
77  return y;
78  }
79 }// end of BL_Squash
80 
81 
82 
83 
84 
85 
86 /// ///////////////////////////////////////////////////////////////////////
87 /// ///////////////////////////////////////////////////////////////////////
88 /// ///////////////////////////////////////////////////////////////////////
89 
90 
91 //====start_of_underformed_wall============================================
92 /// Undeformed wall is a steady, straight 1D line in 2D space
93 /// \f[ x = X_0 + \zeta \f]
94 /// \f[ y = H \f]
95 //=========================================================================
96 class UndeformedWall : public GeomObject
97 {
98 
99 public:
100 
101  /// Constructor: arguments are the starting point and the height
102  /// above y=0.
103  UndeformedWall(const double& x0, const double& h): GeomObject(1,2)
104  {
105  X0=x0;
106  H=h;
107  }
108 
109 
110  /// Position vector at Lagrangian coordinate zeta
111  void position(const Vector<double>& zeta, Vector<double>& r) const
112  {
113  // Position Vector
114  r[0] = zeta[0]+X0;
115  r[1] = H;
116  }
117 
118 
119  /// Parametrised position on object: r(zeta). Evaluated at
120  /// previous timestep. t=0: current time; t>0: previous
121  /// timestep. Calls steady version.
122  void position(const unsigned& t, const Vector<double>& zeta,
123  Vector<double>& r) const
124  {
125  // Use the steady version
126  position(zeta,r);
127 
128  } // end of position
129 
130 
131  /// Posn vector and its 1st & 2nd derivatives
132  /// w.r.t. to coordinates:
133  /// \f$ \frac{dR_i}{d \zeta_\alpha}\f$ = drdzeta(alpha,i).
134  /// \f$ \frac{d^2R_i}{d \zeta_\alpha d \zeta_\beta}\f$ =
135  /// ddrdzeta(alpha,beta,i). Evaluated at current time.
136  void d2position(const Vector<double>& zeta,
137  Vector<double>& r,
138  DenseMatrix<double> &drdzeta,
139  RankThreeTensor<double> &ddrdzeta) const
140  {
141  // Position vector
142  r[0] = zeta[0]+X0;
143  r[1] = H;
144 
145  // Tangent vector
146  drdzeta(0,0)=1.0;
147  drdzeta(0,1)=0.0;
148 
149  // Derivative of tangent vector
150  ddrdzeta(0,0,0)=0.0;
151  ddrdzeta(0,0,1)=0.0;
152 
153  } // end of d2position
154 
155  private :
156 
157  /// x position of the undeformed beam's left end.
158  double X0;
159 
160  /// Height of the undeformed wall above y=0.
161  double H;
162 
163 }; //end_of_undeformed_wall
164 
165 
166 
167 
168 
169 
170 
171 //====start_of_physical_parameters=====================
172 /// Namespace for phyical parameters
173 //======================================================
175 {
176  /// Reynolds number
177  double Re=50.0;
178 
179  /// Womersley = Reynolds times Strouhal
180  double ReSt=50.0;
181 
182  /// Default pressure on the left boundary
183  double P_up=0.0;
184 
185  /// Traction applied on the fluid at the left (inflow) boundary
186  void prescribed_traction(const double& t,
187  const Vector<double>& x,
188  const Vector<double>& n,
189  Vector<double>& traction)
190  {
191  traction.resize(2);
192  traction[0]=P_up;
193  traction[1]=0.0;
194 
195  } //end traction
196 
197  /// Non-dimensional wall thickness. As in Jensen & Heil (2003) paper.
198  double H=1.0e-2;
199 
200  /// 2nd Piola Kirchhoff pre-stress. As in Jensen & Heil (2003) paper.
201  double Sigma0=1.0e3;
202 
203  /// External pressure
204  double P_ext=0.0;
205 
206  /// Load function: Apply a constant external pressure to the wall.
207  /// Note: This is the load without the fluid contribution!
208  /// Fluid load gets added on by FSIWallElement.
209  void load(const Vector<double>& xi, const Vector<double>& x,
210  const Vector<double>& N, Vector<double>& load)
211  {
212  for(unsigned i=0;i<2;i++)
213  {
214  load[i] = -P_ext*N[i];
215  }
216  } //end of load
217 
218 
219  /// Fluid structure interaction parameter: Ratio of stresses used for
220  /// non-dimensionalisation of fluid to solid stresses.
221  double Q=1.0e-5;
222 
223 
224 } // end of namespace
225 
226 
227 
228 
229 //====start_of_problem_class==========================================
230 /// Problem class
231 //====================================================================
232 template <class ELEMENT>
233 class FSICollapsibleChannelProblem : public Problem
234 {
235 
236  public :
237 
238 /// Constructor: The arguments are the number of elements and
239 /// the lengths of the domain.
240  FSICollapsibleChannelProblem(const unsigned& nup,
241  const unsigned& ncollapsible,
242  const unsigned& ndown,
243  const unsigned& ny,
244  const double& lup,
245  const double& lcollapsible,
246  const double& ldown,
247  const double& ly);
248 
249  /// Destructor (empty)
251 
252 
253 #ifdef MACRO_ELEMENT_NODE_UPDATE
254 
255  /// Access function for the specific bulk (fluid) mesh
256  MacroElementNodeUpdateCollapsibleChannelMesh<ELEMENT>* bulk_mesh_pt()
257  {
258  // Upcast from pointer to the Mesh base class to the specific
259  // element type that we're using here.
260  return dynamic_cast<
261  MacroElementNodeUpdateCollapsibleChannelMesh<ELEMENT>*>
262  (Bulk_mesh_pt);
263  }
264 
265 #else
266 
267  /// Access function for the specific bulk (fluid) mesh
268  AlgebraicCollapsibleChannelMesh<ELEMENT>* bulk_mesh_pt()
269  {
270  // Upcast from pointer to the Mesh base class to the specific
271  // element type that we're using here.
272  return dynamic_cast<
273  AlgebraicCollapsibleChannelMesh<ELEMENT>*>
274  (Bulk_mesh_pt);
275  }
276 
277 #endif
278 
279 
280  /// Access function for the wall mesh
281  OneDLagrangianMesh<FSIHermiteBeamElement>* wall_mesh_pt()
282  {
283  return Wall_mesh_pt;
284 
285  } // end of access to wall mesh
286 
287 
288  /// Update the problem specs before solve (empty)
290 
291  /// Update the problem after solve (empty)
293 
294  /// Update before checking Newton convergence: Update the
295  /// nodal positions in the fluid mesh in response to possible
296  /// changes in the wall shape
298  {
299  Bulk_mesh_pt->node_update();
300  }
301 
302  /// Doc the solution
303  void doc_solution(DocInfo& doc_info,ofstream& trace_file);
304 
305  /// Apply initial conditions
306  void set_initial_condition();
307 
308 private :
309 
310  /// Create the prescribed traction elements on boundary b
311  void create_traction_elements(const unsigned &b,
312  Mesh* const &bulk_mesh_pt,
313  Mesh* const &traction_mesh_pt);
314 
315  /// Number of elements in the x direction in the upstream part of the channel
316  unsigned Nup;
317 
318  /// Number of elements in the x direction in the collapsible part of
319  /// the channel
320  unsigned Ncollapsible;
321 
322  /// Number of elements in the x direction in the downstream part of the channel
323  unsigned Ndown;
324 
325  /// Number of elements across the channel
326  unsigned Ny;
327 
328  /// x-length in the upstream part of the channel
329  double Lup;
330 
331  /// x-length in the collapsible part of the channel
332  double Lcollapsible;
333 
334  /// x-length in the downstream part of the channel
335  double Ldown;
336 
337  /// Transverse length
338  double Ly;
339 
340 #ifdef MACRO_ELEMENT_NODE_UPDATE
341 
342  /// Pointer to the "bulk" mesh
343  MacroElementNodeUpdateCollapsibleChannelMesh<ELEMENT>* Bulk_mesh_pt;
344 
345 #else
346 
347  /// Pointer to the "bulk" mesh
348  AlgebraicCollapsibleChannelMesh<ELEMENT>* Bulk_mesh_pt;
349 
350 #endif
351 
352 
353  /// Pointer to the "surface" mesh that applies the traction at the
354  /// inflow
356 
357  /// Pointer to the "wall" mesh
358  OneDLagrangianMesh<FSIHermiteBeamElement>* Wall_mesh_pt;
359 
360  /// Pointer to the left control node
362 
363  /// Pointer to right control node
365 
366  /// Pointer to control node on the wall
368 
369 };//end of problem class
370 
371 
372 
373 
374 //=====start_of_constructor======================================
375 /// Constructor for the collapsible channel problem
376 //===============================================================
377 template <class ELEMENT>
379  const unsigned& nup,
380  const unsigned& ncollapsible,
381  const unsigned& ndown,
382  const unsigned& ny,
383  const double& lup,
384  const double& lcollapsible,
385  const double& ldown,
386  const double& ly)
387 {
388  // Store problem parameters
389  Nup=nup;
390  Ncollapsible=ncollapsible;
391  Ndown=ndown;
392  Ny=ny;
393  Lup=lup;
394  Lcollapsible=lcollapsible;
395  Ldown=ldown;
396  Ly=ly;
397 
398 
399  // Overwrite maximum allowed residual to accomodate bad initial guesses
400  Problem::Max_residuals=1000.0;
401 
402  // Allocate the timestepper -- this constructs the Problem's
403  // time object with a sufficient amount of storage to store the
404  // previous timsteps.
405  add_time_stepper_pt(new BDF<2>);
406 
407  // Geometric object that represents the undeformed wall:
408  // A straight line at height y=ly; starting at x=lup.
409  UndeformedWall* undeformed_wall_pt=new UndeformedWall(lup,ly);
410 
411  //Create the "wall" mesh with FSI Hermite elements
412  Wall_mesh_pt = new OneDLagrangianMesh<FSIHermiteBeamElement>
413  //(2*Ncollapsible+5,Lcollapsible,undeformed_wall_pt);
414  (Ncollapsible,Lcollapsible,undeformed_wall_pt);
415 
416 
417 
418  // Build a geometric object (one Lagrangian, two Eulerian coordinates)
419  // from the wall mesh
420  MeshAsGeomObject* wall_geom_object_pt=
421  new MeshAsGeomObject(Wall_mesh_pt);
422 
423 
424 #ifdef MACRO_ELEMENT_NODE_UPDATE
425 
426  //Build bulk (fluid) mesh
427  Bulk_mesh_pt =
428  new MacroElementNodeUpdateCollapsibleChannelMesh<ELEMENT>
429  (nup, ncollapsible, ndown, ny,
430  lup, lcollapsible, ldown, ly,
431  wall_geom_object_pt,
432  time_stepper_pt());
433 
434  // Set a non-trivial boundary-layer-squash function...
435  Bulk_mesh_pt->bl_squash_fct_pt() = &BL_Squash::squash_fct;
436 
437  // ... and update the nodal positions accordingly
438  Bulk_mesh_pt->node_update();
439 
440 #else
441 
442  //Build bulk (fluid) mesh
443  Bulk_mesh_pt =
444  new AlgebraicCollapsibleChannelMesh<ELEMENT>
445  (nup, ncollapsible, ndown, ny,
446  lup, lcollapsible, ldown, ly,
447  wall_geom_object_pt,
449  time_stepper_pt());
450 
451 #endif
452 
453 
454  // Create "surface mesh" that will contain only the prescribed-traction
455  // elements. The constructor just creates the mesh without
456  // giving it any elements, nodes, etc.
457  Applied_fluid_traction_mesh_pt = new Mesh;
458 
459  // Create prescribed-traction elements from all elements that are
460  // adjacent to boundary 5 (left boundary), but add them to a separate mesh.
461  create_traction_elements(5,Bulk_mesh_pt,Applied_fluid_traction_mesh_pt);
462 
463  // Add the sub meshes to the problem
464  add_sub_mesh(Bulk_mesh_pt);
465  add_sub_mesh(Applied_fluid_traction_mesh_pt);
466  add_sub_mesh(Wall_mesh_pt);
467 
468  // Combine all submeshes into a single Mesh
469  build_global_mesh();
470 
471 
472  // Complete build of fluid mesh
473  //-----------------------------
474 
475  // Loop over the elements to set up element-specific
476  // things that cannot be handled by constructor
477  unsigned n_element=Bulk_mesh_pt->nelement();
478  for(unsigned e=0;e<n_element;e++)
479  {
480  // Upcast from GeneralisedElement to the present element
481  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
482 
483  //Set the Reynolds number
484  el_pt->re_pt() = &Global_Physical_Variables::Re;
485 
486  // Set the Womersley number
487  el_pt->re_st_pt() = &Global_Physical_Variables::ReSt;
488 
489  } // end loop over elements
490 
491 
492 
493  // Apply boundary conditions for fluid
494  //------------------------------------
495 
496  //Pin the velocity on the boundaries
497  //x and y-velocities pinned along boundary 0 (bottom boundary) :
498  unsigned ibound=0;
499  unsigned num_nod= bulk_mesh_pt()->nboundary_node(ibound);
500  for (unsigned inod=0;inod<num_nod;inod++)
501  {
502  for(unsigned i=0;i<2;i++)
503  {
504  bulk_mesh_pt()->boundary_node_pt(ibound, inod)->pin(i);
505  }
506  }
507 
508  //x and y-velocities pinned along boundaries 2, 3, 4 (top boundaries) :
509  for(ibound=2;ibound<5;ibound++)
510  {
511  num_nod= bulk_mesh_pt()->nboundary_node(ibound);
512  for (unsigned inod=0;inod<num_nod;inod++)
513  {
514  for(unsigned i=0;i<2;i++)
515  {
516  bulk_mesh_pt()->boundary_node_pt(ibound, inod)->pin(i);
517  }
518  }
519  }
520 
521  //y-velocity pinned along boundary 1 (right boundary):
522  ibound=1;
523  num_nod= bulk_mesh_pt()->nboundary_node(ibound);
524  for (unsigned inod=0;inod<num_nod;inod++)
525  {
526  bulk_mesh_pt()->boundary_node_pt(ibound, inod)->pin(1);
527  }
528 
529 
530  //y-velocity pinned along boundary 5 (left boundary):
531  ibound=5;
532  num_nod= bulk_mesh_pt()->nboundary_node(ibound);
533  for (unsigned inod=0;inod<num_nod;inod++)
534  {
535  bulk_mesh_pt()->boundary_node_pt(ibound, inod)->pin(1);
536  }
537 //end of pin_velocity
538 
539  // Complete build of applied traction elements
540  //--------------------------------------------
541 
542  // Loop over the traction elements to pass pointer to prescribed
543  // traction function
544  unsigned n_el=Applied_fluid_traction_mesh_pt->nelement();
545  for(unsigned e=0;e<n_el;e++)
546  {
547  // Upcast from GeneralisedElement to NavierStokes traction element
548  NavierStokesTractionElement<ELEMENT> *el_pt =
549  dynamic_cast< NavierStokesTractionElement<ELEMENT>*>(
550  Applied_fluid_traction_mesh_pt->element_pt(e));
551 
552  // Set the pointer to the prescribed traction function
553  el_pt->traction_fct_pt() = &Global_Physical_Variables::prescribed_traction;
554  }
555 
556 
557 
558 
559  // Complete build of wall elements
560  //--------------------------------
561 
562  //Loop over the elements to set physical parameters etc.
563  n_element = wall_mesh_pt()->nelement();
564  for(unsigned e=0;e<n_element;e++)
565  {
566  // Upcast to the specific element type
567  FSIHermiteBeamElement *elem_pt =
568  dynamic_cast<FSIHermiteBeamElement*>(wall_mesh_pt()->element_pt(e));
569 
570  // Set physical parameters for each element:
571  elem_pt->sigma0_pt() = &Global_Physical_Variables::Sigma0;
572  elem_pt->h_pt() = &Global_Physical_Variables::H;
573 
574  // Set the load vector for each element
575  elem_pt->load_vector_fct_pt() = &Global_Physical_Variables::load;
576 
577  // Function that specifies the load ratios
578  elem_pt->q_pt() = &Global_Physical_Variables::Q;
579 
580  // Set the undeformed shape for each element
581  elem_pt->undeformed_beam_pt() = undeformed_wall_pt;
582 
583 
584  // The normal on the wall elements as computed by the FSIHermiteElements
585  // points away from the fluid rather than into the fluid (as assumed
586  // by default)
587  elem_pt->set_normal_pointing_out_of_fluid();
588 
589  } // end of loop over elements
590 
591 
592 
593  // Boundary conditions for wall mesh
594  //----------------------------------
595 
596  // Set the boundary conditions: Each end of the beam is fixed in space
597  // Loop over the boundaries (ends of the beam)
598  for(unsigned b=0;b<2;b++)
599  {
600  // Pin displacements in both x and y directions
601  wall_mesh_pt()->boundary_node_pt(b,0)->pin_position(0);
602  wall_mesh_pt()->boundary_node_pt(b,0)->pin_position(1);
603  }
604 
605 
606 
607 
608  //Choose control nodes
609  //---------------------
610 
611  // Left boundary
612  ibound=5;
613  num_nod= bulk_mesh_pt()->nboundary_node(ibound);
614  unsigned control_nod=num_nod/2;
615  Left_node_pt= bulk_mesh_pt()->boundary_node_pt(ibound, control_nod);
616 
617  // Right boundary
618  ibound=1;
619  num_nod= bulk_mesh_pt()->nboundary_node(ibound);
620  control_nod=num_nod/2;
621  Right_node_pt= bulk_mesh_pt()->boundary_node_pt(ibound, control_nod);
622 
623 
624  // Set the pointer to the control node on the wall
625  Wall_node_pt=wall_mesh_pt()->node_pt(Ncollapsible/2);
626 
627 
628 
629 
630  // Setup FSI
631  //----------
632 
633  // The velocity of the fluid nodes on the wall (fluid mesh boundary 3)
634  // is set by the wall motion -- hence the no-slip condition must be
635  // re-applied whenever a node update is performed for these nodes.
636  // Such tasks may be performed automatically by the auxiliary node update
637  // function specified by a function pointer:
638  ibound=3;
639  num_nod= bulk_mesh_pt()->nboundary_node(ibound);
640  for (unsigned inod=0;inod<num_nod;inod++)
641  {
642  bulk_mesh_pt()->boundary_node_pt(ibound, inod)->
643  set_auxiliary_node_update_fct_pt(
644  FSI_functions::apply_no_slip_on_moving_wall);
645  }
646 
647 
648  // Work out which fluid dofs affect the residuals of the wall elements:
649  // We pass the boundary between the fluid and solid meshes and
650  // pointers to the meshes. The interaction boundary is boundary 3 of the
651  // 2D fluid mesh.
652  FSI_functions::setup_fluid_load_info_for_solid_elements<ELEMENT,2>
653  (this,3,Bulk_mesh_pt,Wall_mesh_pt);
654 
655  // Setup equation numbering scheme
656  cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
657 
658 
659 }//end of constructor
660 
661 
662 
663 
664 //====start_of_doc_solution===================================================
665 /// Doc the solution
666 //============================================================================
667 template <class ELEMENT>
669  ofstream& trace_file)
670 {
671 
672 
673 #ifdef MACRO_ELEMENT_NODE_UPDATE
674 
675  // Doc fsi
676  if (CommandLineArgs::Argc>1)
677  {
678  FSI_functions::doc_fsi<MacroElementNodeUpdateNode>(Bulk_mesh_pt,Wall_mesh_pt,doc_info);
679  }
680 
681 #else
682 
683  // Doc fsi
684  if (CommandLineArgs::Argc>1)
685  {
686  FSI_functions::doc_fsi<AlgebraicNode>(Bulk_mesh_pt,Wall_mesh_pt,doc_info);
687  }
688 
689 #endif
690 
691  ofstream some_file;
692  char filename[100];
693 
694  // Number of plot points
695  unsigned npts;
696  npts=5;
697 
698  // Output fluid solution
699  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
700  doc_info.number());
701  some_file.open(filename);
702  bulk_mesh_pt()->output(some_file,npts);
703  some_file.close();
704 
705  // Document the wall shape
706  sprintf(filename,"%s/beam%i.dat",doc_info.directory().c_str(),
707  doc_info.number());
708  some_file.open(filename);
709  wall_mesh_pt()->output(some_file,npts);
710  some_file.close();
711 
712 
713  // Write trace file
714  trace_file << time_pt()->time() << " "
715  << Wall_node_pt->x(1) << " "
716  << Left_node_pt->value(0) << " "
717  << Right_node_pt->value(0) << " "
719  << std::endl;
720 
721 } // end_of_doc_solution
722 
723 
724 
725 
726 //=====start_of_create_traction_elements======================================
727 /// Create the traction elements
728 //============================================================================
729 template <class ELEMENT>
731  const unsigned &b, Mesh* const &bulk_mesh_pt, Mesh* const &traction_mesh_pt)
732 {
733 
734  // How many bulk elements are adjacent to boundary b?
735  unsigned n_element = bulk_mesh_pt->nboundary_element(b);
736 
737  // Loop over the bulk elements adjacent to boundary b?
738  for(unsigned e=0;e<n_element;e++)
739  {
740  // Get pointer to the bulk element that is adjacent to boundary b
741  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>
742  (bulk_mesh_pt->boundary_element_pt(b,e));
743 
744  //What is the index of the face of element e along boundary b
745  int face_index = bulk_mesh_pt->face_index_at_boundary(b,e);
746 
747  // Build the corresponding prescribed-traction element
748  NavierStokesTractionElement<ELEMENT>* flux_element_pt =
749  new NavierStokesTractionElement<ELEMENT>(bulk_elem_pt,face_index);
750 
751  //Add the prescribed-traction element to the surface mesh
752  traction_mesh_pt->add_element_pt(flux_element_pt);
753 
754  } //end of loop over bulk elements adjacent to boundary b
755 
756 } // end of create_traction_elements
757 
758 
759 
760 
761 //====start_of_apply_initial_condition========================================
762 /// Apply initial conditions
763 //============================================================================
764 template <class ELEMENT>
766 {
767  // Check that timestepper is from the BDF family
768  if (time_stepper_pt()->type()!="BDF")
769  {
770  std::ostringstream error_stream;
771  error_stream << "Timestepper has to be from the BDF family!\n"
772  << "You have specified a timestepper from the "
773  << time_stepper_pt()->type() << " family" << std::endl;
774 
775  throw OomphLibError(error_stream.str(),
776  OOMPH_CURRENT_FUNCTION,
777  OOMPH_EXCEPTION_LOCATION);
778  }
779 
780  // Update the mesh
781  bulk_mesh_pt()->node_update();
782 
783  // Loop over the nodes to set initial guess everywhere
784  unsigned num_nod = bulk_mesh_pt()->nnode();
785  for (unsigned n=0;n<num_nod;n++)
786  {
787  // Get nodal coordinates
788  Vector<double> x(2);
789  x[0]=bulk_mesh_pt()->node_pt(n)->x(0);
790  x[1]=bulk_mesh_pt()->node_pt(n)->x(1);
791 
792  // Assign initial condition: Steady Poiseuille flow
793  bulk_mesh_pt()->node_pt(n)->set_value(0,6.0*(x[1]/Ly)*(1.0-(x[1]/Ly)));
794  bulk_mesh_pt()->node_pt(n)->set_value(1,0.0);
795  }
796 
797  // Assign initial values for an impulsive start
798  bulk_mesh_pt()->assign_initial_values_impulsive();
799 
800 } // end of set_initial_condition
801 
802 
803 
804 
805 
806 //============start_of_main====================================================
807 /// Driver code for a collapsible channel problem with FSI.
808 /// Presence of command line arguments indicates validation run with
809 /// coarse resolution and small number of timesteps.
810 //=============================================================================
811 int main(int argc, char* argv[])
812 {
813 
814  // Store command line arguments
815  CommandLineArgs::setup(argc,argv);
816 
817  // Reduction in resolution for validation run?
818  unsigned coarsening_factor=1;
819  if (CommandLineArgs::Argc>1)
820  {
821  coarsening_factor=4;
822  }
823 
824  // Number of elements in the domain
825  unsigned nup=20/coarsening_factor;
826  unsigned ncollapsible=40/coarsening_factor;
827  unsigned ndown=40/coarsening_factor;
828  unsigned ny=16/coarsening_factor;
829 
830  // Length of the domain
831  double lup=5.0;
832  double lcollapsible=10.0;
833  double ldown=10.0;
834  double ly=1.0;
835 
836  // Set external pressure (on the wall stiffness scale).
838 
839  // Pressure on the left boundary: This is consistent with steady
840  // Poiseuille flow
841  Global_Physical_Variables::P_up=12.0*(lup+lcollapsible+ldown);
842 
843 #ifdef MACRO_ELEMENT_NODE_UPDATE
844 
845 #ifdef TAYLOR_HOOD
846 
847  // Build the problem with QTaylorHoodElements
849  <MacroElementNodeUpdateElement<QTaylorHoodElement<2> > >
850  problem(nup, ncollapsible, ndown, ny,
851  lup, lcollapsible, ldown, ly);
852 
853 #else
854 
855  // Build the problem with QCrouzeixRaviartElements
857  <MacroElementNodeUpdateElement<QCrouzeixRaviartElement<2> > >
858  problem(nup, ncollapsible, ndown, ny,
859  lup, lcollapsible, ldown, ly);
860 
861 #endif
862 
863 #else
864 
865 #ifdef TAYLOR_HOOD
866 
867  // Build the problem with QCrouzeixRaviartElements
869  <AlgebraicElement<QTaylorHoodElement<2> > >
870  problem(nup, ncollapsible, ndown, ny,
871  lup, lcollapsible, ldown, ly);
872 
873 #else
874 
875  // Build the problem with QCrouzeixRaviartElements
877  <AlgebraicElement<QCrouzeixRaviartElement<2> > >
878  problem(nup, ncollapsible, ndown, ny,
879  lup, lcollapsible, ldown, ly);
880 
881 #endif
882 
883 #endif
884  // Timestep. Note: Preliminary runs indicate that the period of
885  // the oscillation is about 1 so this gives us 40 steps per period.
886  double dt=1.0/40.0;
887 
888  // Initial time for the simulation
889  double t_min=0.0;
890 
891  // Maximum time for simulation
892  double t_max=3.5;
893 
894  // Initialise timestep
895  problem.time_pt()->time()=t_min;
896  problem.initialise_dt(dt);
897 
898  // Apply initial condition
899  problem.set_initial_condition();
900 
901  //Set output directory
902  DocInfo doc_info;
903  doc_info.set_directory("RESLT");
904 
905  // Open a trace file
906  ofstream trace_file;
907  char filename[100];
908  sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
909  trace_file.open(filename);
910 
911  // Output the initial solution
912  problem.doc_solution(doc_info, trace_file);
913 
914  // Increment step number
915  doc_info.number()++;
916 
917  // Find number of timesteps (reduced for validation)
918  unsigned nstep = unsigned((t_max-t_min)/dt);
919  if (CommandLineArgs::Argc>1)
920  {
921  nstep=3;
922  }
923 
924  // Timestepping loop
925  for (unsigned istep=0;istep<nstep;istep++)
926  {
927  // Solve the problem
928  problem.unsteady_newton_solve(dt);
929 
930  // Outpt the solution
931  problem.doc_solution(doc_info, trace_file);
932 
933  // Step number
934  doc_info.number()++;
935  }
936 
937 
938  // Close trace file.
939  trace_file.close();
940 
941 }//end of main
942 
Mesh * Applied_fluid_traction_mesh_pt
Pointer to the "surface" mesh that applies the traction at the inflow.
double Ldown
x-length in the downstream part of the channel
Node * Left_node_pt
Pointer to the left control node.
Node * Right_node_pt
Pointer to right control node.
unsigned Nup
Number of elements in the x direction in the upstream part of the channel.
OneDLagrangianMesh< FSIHermiteBeamElement > * Wall_mesh_pt
Pointer to the "wall" mesh.
unsigned Ncollapsible
Number of elements in the x direction in the collapsible part of the channel.
void actions_after_newton_solve()
Update the problem after solve (empty)
unsigned Ny
Number of elements across the channel.
MacroElementNodeUpdateCollapsibleChannelMesh< ELEMENT > * bulk_mesh_pt()
Access function for the specific bulk (fluid) mesh.
unsigned Ndown
Number of elements in the x direction in the downstream part of the channel.
AlgebraicCollapsibleChannelMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the "bulk" mesh.
double Lup
x-length in the upstream part of the channel
AlgebraicCollapsibleChannelMesh< ELEMENT > * bulk_mesh_pt()
Access function for the specific bulk (fluid) mesh.
MacroElementNodeUpdateCollapsibleChannelMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the "bulk" mesh.
~FSICollapsibleChannelProblem()
Destructor (empty)
void actions_before_newton_convergence_check()
Update before checking Newton convergence: Update the nodal positions in the fluid mesh in response t...
void actions_before_newton_solve()
Update the problem specs before solve (empty)
Node * Wall_node_pt
Pointer to control node on the wall.
double Lcollapsible
x-length in the collapsible part of the channel
OneDLagrangianMesh< FSIHermiteBeamElement > * wall_mesh_pt()
Access function for the wall mesh.
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.
void set_initial_condition()
Apply initial conditions.
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.
void doc_solution(DocInfo &doc_info, ofstream &trace_file)
Doc the solution.
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
double H
Height of the undeformed wall above y=0.
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,...
double X0
x position of the undeformed beam's left end.
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...
void position(const Vector< double > &zeta, Vector< double > &r) const
Position vector at Lagrangian coordinate zeta.
UndeformedWall(const double &x0, const double &h)
Constructor: arguments are the starting point and the height above y=0.
int main(int argc, char *argv[])
Driver code for a collapsible channel problem with FSI. Presence of command line arguments indicates ...
Namespace to define the mapping [0,1] -> [0,1] that re-distributes nodal points across the channel wi...
double squash_fct(const double &s)
Mapping [0,1] -> [0,1] that re-distributes nodal points across the channel width.
double Delta
Boundary layer width.
double Fract_in_BL
Fraction of points in boundary layer.
Namespace for phyical parameters.
double P_ext
External pressure.
double ReSt
Womersley = Reynolds times Strouhal.
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.
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...
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.
double H
Non-dimensional wall thickness. As in Jensen & Heil (2003) paper.