fsi_osc_ring.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 // Driver for 2D Navier Stokes flow interacting with an elastic ring
27 
28 // Oomph-lib include files
29 #include "generic.h"
30 #include "navier_stokes.h"
31 #include "beam.h"
32 
33 //Need to include templated meshes, so that all functions
34 //are instantiated for our particular element types.
35 #include "meshes/quarter_circle_sector_mesh.h"
36 #include "meshes/one_d_lagrangian_mesh.h"
37 
38 using namespace std;
39 
40 using namespace oomph;
41 
42 //===========================================================================
43 /// Namespace for physical parameters
44 //===========================================================================
46 {
47 
48  // Independent parameters:
49  //------------------------
50 
51  /// Square of Womersly number (a frequency parameter)
52  double Alpha_sq=50.0;
53 
54  /// Density ratio of the solid and the fluid
55  double Density_ratio=1.0;
56 
57  /// External Pressure
58  double Pext=0.0;
59 
60  /// Poisson ratio
61  double Nu=0.49;
62 
63  /// Nondimensional thickness of the beam
64  double H=0.05;
65 
66  /// Perturbation pressure
67  double Pcos=0.0;
68 
69 
70  // Dependent parameters:
71  //----------------------
72 
73  /// Reynolds number
74  double Re;
75 
76  /// Reynolds x Strouhal number
77  double ReSt;
78 
79  /// Timescale ratio (non-dimensation density)
80  double Lambda_sq;
81 
82  /// Stress ratio
83  double Q;
84 
85  /// Set the parameters that are used in the code as a function
86  /// of the Womersley number, the density ratio and H
87  void set_params()
88  {
89  cout << "\n\n======================================================" <<std::endl;
90  cout << "\nSetting parameters. \n\n";
91  cout << "Prescribed: Square of Womersley number: Alpha_sq = "
92  << Alpha_sq << std::endl;
93  cout << " Density ratio: Density_ratio = "
94  << Density_ratio << std::endl;
95  cout << " Wall thickness: H = "
96  << H << std::endl;
97  cout << " Poisson ratio: Nu = "
98  << Nu << std::endl;
99  cout << " Pressure perturbation: Pcos = "
100  << Pcos << std::endl;
101 
102 
103  Q=1.0/12.0*pow(H,3)/Alpha_sq;
104  cout << "\nDependent: Stress ratio: Q = "
105  << Q << std::endl;
106 
107  Lambda_sq=1.0/12.0*pow(H,3)*Density_ratio;
108  cout << " Timescale ratio: Lambda_sq = "
109  << Lambda_sq << std::endl;
110 
111  Re=Alpha_sq;
112  cout << " Reynolds number: Re = "
113  << Re << std::endl;
114 
115  ReSt=Re;
116  cout << " Womersley number: ReSt = "
117  << ReSt << std::endl;
118  cout << "\n======================================================\n\n"
119  <<std::endl;
120 
121  }
122 
123  /// Non-FSI load function, a constant external pressure plus
124  /// a (small) sinusoidal perturbation of wavenumber two.
125  void pcos_load(const Vector<double>& xi, const Vector<double> &x,
126  const Vector<double>& N, Vector<double>& load)
127  {
128  for(unsigned i=0;i<2;i++)
129  {load[i] = (Pext - Pcos*cos(2.0*xi[0]))*N[i];}
130  }
131 
132 }
133 
134 
135 //======================================================================
136 /// FSI Ring problem: a fluid-structure interaction problem in which
137 /// a viscous fluid bounded by an initially circular beam is set into motion
138 /// by a small sinusoidal perturbation of the beam (the domain boundary).
139 //======================================================================
140 class FSIRingProblem : public Problem
141 {
142  /// There are very few element types that will work for this problem.
143  /// Rather than passing the element type as a template parameter to the
144  /// problem, we choose instead to use a typedef to specify the
145  /// particular element fluid used.
146  typedef AlgebraicElement<RefineableQCrouzeixRaviartElement<2> > FLUID_ELEMENT;
147 
148  /// Typedef to specify the solid element used
149  typedef FSIHermiteBeamElement SOLID_ELEMENT;
150 
151 public:
152 
153  /// Constructor: Number of elements in wall mesh, amplitude of the
154  /// initial wall deformation, amplitude of pcos perturbation and its duration.
155  FSIRingProblem(const unsigned &nelement_wall,
156  const double& eps_ampl, const double& pcos_initial,
157  const double& pcos_duration);
158 
159  /// Update after solve (empty)
161 
162  /// Update before solve (empty)
164 
165  /// Update the problem specs before checking Newton
166  /// convergence
168  {
169  // Update the fluid mesh -- auxiliary update function for algebraic
170  // nodes automatically updates no slip condition.
171  Fluid_mesh_pt->node_update();
172  }
173 
174  /// Update the problem specs after adaptation:
176  {
177  // The functions used to update the no slip boundary conditions
178  // must be set on any new nodes that have been created during the
179  // mesh adaptation process.
180  // There is no mechanism by which auxiliary update functions
181  // are copied to newly created nodes.
182  // (because, unlike boundary conditions, they don't occur exclusively
183  // at boundaries)
184 
185  // The no-slip boundary is boundary 1 of the mesh
186  // Loop over the nodes on this boundary and reset the auxilliary
187  // node update function
188  unsigned n_node = Fluid_mesh_pt->nboundary_node(1);
189  for (unsigned n=0;n<n_node;n++)
190  {
191  Fluid_mesh_pt->boundary_node_pt(1,n)->set_auxiliary_node_update_fct_pt(
192  FSI_functions::apply_no_slip_on_moving_wall);
193  }
194 
195  // (Re-)setup fsi: Work out which fluid dofs affect wall elements
196  // the correspondance between wall dofs and fluid elements is handled
197  // during the remeshing, but the "reverse" association must be done
198  // separately.
199  // We need to set up the interaction every time because the fluid element
200  // adjacent to a given solid element's integration point may have changed
201  // We pass the boundary between the fluid and solid meshes and pointers
202  // to the meshes. The interaction boundary is boundary 1 of the 2D
203  // fluid mesh.
204  FSI_functions::setup_fluid_load_info_for_solid_elements<FLUID_ELEMENT,2>
205  (this,1,Fluid_mesh_pt,Wall_mesh_pt);
206  }
207 
208  /// Doc solution: Pass number of timestep, i (we append to tracefile
209  /// after every timestep but do a full doc only at certain intervals),
210  /// DocInfo object and tracefile
211  void doc_solution(const unsigned& i, DocInfo& doc_info, ofstream& trace_file);
212 
213  /// Do dynamic run
214  void dynamic_run();
215 
216 private:
217 
218  /// Setup initial condition for both domains
219  void set_initial_condition();
220 
221  /// Setup initial condition for wall
222  void set_wall_initial_condition();
223 
224  /// Setup initial condition for fluid
225  void set_fluid_initial_condition();
226 
227  /// Element used for documenting displacement
229 
230  /// Pointer to wall mesh
231  OneDLagrangianMesh<SOLID_ELEMENT> *Wall_mesh_pt;
232 
233  /// Pointer to fluid mesh
234  AlgebraicRefineableQuarterCircleSectorMesh<FLUID_ELEMENT> *Fluid_mesh_pt;
235 
236  /// Pointer to geometric object that represents the undeformed wall shape
237  GeomObject* Undef_geom_pt;
238 
239  /// Pointer to wall timestepper
240  Newmark<2>* Wall_time_stepper_pt;
241 
242  /// Pointer to fluid timestepper
244 
245  /// Pointer to node on coarsest mesh on which velocity is traced
247 
248  /// Amplitude of initial deformation
249  double Eps_ampl;
250 
251  /// Initial pcos
252  double Pcos_initial;
253 
254  /// Duration of initial pcos
256 
257 };
258 
259 
260 //===============================================================
261 /// Setup initial condition: When we're done here, all variables
262 /// represent the state at the initial time.
263 //===============================================================
265 {
266 
267  cout << "Setting wall ic" << std::endl;
268  set_wall_initial_condition();
269 
270  cout << "Setting fluid ic" << std::endl;
271  set_fluid_initial_condition();
272 
273 }
274 
275 
276 //===============================================================
277 /// Setup initial condition for fluid: Impulsive start
278 //===============================================================
280 {
281 
282  // Update fluid domain: Careful!!! This also applies the no slip conditions
283  // on all nodes on the wall! Since the wall might have moved since
284  // we created the mesh; we're therefore imposing a nonzero
285  // velocity on these nodes. Must wipe this afterwards (done
286  // by setting *all* velocities to zero) otherwise we get
287  // an impulsive start from a very bizarre initial velocity
288  // field! [Yes, it took me a while to figure this out...]
289  Fluid_mesh_pt->node_update();
290 
291  // Assign initial values for the velocities;
292  // pressures don't carry a time history and can be left alone.
293 
294  //Find number of nodes in fluid mesh
295  unsigned n_node = Fluid_mesh_pt->nnode();
296 
297  // Loop over the nodes to set initial guess everywhere
298  for(unsigned n=0;n<n_node;n++)
299  {
300  // Loop over velocity directions: Impulsive initial start from
301  // zero velocity!
302  for(unsigned i=0;i<2;i++)
303  {
304  Fluid_mesh_pt->node_pt(n)->set_value(i,0.0);
305  }
306  }
307 
308  // Do an impulsive start with the assigned velocity field
309  Fluid_mesh_pt->assign_initial_values_impulsive();
310 
311 }
312 
313 
314 //===============================================================
315 /// Setup initial condition: Impulsive start either from
316 /// deformed or undeformed wall shape.
317 //===============================================================
319 {
320 
321  // Geometric object that specifies the initial conditions:
322  // A ring that is bucked in a 2-lobed mode
323  GeomObject* ic_geom_object_pt=
324  new PseudoBucklingRing(Eps_ampl,Global_Physical_Variables::H,2,2,
325  Wall_time_stepper_pt);
326 
327  // Assign period of oscillation of the geometric object
328  static_cast<PseudoBucklingRing*>(ic_geom_object_pt)->set_T(1.0);
329 
330  //Set initial time (to deform wall into max. amplitude)
331  double time=0.25;
332 
333  // Assign initial radius of the object
334  static_cast<PseudoBucklingRing*>(ic_geom_object_pt)->set_R_0(1.00);
335 
336  // Setup object that specifies the initial conditions:
337  SolidInitialCondition* IC_pt = new SolidInitialCondition(ic_geom_object_pt);
338 
339  // Assign values of positional data of all elements on wall mesh
340  // so that the wall deforms into the shape specified by IC object.
341  SolidMesh::Solid_IC_problem.set_static_initial_condition(
342  this,Wall_mesh_pt,IC_pt,time);
343 
344 }
345 
346 
347 //========================================================================
348 /// Document solution: Pass number of timestep, i; we append to trace file
349 /// at every timestep and do a full doc only after a certain number
350 /// of steps.
351 //========================================================================
352 void FSIRingProblem::doc_solution(const unsigned& i,
353  DocInfo& doc_info, ofstream& trace_file)
354 {
355 
356  // Full doc every nskip steps
357  unsigned nskip=1; // ADJUST
358 
359  // If we at an integer multiple of nskip, full documentation.
360  if (i%nskip==0)
361  {
362  doc_info.enable_doc();
363  cout << "Full doc step " << doc_info.number()
364  << " for time " << time_stepper_pt()->time_pt()->time() << std::endl;
365  }
366  //Otherwise, just output the trace file
367  else
368  {
369  doc_info.disable_doc();
370  cout << "Only trace for time "
371  << time_stepper_pt()->time_pt()->time() << std::endl;
372  }
373 
374 
375  // If we are at a full documentation step, output the fluid solution
376  if (doc_info.is_doc_enabled())
377  {
378  //Variables used in the output file.
379  ofstream some_file; char filename[100];
380  //Construct the output filename from the doc_info number and the
381  //output directory
382  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
383  doc_info.number());
384  //Open the output file
385  some_file.open(filename);
386  /// Output the solution using 5x5 plot points
387  Fluid_mesh_pt->output(some_file,5);
388  //Close the output file
389  some_file.close();
390  }
391 
392  //Temporary vector to give the local coordinate at which to document
393  //the wall displacment
394  Vector<double> s(1,1.0);
395  // Write to the trace file:
396  trace_file << time_pt()->time()
397  //Document the displacement at the end of the the chosen element
398  << " " << Doc_displacement_elem_pt->interpolated_x(s,1)
399  << " " << Veloc_trace_node_pt->x(0)
400  << " " << Veloc_trace_node_pt->x(1)
401  << " " << Veloc_trace_node_pt->value(0)
402  << " " << Veloc_trace_node_pt->value(1)
403  << " " << Fluid_mesh_pt->nelement()
404  << " " << ndof()
405  << " " << Fluid_mesh_pt->nrefinement_overruled()
406  << " " << Fluid_mesh_pt->max_error()
407  << " " << Fluid_mesh_pt->min_error()
408  << " " << Fluid_mesh_pt->max_permitted_error()
409  << " " << Fluid_mesh_pt->min_permitted_error()
410  << " " << Fluid_mesh_pt->max_keep_unrefined();
411 
412  // Output the number of the corresponding full documentation
413  // file number (or -1 if no full doc was made)
414  if (doc_info.is_doc_enabled())
415  {trace_file << " " <<doc_info.number() << " ";}
416  else {trace_file << " " <<-1 << " ";}
417 
418  //End the trace file
419  trace_file << std::endl;
420 
421  // Increment counter for full doc
422  if (doc_info.is_doc_enabled()) {doc_info.number()++;}
423 }
424 
425 //======================================================================
426 /// Constructor for FSI ring problem. Pass number of wall elements
427 /// and length of wall (in Lagrangian coordinates) amplitude of
428 /// initial deformation, pcos perturbation and duration.
429 //======================================================================
431  const double& eps_ampl, const double& pcos_initial,
432  const double& pcos_duration) :
433  Eps_ampl(eps_ampl), Pcos_initial(pcos_initial),
434  Pcos_duration(pcos_duration)
435 {
436  //-----------------------------------------------------------
437  // Create timesteppers
438  //-----------------------------------------------------------
439 
440  // Allocate the wall timestepper and add it to the problem's vector
441  // of timesteppers
442  Wall_time_stepper_pt = new Newmark<2>;
443  add_time_stepper_pt(Wall_time_stepper_pt);
444 
445  // Allocate the fluid timestepper and add it to the problem's Vector
446  // of timesteppers
447  Fluid_time_stepper_pt = new BDF<2>;
448  add_time_stepper_pt(Fluid_time_stepper_pt);
449 
450  //----------------------------------------------------------
451  // Create the wall mesh
452  //----------------------------------------------------------
453 
454  // Undeformed wall is an elliptical ring
455  Undef_geom_pt = new Ellipse(1.0,1.0);
456 
457  //Length of wall in Lagrangian coordinates
458  double L = 2.0*atan(1.0);
459 
460  //Now create the (Lagrangian!) mesh
461  Wall_mesh_pt = new
462  OneDLagrangianMesh<SOLID_ELEMENT>(N,L,Undef_geom_pt,Wall_time_stepper_pt);
463 
464  //----------------------------------------------------------
465  // Set the boundary conditions for wall mesh (problem)
466  //----------------------------------------------------------
467 
468  // Bottom boundary: (Boundary 0)
469  // No vertical displacement
470  Wall_mesh_pt->boundary_node_pt(0,0)->pin_position(1);
471  // Zero slope: Pin type 1 dof for displacement direction 0
472  Wall_mesh_pt->boundary_node_pt(0,0)->pin_position(1,0);
473 
474  // Top boundary: (Boundary 1)
475  // No horizontal displacement
476  Wall_mesh_pt->boundary_node_pt(1,0)->pin_position(0);
477  // Zero slope: Pin type 1 dof for displacement direction 1
478  Wall_mesh_pt->boundary_node_pt(1,0)->pin_position(1,1);
479 
480 
481  //-----------------------------------------------------------
482  // Create the fluid mesh:
483  //-----------------------------------------------------------
484 
485  // Fluid mesh is suspended from wall between the following Lagrangian
486  // coordinates:
487  double xi_lo=0.0;
488  double xi_hi=L;
489 
490  // Fractional position of dividing line for two outer blocks in mesh
491  double fract_mid=0.5;
492 
493  //Create a geometric object that represents the wall geometry from the
494  //wall mesh (one Lagrangian, two Eulerian coordinates).
495  MeshAsGeomObject *wall_mesh_as_geometric_object_pt
496  = new MeshAsGeomObject(Wall_mesh_pt);
497 
498  // Build fluid mesh using the wall mesh as a geometric object
499  Fluid_mesh_pt = new AlgebraicRefineableQuarterCircleSectorMesh<FLUID_ELEMENT >
500  (wall_mesh_as_geometric_object_pt,
501  xi_lo,fract_mid,xi_hi,Fluid_time_stepper_pt);
502 
503  // Set the error estimator
504  Z2ErrorEstimator* error_estimator_pt=new Z2ErrorEstimator;
505  Fluid_mesh_pt->spatial_error_estimator_pt()=error_estimator_pt;
506 
507  // Extract pointer to node at center of mesh
508  unsigned nnode=Fluid_mesh_pt->finite_element_pt(0)->nnode();
509  Veloc_trace_node_pt=Fluid_mesh_pt->finite_element_pt(0)->node_pt(nnode-1);
510 
511  //-------------------------------------------------------
512  // Set the fluid boundary conditions
513  //-------------------------------------------------------
514 
515  // Bottom boundary (boundary 0):
516  {
517  unsigned n_node = Fluid_mesh_pt->nboundary_node(0);
518  for (unsigned n=0;n<n_node;n++)
519  {
520  // Pin vertical velocity
521  Fluid_mesh_pt->boundary_node_pt(0,n)->pin(1);
522  }
523  }
524 
525  // Ring boundary (boundary 1):
526  // No slip; this also implies that the velocity needs
527  // to be updated in response to wall motion
528  {
529  unsigned n_node = Fluid_mesh_pt->nboundary_node(1);
530  for (unsigned n=0;n<n_node;n++)
531  {
532  // Which node are we dealing with?
533  Node* node_pt=Fluid_mesh_pt->boundary_node_pt(1,n);
534 
535  // Set auxiliary update function pointer
536  node_pt->set_auxiliary_node_update_fct_pt(
537  FSI_functions::apply_no_slip_on_moving_wall);
538 
539  // Pin both velocities
540  for(unsigned i=0;i<2;i++) {node_pt->pin(i);}
541  }
542  }
543 
544  // Left boundary (boundary 2):
545  {
546  unsigned n_node = Fluid_mesh_pt->nboundary_node(2);
547  for (unsigned n=0;n<n_node;n++)
548  {
549  // Pin horizontal velocity
550  Fluid_mesh_pt->boundary_node_pt(2,n)->pin(0);
551  }
552  }
553 
554 
555  //--------------------------------------------------------
556  // Add the submeshes and build global mesh
557  // -------------------------------------------------------
558 
559  // Wall mesh
560  add_sub_mesh(Wall_mesh_pt);
561 
562  //Fluid mesh
563  add_sub_mesh(Fluid_mesh_pt);
564 
565  // Combine all submeshes into a single Mesh
566  build_global_mesh();
567 
568 
569  //----------------------------------------------------------
570  // Finish problem setup
571  // ---------------------------------------------------------
572 
573  //Find number of elements in fluid mesh
574  unsigned n_element = Fluid_mesh_pt->nelement();
575 
576  // Loop over the fluid elements to set up element-specific
577  // things that cannot be handled by constructor
578  for(unsigned e=0;e<n_element;e++)
579  {
580  // Upcast from FiniteElement to the present element
581  FLUID_ELEMENT *el_pt
582  = dynamic_cast<FLUID_ELEMENT*>(Fluid_mesh_pt->element_pt(e));
583 
584  //Set the Reynolds number, etc
585  el_pt->re_pt() = &Global_Physical_Variables::Re;
586  el_pt->re_st_pt() = &Global_Physical_Variables::ReSt;
587 
588 
589  el_pt->evaluate_shape_derivs_by_direct_fd();
590 
591 // el_pt->evaluate_shape_derivs_by_chain_rule();
592 // el_pt->enable_always_evaluate_dresidual_dnodal_coordinates_by_fd();
593 
594 // if (e==0)
595 // {
596 // el_pt->disable_always_evaluate_dresidual_dnodal_coordinates_by_fd();
597 // }
598 // else
599 // {
600 // el_pt->enable_always_evaluate_dresidual_dnodal_coordinates_by_fd();
601 // }
602 
603  //el_pt->evaluate_shape_derivs_by_direct_fd();
604 
605  }
606 
607 
608  //Loop over the solid elements to set physical parameters etc.
609  unsigned n_wall_element = Wall_mesh_pt->nelement();
610  for(unsigned e=0;e<n_wall_element;e++)
611  {
612  //Cast to proper element type
613  SOLID_ELEMENT *el_pt = dynamic_cast<SOLID_ELEMENT*>(
614  Wall_mesh_pt->element_pt(e));
615 
616  //Set physical parameters for each element:
617  el_pt->h_pt() = &Global_Physical_Variables::H;
618  el_pt->lambda_sq_pt() = &Global_Physical_Variables::Lambda_sq;
619 
620  //Function that specifies the external load Vector
621  el_pt->load_vector_fct_pt() = &Global_Physical_Variables::pcos_load;
622 
623  // Function that specifies the load ratios
624  el_pt->q_pt() = &Global_Physical_Variables::Q;
625 
626  //Assign the undeformed beam shape
627  el_pt->undeformed_beam_pt() = Undef_geom_pt;
628  }
629 
630  // Establish control displacment: (even if no displacement control is applied
631  // we still want to doc the displacement at the same point)
632 
633  // Choose element: (This is the last one)
634  Doc_displacement_elem_pt = dynamic_cast<SOLID_ELEMENT*>(
635  Wall_mesh_pt->element_pt(n_wall_element-1));
636 
637  // Setup fsi: Work out which fluid dofs affect the wall elements
638  // the correspondance between wall dofs and fluid elements is handled
639  // during the remeshing, but the "reverse" association must be done
640  // separately.
641  // We pass the boundary between the fluid and solid meshes and pointers
642  // to the meshes. The interaction boundary is boundary 1 of the
643  // 2D fluid mesh.
644  FSI_functions::setup_fluid_load_info_for_solid_elements<FLUID_ELEMENT,2>
645  (this,1,Fluid_mesh_pt,Wall_mesh_pt);
646 
647  // Do equation numbering
648  cout << "# of dofs " << assign_eqn_numbers() << std::endl;
649 
650 }
651 
652 
653 //=========================================================================
654 /// Solver loop to perform unsteady run
655 //=========================================================================
657 {
658  // Setup documentation
659  //---------------------------------------------------------------
660 
661  /// Label for output
662  DocInfo doc_info;
663 
664  // Output directory
665  doc_info.set_directory("RESLT");
666 
667  // Step number
668  doc_info.number()=0;
669 
670  //Open a trace file
671  ofstream trace_file("RESLT/trace_ring.dat");
672 
673  // Write header for trace file
674  trace_file << "VARIABLES=\"time\",\"V_c_t_r_l\"";
675  trace_file << ",\"x<sub>1</sub><sup>(track)</sup>\"";
676  trace_file << ",\"x<sub>2</sub><sup>(track)</sup>\"";
677  trace_file << ",\"u<sub>1</sub><sup>(track)</sup>\"";
678  trace_file << ",\"u<sub>2</sub><sup>(track)</sup>\"";
679  trace_file << ",\"N<sub>element</sub>\"";
680  trace_file << ",\"N<sub>dof</sub>\"";
681  trace_file << ",\"# of under-refined elements\"";
682  trace_file << ",\"max. error\"";
683  trace_file << ",\"min. error\"";
684  trace_file << ",\"max. permitted error\"";
685  trace_file << ",\"min. permitted error\"";
686  trace_file << ",\"max. permitted # of unrefined elements\"";
687  trace_file << ",\"doc number\"";
688  trace_file << std::endl;
689 
690 
691  // Initialise timestepping
692  // -------------------------------------------------------------
693 
694  // Number of steps
695  unsigned nstep=300;
696 
697  // Nontrivial command line input: validation: only do three steps
698  if (CommandLineArgs::Argc>1)
699  {
700  nstep=1;
701  cout << "Only doing nstep steps for validation: " << nstep << std::endl;
702  }
703 
704  // Set initial timestep
705  double dt=0.004;
706 
707  // Set initial value for dt -- also assigns weights to the timesteppers
708  initialise_dt(dt);
709 
710  // Set physical parameters
711  //---------------------------------------------------------
712  using namespace Global_Physical_Variables;
713 
714  // Set Womersley number
715  Alpha_sq=100.0; // 50.0; // ADJUST
716 
717  // Set density ratio
718  Density_ratio=10.0; // 0.0; ADJUST
719 
720  // Wall thickness
721  H=1.0/20.0;
722 
723  // Set external pressure
724  Pext=0.0;
725 
726  /// Perturbation pressure
728 
729  // Assign/doc corresponding computational parameters
730  set_params();
731 
732 
733  // Refine uniformly and assign initial conditions
734  //--------------------------------------------------------------
735 
736  // Refine the problem uniformly
737  refine_uniformly();
738  refine_uniformly();
739 
740  // This sets up the solution at the initial time
742 
743  // Set targets for spatial adptivity
744  //---------------------------------------------------------------
745 
746  // Max. and min. error for adaptive refinement/unrefinement
747  Fluid_mesh_pt->max_permitted_error()=1.0e-2;
748  Fluid_mesh_pt->min_permitted_error()=1.0e-3;
749 
750  // Don't allow refinement to drop under given level
751  Fluid_mesh_pt->min_refinement_level()=2;
752 
753  // Don't allow refinement beyond given level
754  Fluid_mesh_pt->max_refinement_level()=6;
755 
756  // Don't bother adapting the mesh if no refinement is required
757  // and if less than ... elements are to be merged.
758  Fluid_mesh_pt->max_keep_unrefined()=20;
759 
760  // Doc refinement targets
761  Fluid_mesh_pt->doc_adaptivity_targets(cout);
762 
763 
764  // Do the timestepping
765  //----------------------------------------------------------------
766 
767  // Reset initial conditions after refinement for first step only
768  bool first=true;
769 
770  //Output initial data
771  doc_solution(0,doc_info,trace_file);
772 
773 
774 // {
775 // unsigned nel=Fluid_mesh_pt->nelement();
776 // for (unsigned e=0;e<nel;e++)
777 // {
778 // std::cout << "\n\nEl: " << e << std::endl << std::endl;
779 // FiniteElement* el_pt=Fluid_mesh_pt->finite_element_pt(e);
780 // unsigned n_dof=el_pt->ndof();
781 // Vector<double> residuals(n_dof);
782 // DenseDoubleMatrix jacobian(n_dof,n_dof);
783 // el_pt->get_jacobian(residuals,jacobian);
784 // }
785 // exit(0);
786 // }
787 
788  // Time integration loop
789  for(unsigned i=1;i<=nstep;i++)
790  {
791  // Switch doc off during solve
792  doc_info.disable_doc();
793 
794  // Solve
795  unsigned max_adapt=1;
796  unsteady_newton_solve(dt,max_adapt,first);
797 
798  // Now we've done the first step
799  first=false;
800 
801  // Doc solution
802  doc_solution(i,doc_info,trace_file);
803 
804  /// Switch off perturbation pressure
805  if (time_pt()->time()>Pcos_duration) {Pcos=0.0;}
806 
807  }
808 
809 }
810 
811 
812 //=====================================================================
813 /// Driver for fsi ring test problem
814 //=====================================================================
815 int main(int argc, char* argv[])
816 {
817 
818  // Store command line arguments
819  CommandLineArgs::setup(argc,argv);
820 
821  // Number of elements
822  unsigned nelem = 13;
823 
824  /// Perturbation pressure
825  double pcos_initial=1.0e-6; // ADJUST
826 
827  /// Duration of initial pcos perturbation
828  double pcos_duration=0.3; // ADJUST
829 
830  /// Amplitude of initial deformation
831  double eps_ampl=0.0; // ADJUST
832 
833  //Set up the problem
834  FSIRingProblem problem(nelem,eps_ampl,pcos_initial,pcos_duration);
835 
836  // Do parameter study
837  problem.dynamic_run();
838 
839 }
840 
841 
842 
843 
844 
845 
FSI Ring problem: a fluid-structure interaction problem in which a viscous fluid bounded by an initia...
double Pcos_initial
Initial pcos.
OneDLagrangianMesh< SOLID_ELEMENT > * Wall_mesh_pt
Pointer to wall mesh.
GeomObject * Undef_geom_pt
Pointer to geometric object that represents the undeformed wall shape.
AlgebraicElement< RefineableQCrouzeixRaviartElement< 2 > > FLUID_ELEMENT
There are very few element types that will work for this problem. Rather than passing the element typ...
Node * Veloc_trace_node_pt
Pointer to node on coarsest mesh on which velocity is traced.
void set_initial_condition()
Setup initial condition for both domains.
AlgebraicRefineableQuarterCircleSectorMesh< FLUID_ELEMENT > * Fluid_mesh_pt
Pointer to fluid mesh.
void set_wall_initial_condition()
Setup initial condition for wall.
void actions_after_newton_solve()
Update after solve (empty)
Newmark< 2 > * Wall_time_stepper_pt
Pointer to wall timestepper.
void doc_solution(const unsigned &i, DocInfo &doc_info, ofstream &trace_file)
Doc solution: Pass number of timestep, i (we append to tracefile after every timestep but do a full d...
BDF< 2 > * Fluid_time_stepper_pt
Pointer to fluid timestepper.
double Eps_ampl
Amplitude of initial deformation.
SOLID_ELEMENT * Doc_displacement_elem_pt
Element used for documenting displacement.
FSIRingProblem(const unsigned &nelement_wall, const double &eps_ampl, const double &pcos_initial, const double &pcos_duration)
Constructor: Number of elements in wall mesh, amplitude of the initial wall deformation,...
double Pcos_duration
Duration of initial pcos.
FSIHermiteBeamElement SOLID_ELEMENT
Typedef to specify the solid element used.
void actions_before_newton_solve()
Update before solve (empty)
void dynamic_run()
Do dynamic run.
void set_fluid_initial_condition()
Setup initial condition for fluid.
void actions_after_adapt()
Update the problem specs after adaptation:
void actions_before_newton_convergence_check()
Update the problem specs before checking Newton convergence.
int main(int argc, char *argv[])
Driver for fsi ring test problem.
Namespace for physical parameters.
Definition: fsi_osc_ring.cc:46
double Pext
External Pressure.
Definition: fsi_osc_ring.cc:58
double Alpha_sq
Square of Womersly number (a frequency parameter)
Definition: fsi_osc_ring.cc:52
double ReSt
Reynolds x Strouhal number.
Definition: fsi_osc_ring.cc:77
double Nu
Poisson ratio.
Definition: fsi_osc_ring.cc:61
void set_params()
Set the parameters that are used in the code as a function of the Womersley number,...
Definition: fsi_osc_ring.cc:87
double Q
Stress ratio.
Definition: fsi_osc_ring.cc:83
double Lambda_sq
Timescale ratio (non-dimensation density)
Definition: fsi_osc_ring.cc:80
void pcos_load(const Vector< double > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
Non-FSI load function, a constant external pressure plus a (small) sinusoidal perturbation of wavenum...
double Pcos
Perturbation pressure.
Definition: fsi_osc_ring.cc:67
double Re
Reynolds number.
Definition: fsi_osc_ring.cc:74
double Density_ratio
Density ratio of the solid and the fluid.
Definition: fsi_osc_ring.cc:55
double H
Nondimensional thickness of the beam.
Definition: fsi_osc_ring.cc:64