turek_flag_non_fsi.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 //------------------------------------------------------------------------
27 // Flow past cylinder with flag -- first draft developed by
28 // Stefan Kollmannsberger and students Iason Papaioannou and
29 // Orkun Oezkan Doenmez.
30 //
31 // http://www.inf.bauwesen.tu-muenchen.de/~kollmannsberger/
32 //
33 // Completed by Floraine Cordier.
34 //------------------------------------------------------------------------
35 
36 // Generic oomph-lib stuff
37 #include "generic.h"
38 
39 // Navier Stokes
40 #include "navier_stokes.h"
41 
42 // The mesh
43 #include "meshes/cylinder_with_flag_mesh.h"
44 
45 using namespace std;
46 
47 using namespace oomph;
48 
49 //#define USE_MACRO_ELEMENTS
50 
51 
52 //====start_of_global_parameters==========================================
53 /// Global parameters
54 //========================================================================
56 {
57  /// Reynolds number
58  double Re=100.0;
59 }
60 
61 
62 
63 //====start_of_flag_definition===========================================
64 /// Namespace for definition of flag boundaries
65 //=======================================================================
66 namespace Flag_definition
67 {
68 
69  /// Period of prescribed flag oscillation
70  double Period=10.0;
71 
72  /// Height of flag
73  double H=0.2;
74 
75  /// Length of flag
76  double L=3.5;
77 
78  /// x position of centre of cylinder
79  double Centre_x=2.0;
80 
81  /// y position of centre of cylinder
82  double Centre_y=2.0;
83 
84  /// Radius of cylinder
85  double Radius=0.5;
86 
87  /// Amplitude of tip deflection
88  double Amplitude=0.33;
89 
90  /// Pointer to the global time object
91  Time* Time_pt=0;
92 
93 
94 
95  /// Time-dependent vector to upper tip of the "flag"
96  Vector<double> upper_tip(const double& t)
97  {
98  double tmp_ampl=Amplitude;
99  if (t<=0.0) tmp_ampl=0.0;
100  Vector<double> uppertip(2);
101  uppertip[0]= Centre_x+Radius*sqrt(1.0-H*H/(4.0*Radius*Radius))+L;
102  uppertip[1]= Centre_y+0.5*H-
103  tmp_ampl*sin(2.0*MathematicalConstants::Pi*t/Period);
104  return uppertip;
105  }
106 
107  /// Time-dependent vector to bottom tip of the "flag"
108  Vector<double> lower_tip(const double& t)
109  {
110  double tmp_ampl=Amplitude;
111  if (t<=0.0) tmp_ampl=0.0;
112  Vector<double> lowertip(2);
113  lowertip[0]= Centre_x+Radius*sqrt(1.0-H*H/(4.0*Radius*Radius))+L;
114  lowertip[1]= Centre_y-0.5*H-
115  tmp_ampl*sin(2.0*MathematicalConstants::Pi*t/Period);
116  return lowertip;
117  } // end of bottom tip
118 
119 
120 
121 
122 
123 //-----start_of_top_of_flag--------------------------------------
124 /// GeomObject that defines the upper boundary of the flag
125 //---------------------------------------------------------------
126  class TopOfFlag : public GeomObject
127  {
128 
129  public:
130 
131  /// Constructor: It's a 1D object in 2D space
132  TopOfFlag() : GeomObject(1,2) {}
133 
134  /// Destructor (empty)
136 
137  /// Return the position along the top of the flag (xi[0] varies
138  /// between 0 and Lx)
139  void position(const unsigned& t,const Vector<double> &xi, Vector<double> &r)
140  const
141  {
142  // Compute position of fixed point on the cylinder
143  Vector<double> point_on_circle(2);
144  point_on_circle[0]=Centre_x+Radius*sqrt(1.0-H*H/(4.0*Radius*Radius));
145  point_on_circle[1]=Centre_y+H/2.0;
146 
147  r[0] = point_on_circle[0]+xi[0]/L*(upper_tip(Time_pt->time(t))[0]-
148  point_on_circle[0]);
149 
150  r[1] = point_on_circle[1]+xi[0]/L*(upper_tip(Time_pt->time(t))[1]-
151  point_on_circle[1])+
152  1.0/3.0*sin((r[0]-point_on_circle[0])/
153  (upper_tip(Time_pt->time(t))[0]-
154  point_on_circle[0])*MathematicalConstants::Pi)
155  *sin(2.0* MathematicalConstants::Pi*Time_pt->time(t)/Period);
156 
157  }
158 
159 
160  /// Current position
161  void position(const Vector<double> &xi, Vector<double> &r) const
162  {
163  return position(0,xi,r);
164  }
165 
166  /// Number of geometric Data in GeomObject: None.
167  unsigned ngeom_data() const {return 0;}
168 
169 };
170 
171 
172 
173 //-----start_of_bottom_of_flag-----------------------------------
174 /// GeomObject that defines the lower boundary of the flag
175 //---------------------------------------------------------------
176  class BottomOfFlag : public GeomObject
177  {
178 
179  public:
180 
181  /// Constructor:
182  BottomOfFlag() : GeomObject(1,2) {}
183 
184  /// Destructor (empty)
186 
187 
188  /// Return the position along the bottom of the flag (xi[0] varies
189  /// between 0 and Lx)
190  void position(const unsigned& t,const Vector<double> &xi, Vector<double> &r)
191  const
192  {
193  // Compute position of fixed point on the cylinder
194  Vector<double> point_on_circle(2);
195  point_on_circle[0]=Centre_x+Radius*sqrt(1.0-H*H/(4.0*Radius*Radius));
196  point_on_circle[1]=Centre_y-H/2.0;
197 
198  r[0] = point_on_circle[0]+ xi[0]/L*(lower_tip(Time_pt->time(t))[0]-
199  point_on_circle[0]);
200 
201  r[1] = point_on_circle[1]+ xi[0]/L*(lower_tip(Time_pt->time(t))[1]-
202  point_on_circle[1])+
203  1.0/3.0*sin((r[0]-point_on_circle[0])/
204  (lower_tip(Time_pt->time(t))[0]-
205  point_on_circle[0])*MathematicalConstants::Pi)
206  *sin(2.0*MathematicalConstants::Pi*Time_pt->time(t)/Period);
207  }
208 
209 
210  /// Current position
211  void position(const Vector<double> &xi, Vector<double> &r) const
212  {
213  return position(0,xi,r);
214  }
215 
216  /// Number of geometric Data in GeomObject: None.
217  unsigned ngeom_data() const {return 0;}
218 
219 };
220 
221 
222 
223 
224  //-----start_of_tip_of_flag--------------------------------------
225  /// GeomObject that defines the tip of the flag
226  //---------------------------------------------------------------
227  class TipOfFlag : public GeomObject
228  {
229 
230  public:
231 
232  /// Constructor
233  TipOfFlag() : GeomObject(1,2) {}
234 
235  /// Destructor (empty)
237 
238  /// Return the position
239  /// This object links the tips of the top and bottom by a straight line
240  /// whilst xi[0] goes from -H/2 to H/2.
241  void position(const unsigned& t,const Vector<double> &xi, Vector<double> &r)
242  const
243  {
244  Vector<double> point_top(upper_tip(Time_pt->time(t)));
245  Vector<double> point_bottom(lower_tip(Time_pt->time(t)));
246 
247  r[1]= point_bottom[1]+(xi[0]+H/2.0)/H*(point_top[1]-point_bottom[1]);
248  r[0]= point_bottom[0]+(xi[0]+H/2.0)/H*(point_top[0]-point_bottom[0]);
249  }
250 
251  /// Current position.
252  void position(const Vector<double> &xi, Vector<double> &r) const
253  {
254  return position(0,xi,r);
255  }
256 
257  /// Number of geometric Data in GeomObject: None.
258  unsigned ngeom_data() const {return 0;}
259 
260  };
261 
262 
263  /// Pointer to GeomObject that bounds the upper edge of the flag
265 
266  /// Pointer to GeomObject that bounds the bottom edge of the flag
268 
269  /// Pointer to GeomObject that bounds the tip edge of the flag
271 
272  /// Pointer to GeomObject of type Circle that defines the
273  /// central cylinder.
274  Circle* Cylinder_pt=0;
275 
276  /// Create all GeomObjects needed to define the cylinder and the flag
277  void setup(Time* time_pt)
278  {
279  // Assign pointer to global time object
280  Time_pt=time_pt;
281 
282  // Create GeomObject of type Circle that defines the
283  // central cylinder.
284  Cylinder_pt=new Circle(Centre_x,Centre_y,Radius);
285 
286  /// Create GeomObject that bounds the upper edge of the flag
288 
289  /// Create GeomObject that bounds the bottom edge of the flag
291 
292  /// Create GeomObject that bounds the tip edge of the flag
294 
295  }
296 
297 } // end of namespace that specifies the flag
298 
299 
300 
301 /// /////////////////////////////////////////////////////////////////
302 /// /////////////////////////////////////////////////////////////////
303 /// /////////////////////////////////////////////////////////////////
304 
305 
306 
307 //======= start_of_problem_class=====================================
308 /// Flow around a cylinder with flag
309 //===================================================================
310 template<class ELEMENT>
311 class TurekNonFSIProblem : public Problem
312 {
313 
314 public:
315 
316  /// Constructor: Pass length and height of domain.
317  TurekNonFSIProblem(const double &length,
318  const double &height);
319 
320  /// Update the problem specs after solve (empty)
322 
323  /// Update the problem specs before solve (empty)
325 
326  /// After adaptation: Unpin pressures and pin redudant pressure dofs.
327  void actions_after_adapt();
328 
329  /// Update the velocity boundary condition on the flag
330  void actions_before_implicit_timestep();
331 
332 
333 #ifdef USE_MACRO_ELEMENTS
334 
335  /// Access function for the specific mesh
336  RefineableCylinderWithFlagMesh<ELEMENT>* mesh_pt()
337  {
338  return dynamic_cast<RefineableCylinderWithFlagMesh<ELEMENT>*>
339  (Problem::mesh_pt());
340  }
341 
342 #else
343 
344  /// Access function for the specific mesh
345  RefineableAlgebraicCylinderWithFlagMesh<ELEMENT>* mesh_pt()
346  {
347  return dynamic_cast<RefineableAlgebraicCylinderWithFlagMesh<ELEMENT>*>
348  (Problem::mesh_pt());
349  }
350 
351 #endif
352 
353  /// Doc the solution
354  void doc_solution(DocInfo& doc_info);
355 
356 
357 private:
358 
359  /// Height of channel
360  double Height;
361 
362 };//end_of_problem_class
363 
364 
365 
366 
367 //=====start_of_constructor===============================================
368 /// Constructor: Pass length and height of domain.
369 //========================================================================
370 template<class ELEMENT>
372  const double &length, const double &height) : Height(height)
373 {
374 
375  // Bump up Newton solver parameters to allow for crappy initial guesses
376  Max_residuals=100.0;;
377  Max_newton_iterations=50;
378 
379  // Allocate the timestepper
380  add_time_stepper_pt(new BDF<2>);
381 
382  // Setup flag/cylinder geometry
383  Flag_definition::setup(time_pt());
384 
385 #ifdef USE_MACRO_ELEMENTS
386 
387 
388 
389  // Build mesh (with Domain/MacroElement-based node update)
390  Problem::mesh_pt()=new RefineableCylinderWithFlagMesh<ELEMENT>
395  length, height,
401  time_stepper_pt());
402 
403  std::cout << "Using Domain/MacroElement-based node update" << std::endl;
404 
405 #else
406 
407  // Build mesh (with AlgebraicMesh-based node update)
408  Problem::mesh_pt()=new RefineableAlgebraicCylinderWithFlagMesh<ELEMENT>
413  length, height,
419  time_stepper_pt());
420 
421  std::cout << "Using AlgebraicMesh-based node update" << std::endl;
422 
423 #endif
424 
425  // Refine uniformly
426  mesh_pt()->refine_uniformly();
427  mesh_pt()->refine_uniformly();
428 
429  // Set error estimator
430  Z2ErrorEstimator* error_estimator_pt=new Z2ErrorEstimator;
431  mesh_pt()->spatial_error_estimator_pt()=error_estimator_pt;
432 
433  // Set the boundary conditions for this problem: All nodes are
434  // free by default -- just pin the ones that have Dirichlet conditions
435  // here.
436 
437  //Pin both velocities at all boundaries apart from outflow
438  unsigned num_bound = mesh_pt()->nboundary();
439  for(unsigned ibound=0;ibound<num_bound;ibound++)
440  {
441  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
442  for (unsigned inod=0;inod<num_nod;inod++)
443  {
444  // Parallel, axially traction free outflow at downstream end
445  if (ibound != 1)
446  {
447  mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
448  mesh_pt()->boundary_node_pt(ibound,inod)->pin(1);
449  }
450  else
451  {
452  mesh_pt()->boundary_node_pt(ibound,inod)->pin(1);
453  }
454  }
455  } // done bc
456 
457 
458  // Pass pointer to Reynolds number to elements
459  unsigned nelem=mesh_pt()->nelement();
460  for (unsigned e=0;e<nelem;e++)
461  {
462  // Upcast from GeneralisedElement to the present element
463  ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
464 
465  //Set the Reynolds number
466  el_pt->re_pt() = &Global_Parameters::Re;
467 
468  //Set the Womersley number (assuming St=1)
469  el_pt->re_st_pt() = &Global_Parameters::Re;
470  }
471 
472 
473  // Setup Poiseuille flow along boundary 3
474  unsigned ibound=3;
475  unsigned num_nod= mesh_pt()->nboundary_node(ibound);
476  for (unsigned inod=0;inod<num_nod;inod++)
477  {
478  double ycoord = mesh_pt()->boundary_node_pt(ibound,inod)->x(1);
479 
480  // Set Poiseuille velocity
481  double uy = 6.0*ycoord/Height*(1.0-ycoord/Height);
482 
483  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(0,uy);
484  mesh_pt()->boundary_node_pt(ibound,inod)->set_value(1,0.0);
485  }
486 
487 
488  // Pin redudant pressure dofs
489  RefineableNavierStokesEquations<2>::
490  pin_redundant_nodal_pressures(mesh_pt()->element_pt());
491 
492  // Assign equations numbers
493  cout <<"Number of equations: " << assign_eqn_numbers() << endl;
494 
495 } // end of constructor
496 
497 
498 
499 //=====actions_after_adapt================================================
500 /// Actions after adapt
501 //========================================================================
502 template <class ELEMENT>
504 {
505  // Unpin all pressure dofs
506  RefineableNavierStokesEquations<2>::
507  unpin_all_pressure_dofs(mesh_pt()->element_pt());
508 
509  // Pin redundant pressure dofs
510  RefineableNavierStokesEquations<2>::
511  pin_redundant_nodal_pressures(mesh_pt()->element_pt());
512 
513 } // end_of_actions_after_adapt
514 
515 
516 //==== start_of_actions_before_implicit_timestep==========================
517 /// Actions before implicit timestep: Update velocity boundary conditions
518 //========================================================================
519 template <class ELEMENT>
521 {
522 
523  // Update the domain shape
524  mesh_pt()->node_update();
525 
526  // Moving leaflet: No slip; this implies that the velocity needs
527  // to be updated in response to flag motion
528  for( unsigned ibound=5;ibound<8;ibound++)
529  {
530  unsigned num_nod=mesh_pt()->nboundary_node(ibound);
531  for (unsigned inod=0;inod<num_nod;inod++)
532  {
533  // Which node are we dealing with?
534  Node* node_pt=mesh_pt()->boundary_node_pt(ibound,inod);
535 
536  // Apply no slip
537  FSI_functions::apply_no_slip_on_moving_wall(node_pt);
538  }
539  }
540 } //end_of_actions_before_implicit_timestep
541 
542 
543 //==start_of_doc_solution=================================================
544 /// Doc the solution
545 //========================================================================
546 template<class ELEMENT>
548 {
549 
550  ofstream some_file;
551  char filename[100];
552 
553  // Number of plot points
554  unsigned npts;
555  npts=5;
556 
557  // Output solution
558  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
559  doc_info.number());
560  some_file.open(filename);
561  mesh_pt()->output(some_file,npts);
562  some_file.close();
563 
564 } // end_of_doc_solution
565 
566 
567 /// ////////////////////////////////////////////////////////////////////
568 /// ////////////////////////////////////////////////////////////////////
569 /// ////////////////////////////////////////////////////////////////////
570 
571 
572 //======start_of_main==================================================
573 /// Driver code -- pass a command line argument if you want to run
574 /// the code in validation mode where it only performs a few steps
575 //=====================================================================
576 int main(int argc, char* argv[])
577 {
578 
579  // Store command line arguments
580  CommandLineArgs::setup(argc,argv);
581 
582  // Set up doc info
583  DocInfo doc_info;
584  doc_info.set_directory("RESLT");
585  doc_info.number()=0;
586 
587  // Length and height of domain
588  double length=25.0;
589  double height=4.1;
590 
591 #ifdef USE_MACRO_ELEMENTS
592 
593  // Create Problem with Domain/MacroElement-based node update
595  problem(length,height);
596 
597 #else
598 
599  // Create Problem with AlgebraicMesh-based node update
601  <AlgebraicElement<RefineableQCrouzeixRaviartElement<2> > >
602  problem(length, height);
603 
604 #endif
605 
606 
607 // Number of timesteps per period
608  unsigned nsteps_per_period=40;
609 
610  // Number of periods
611  unsigned nperiod=3;
612 
613  // Number of timesteps (reduced for validation)
614  unsigned nstep=nsteps_per_period*nperiod;
615  if (CommandLineArgs::Argc>1)
616  {
617  nstep=2;
618  // Also reduce the Reynolds number to reduce the mesh refinement
620  }
621 
622  //Timestep:
623  double dt=Flag_definition::Period/double(nsteps_per_period);
624 
625  /// Initialise timestep
626  problem.initialise_dt(dt);
627 
628  // Solve adaptively with up to max_adapt rounds of refinement (fewer if
629  // run during self-test)
630  unsigned max_adapt=3;
631  if (CommandLineArgs::Argc>1)
632  {
633  max_adapt=1;
634  }
635 
636  /// Output intial guess for steady Newton solve
637  problem.doc_solution(doc_info);
638  doc_info.number()++;
639 
640  // Do steady solve first -- this also sets the history values
641  // to those corresponding to an impulsive start from the
642  // steady solution
643  problem.steady_newton_solve(max_adapt);
644 
645  /// Output steady solution = initial condition for subsequent unsteady solve
646  problem.doc_solution(doc_info);
647  doc_info.number()++;
648 
649  /// Reduce the max number of adaptations for time-dependent simulation
650  max_adapt=1;
651 
652  // We don't want to re-assign the initial condition after the mesh
653  // adaptation
654  bool first=false;
655 
656 // Timestepping loop
657  for (unsigned istep=0;istep<nstep;istep++)
658  {
659  // Solve the problem
660  problem.unsteady_newton_solve(dt,max_adapt,first);
661 
662  // Output the solution
663  problem.doc_solution(doc_info);
664 
665  // Step number
666  doc_info.number()++;
667  }
668 
669 }//end of main
670 
671 
672 
GeomObject that defines the lower boundary of the flag.
void position(const unsigned &t, const Vector< double > &xi, Vector< double > &r) const
Return the position along the bottom of the flag (xi[0] varies between 0 and Lx)
unsigned ngeom_data() const
Number of geometric Data in GeomObject: None.
void position(const Vector< double > &xi, Vector< double > &r) const
Current position.
GeomObject that defines the tip of the flag.
unsigned ngeom_data() const
Number of geometric Data in GeomObject: None.
~TipOfFlag()
Destructor (empty)
void position(const Vector< double > &xi, Vector< double > &r) const
Current position.
void position(const unsigned &t, const Vector< double > &xi, Vector< double > &r) const
Return the position This object links the tips of the top and bottom by a straight line whilst xi[0] ...
GeomObject that defines the upper boundary of the flag.
void position(const Vector< double > &xi, Vector< double > &r) const
Current position.
TopOfFlag()
Constructor: It's a 1D object in 2D space.
unsigned ngeom_data() const
Number of geometric Data in GeomObject: None.
void position(const unsigned &t, const Vector< double > &xi, Vector< double > &r) const
Return the position along the top of the flag (xi[0] varies between 0 and Lx)
~TopOfFlag()
Destructor (empty)
///////////////////////////////////////////////////////////////// ///////////////////////////////////...
void actions_after_newton_solve()
Update the problem specs after solve (empty)
void actions_before_implicit_timestep()
Update the velocity boundary condition on the flag.
void doc_solution(DocInfo &doc_info)
Doc the solution.
RefineableAlgebraicCylinderWithFlagMesh< ELEMENT > * mesh_pt()
Access function for the specific mesh.
RefineableCylinderWithFlagMesh< ELEMENT > * mesh_pt()
Access function for the specific mesh.
double Height
Height of channel.
void actions_after_adapt()
After adaptation: Unpin pressures and pin redudant pressure dofs.
TurekNonFSIProblem(const double &length, const double &height)
Constructor: Pass length and height of domain.
void actions_before_newton_solve()
Update the problem specs before solve (empty)
Namespace for definition of flag boundaries.
double Centre_y
y position of centre of cylinder
TipOfFlag * Tip_flag_pt
Pointer to GeomObject that bounds the tip edge of the flag.
double Period
Period of prescribed flag oscillation.
double Centre_x
x position of centre of cylinder
void setup(Time *time_pt)
Create all GeomObjects needed to define the cylinder and the flag.
Vector< double > upper_tip(const double &t)
Time-dependent vector to upper tip of the "flag".
double H
Height of flag.
Circle * Cylinder_pt
Pointer to GeomObject of type Circle that defines the central cylinder.
Vector< double > lower_tip(const double &t)
Time-dependent vector to bottom tip of the "flag".
double Radius
Radius of cylinder.
double L
Length of flag.
double Amplitude
Amplitude of tip deflection.
Time * Time_pt
Pointer to the global time object.
BottomOfFlag * Bottom_flag_pt
Pointer to GeomObject that bounds the bottom edge of the flag.
TopOfFlag * Top_flag_pt
Pointer to GeomObject that bounds the upper edge of the flag.
Global parameters.
double Re
Reynolds number.
int main(int argc, char *argv[])
//////////////////////////////////////////////////////////////////// ////////////////////////////////...