channel_with_leaflet.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-2022 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// Generic oomph-lib includes
27#include "generic.h"
28
29// Generic oomph-lib includes
30#include "navier_stokes.h"
31
32//Include the mesh
33#include "meshes/channel_with_leaflet_mesh.h"
34
35using namespace std;
36using namespace oomph;
37
38
39//===start_of_leaflet_class===========================================
40/// GeomObject representing a vertical leaflet that performs
41/// bending and stretching oscillations.
42//====================================================================
43class Leaflet : public GeomObject
44{
45
46public:
47
48 /// Constructor: Pass length (in Lagrangian coordinates),
49 /// the amplitude of the horizontal and vertical deflection of the tip,
50 /// the x-coordinate of the origin and the period of the oscillation.
51 /// Passes the number of Lagrangian and Eulerian coordinates to the
52 /// constructor of the GeomObject base class.
53 Leaflet(const double& length, const double& d_x,const double& d_y,
54 const double& x_0, const double& period, Time* time_pt)
55 : GeomObject(1,2), Length(length), D_x(d_x), D_y(d_y), X_0(x_0),
56 T(period),Time_pt(time_pt) {}
57
58
59 /// Destructor -- emtpy
60 virtual ~Leaflet(){}
61
62 /// Position vector, r, to the point identified by
63 /// its 1D Lagrangian coordinate, xi (passed as a 1D Vector) at discrete time
64 /// level t (t=0: present; t>0: previous).
65 void position(const unsigned& t, const Vector<double>& xi,
66 Vector<double>& r) const
67 {
68 using namespace MathematicalConstants;
69
70 //Position
71 r[0] = X_0 + D_x*xi[0]*xi[0]/Length/Length*sin(2.0*Pi*Time_pt->time(t)/T);
72 r[1] = xi[0]*(1.0+D_y/Length*0.5*(1.0-cos(4.0*Pi*Time_pt->time(t)/T)));
73 }
74
75 /// Steady version: Get current shape
76 void position(const Vector<double>& xi, Vector<double>& r) const
77 {
78 position(0,xi,r);
79 }
80
81 /// Number of geometric Data in GeomObject: None.
82 unsigned ngeom_data() const {return 0;}
83
84 /// Length of the leaflet
85 double length() { return Length; }
86
87 /// Amplitude of horizontal tip displacement
88 double& d_x() {return D_x;}
89
90 /// Amplitude of vertical tip displacement
91 double d_y() {return D_y;}
92
93 /// x-coordinate of leaflet origin
94 double x_0() {return X_0;}
95
96private :
97
98 /// Length in terms of Lagrangian coordinates
99 double Length;
100
101 /// Horizontal displacement of tip
102 double D_x;
103
104 /// Vertical displacement of tip
105 double D_y;
106
107/// Origin
108 double X_0;
109
110 /// Period of the oscillations
111 double T;
112
113 /// Pointer to the global time object
114 Time* Time_pt;
115
116}; //end_of_the_GeomObject
117
118
119
120/// ////////////////////////////////////////////////////////////////////
121/// ///////////////////////////////////////////////////////////////////
122/// ////////////////////////////////////////////////////////////////////
123
124
125//==start_of_global_parameters=======================================
126/// Global parameters
127//===================================================================
129{
130
131 /// Reynolds number
132 double Re=20.0;
133
134} // end_of_namespace
135
136
137
138/// ////////////////////////////////////////////////////////////////////
139/// ///////////////////////////////////////////////////////////////////
140/// ////////////////////////////////////////////////////////////////////
141
142
143//==start_of_problem_class===========================================
144/// Problem class
145//===================================================================
146template<class ELEMENT>
147class ChannelWithLeafletProblem : public Problem
148{
149
150public:
151
152 /// Constructor: Pass the length of the domain at the left
153 /// of the leaflet lleft,the length of the domain at the right of the
154 /// leaflet lright,the height of the leaflet hleaflet, the total height
155 /// of the domain htot, the number of macro-elements at the left of the
156 /// leaflet nleft, the number of macro-elements at the right of the
157 /// leaflet nright, the number of macro-elements under hleaflet ny1,
158 /// the number of macro-elements above hleaflet ny2,the x-displacement
159 /// of the leaflet d_x,the y-displacement of the leaflet d_y,the abscissa
160 /// of the origin of the leaflet x_0, the period of the moving leaflet.
161 ChannelWithLeafletProblem(const double& l_left,
162 const double& l_right, const double& h_leaflet,
163 const double& h_tot,
164 const unsigned& nleft, const unsigned& nright,
165 const unsigned& ny1, const unsigned& ny2,
166 const double& d_x,const double& d_y,
167 const double& x_0, const double& period);
168
169 /// Destructor (empty)
171
172 /// Overloaded access function to specific mesh
173 RefineableAlgebraicChannelWithLeafletMesh<ELEMENT>* mesh_pt()
174 {
175 // Upcast from pointer to the Mesh base class to the specific
176 // element type that we're using here.
177 return dynamic_cast<RefineableAlgebraicChannelWithLeafletMesh<ELEMENT>*>(
178 Problem::mesh_pt());
179 }
180
181 /// Update after solve (empty)
183
184 /// Update before solve (empty)
186
187 /// Actions after adaptation: Pin redundant pressure dofs
188 void actions_after_adapt();
189
190 /// Update the velocity boundary condition on the moving leaflet
192
193 /// Doc the solution
194 void doc_solution(DocInfo& doc_info);
195
196private:
197
198 /// Pointer to the GeomObject
199 GeomObject* Leaflet_pt;
200
201};
202
203
204
205
206//==start_of_constructor=================================================
207/// Constructor
208//=======================================================================
209template <class ELEMENT>
211 const double& l_left,
212 const double& l_right, const double& h_leaflet,
213 const double& h_tot,
214 const unsigned& nleft, const unsigned& nright,
215 const unsigned& ny1, const unsigned& ny2,
216 const double& d_x,const double& d_y,
217 const double& x_0, const double& period)
218{
219 // Allocate the timestepper
220 add_time_stepper_pt(new BDF<2>);
221
222 //Create the geometric object that represents the leaflet
223 Leaflet_pt = new Leaflet(h_leaflet, d_x, d_y, x_0, period, time_pt()) ;
224
225//Build the mesh
226Problem::mesh_pt()=new RefineableAlgebraicChannelWithLeafletMesh<ELEMENT>(
227 Leaflet_pt,
228 l_left, l_right,
229 h_leaflet,
230 h_tot,nleft,
231 nright,ny1,ny2,
232 time_stepper_pt());
233
234
235 // Set error estimator
236 Z2ErrorEstimator* error_estimator_pt=new Z2ErrorEstimator;
237 dynamic_cast<RefineableAlgebraicChannelWithLeafletMesh<ELEMENT>*>(mesh_pt())->
238 spatial_error_estimator_pt()=error_estimator_pt;
239
240
241 // Loop over the elements to set up element-specific
242 // things that cannot be handled by constructor
243 unsigned n_element = mesh_pt()->nelement();
244 for(unsigned e=0;e<n_element;e++)
245 {
246 // Upcast from GeneralisedElement to the present element
247 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
248
249 //Set the Reynolds number
250 el_pt->re_pt() = &Global_Physical_Variables::Re;
251
252 // Set the Womersley number (product of Reynolds and Strouhal).
253 // We're assuming a Strouhal number of one, corresponding to
254 // a non-dimensionalisation of time on the flow's natural timescale.
255 el_pt->re_st_pt() = &Global_Physical_Variables::Re;
256
257 } // end loop over elements
258
259
260 //Pin the boundary nodes
261 unsigned num_bound = mesh_pt()->nboundary();
262 for(unsigned ibound=0;ibound<num_bound;ibound++)
263 {
264 unsigned num_nod= mesh_pt()->nboundary_node(ibound);
265 for (unsigned inod=0;inod<num_nod;inod++)
266 {
267 mesh_pt()->boundary_node_pt(ibound,inod)->pin(1);
268 //do not pin the x velocity of the outflow
269 if( ibound != 1)
270 {
271 mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
272 }
273 }
274 }
275
276 // Setup parabolic flow along the inflow boundary 3
277 unsigned ibound=3;
278 unsigned num_nod= mesh_pt()->nboundary_node(ibound);
279 for (unsigned inod=0;inod<num_nod;inod++)
280 {
281 double ycoord = mesh_pt()->boundary_node_pt(ibound,inod)->x(1);
282 double uy = 6.0*(ycoord/h_tot)*(1.0-(ycoord/h_tot));
283
284 mesh_pt()->boundary_node_pt(ibound,inod)->set_value(0,uy);
285 mesh_pt()->boundary_node_pt(ibound,inod)->set_value(1,0.0);
286 }// end of setup boundary condition
287
288 // Pin redudant pressure dofs
289 RefineableNavierStokesEquations<2>::
290 pin_redundant_nodal_pressures(Problem::mesh_pt()->element_pt());
291
292 // Setup equation numbering scheme
293 cout << "Number of equations: " << assign_eqn_numbers() << std::endl;
294
295}//end of constructor
296
297
298
299
300//=====start_of_actions_before_implicit_timestep=========================
301/// Actions before implicit timestep: Update domain shape and
302/// the velocity boundary conditions
303//=======================================================================
304template <class ELEMENT>
306{
307 // Update the domain shape
308 mesh_pt()->node_update();
309
310 // Moving leaflet: No slip; this implies that the velocity needs
311 // to be updated in response to leaflet motion
312 for( unsigned ibound=4;ibound<6;ibound++)
313 {
314 unsigned num_nod=mesh_pt()->nboundary_node(ibound);
315 for (unsigned inod=0;inod<num_nod;inod++)
316 {
317 // Which node are we dealing with?
318 Node* node_pt=mesh_pt()->boundary_node_pt(ibound,inod);
319
320 // Apply no slip
321 FSI_functions::apply_no_slip_on_moving_wall(node_pt);
322 }
323 }
324} //end_of_actions_before_implicit_timestep
325
326
327
328//==========start_of_actions_after_adaptation============================
329// Actions after adaptation: Pin redundant pressure dofs
330//=======================================================================
331template<class ELEMENT>
333{
334 // Unpin all pressure dofs
335 RefineableNavierStokesEquations<2>::
336 unpin_all_pressure_dofs(mesh_pt()->element_pt());
337
338 // Pin redundant pressure dofs
339 RefineableNavierStokesEquations<2>::
340 pin_redundant_nodal_pressures(mesh_pt()->element_pt());
341
342} // end_of_actions_after_adapt
343
344
345
346//==start_of_doc_solution=================================================
347/// Doc the solution
348//========================================================================
349template<class ELEMENT>
351{
352
353 ofstream some_file;
354 char filename[100];
355
356 // Number of plot points
357 unsigned npts;
358 npts=5;
359
360 // Output solution
361 sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
362 doc_info.number());
363 some_file.open(filename);
364 mesh_pt()->output(some_file,npts);
365 some_file.close();
366
367 // Output boundaries
368 sprintf(filename,"%s/boundaries%i.dat",doc_info.directory().c_str(),
369 doc_info.number());
370 some_file.open(filename);
371 mesh_pt()->output_boundaries(some_file);
372 some_file.close();
373
374} // end_of_doc_solution
375
376
377/// ////////////////////////////////////////////////////////////////////
378/// ///////////////////////////////////////////////////////////////////
379/// ////////////////////////////////////////////////////////////////////
380
381
382//======start_of_main==================================================
383/// Driver code -- pass a command line argument if you want to run
384/// the code in validation mode where it only performs a few steps
385//=====================================================================
386int main(int argc, char* argv[])
387{
388
389 // Store command line arguments
390 CommandLineArgs::setup(argc,argv);
391
392 // Set up doc info
393 DocInfo doc_info;
394 doc_info.set_directory("RESLT");
395 doc_info.number()=0;
396
397 // Parameters for the leaflet
398 //----------------------------
399
400 // Height
401 double h_leaflet = 0.5;
402
403
404 // Tip deflection
405 double d_x = 0.25;
406 double d_y = -0.05;
407
408 // x-positon of root
409 double x_0 = 3.0;
410
411 // Period of the oscillation on the natural timescale of the flow
412 double period = 20.0;
413
414
415 //Parameters for the domain
416 //-------------------------
417
418 // Length of the mesh to right and left of the leaflet
419 double l_left =2.0;
420 double l_right= 3.0;
421
422 // Total height of domain (unity because lengths have been scaled on it)
423 double h_tot=1.0;
424
425 // Initial number of element rows/columns in various mesh regions
426 unsigned nleft=8;
427 unsigned nright=12;
428 unsigned ny1=2;
429 unsigned ny2=2;
430
431 //Build the problem
433 problem(l_left, l_right, h_leaflet,
434 h_tot,nleft, nright,ny1,ny2,
435 d_x, d_y, x_0,
436 period);
437
438
439 // Number of timesteps per period
440 unsigned nsteps_per_period=40;
441
442 // Number of periods
443 unsigned nperiod=3;
444
445 // Number of timesteps (reduced for validation)
446 unsigned nstep=nsteps_per_period*nperiod;
447 if (CommandLineArgs::Argc>1)
448 {
449 nstep=3;
450 }
451
452 //Timestep:
453 double dt=period/double(nsteps_per_period);
454
455 /// Initialise timestep
456 problem.initialise_dt(dt);
457
458
459 /// Set max. number of adaptations (reduced for validation)
460 unsigned max_adapt=5;
461 if (CommandLineArgs::Argc>1)
462 {
463 max_adapt=2;
464 }
465
466 // Do steady solve first -- this also sets the history values
467 // to those corresponding to an impulsive start from the
468 // steady solution
469 problem.steady_newton_solve(max_adapt);
470
471 /// Output steady solution
472 problem.doc_solution(doc_info);
473 doc_info.number()++;
474
475
476 /// Reduce the max number of adaptations for time-dependent simulation
477 max_adapt=1;
478
479 // We don't want to re-assign the initial condition
480 bool first=false;
481
482// Timestepping loop
483 for (unsigned istep=0;istep<nstep;istep++)
484 {
485 // Solve the problem
486 problem.unsteady_newton_solve(dt,max_adapt,first);
487
488 // Output the solution
489 problem.doc_solution(doc_info);
490
491 // Step number
492 doc_info.number()++;
493
494 }
495
496
497}//end of main
498
499
int main(int argc, char *argv[])
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void actions_after_newton_solve()
Update after solve (empty)
ChannelWithLeafletProblem(const double &l_left, const double &l_right, const double &h_leaflet, const double &h_tot, const unsigned &nleft, const unsigned &nright, const unsigned &ny1, const unsigned &ny2, const double &d_x, const double &d_y, const double &x_0, const double &period)
Constructor: Pass the length of the domain at the left of the leaflet lleft,the length of the domain ...
void actions_before_newton_solve()
Update before solve (empty)
~ChannelWithLeafletProblem()
Destructor (empty)
void actions_after_adapt()
Actions after adaptation: Pin redundant pressure dofs.
GeomObject * Leaflet_pt
Pointer to the GeomObject.
void actions_before_implicit_timestep()
Update the velocity boundary condition on the moving leaflet.
RefineableAlgebraicChannelWithLeafletMesh< ELEMENT > * mesh_pt()
Overloaded access function to specific mesh.
void doc_solution(DocInfo &doc_info)
Doc the solution.
GeomObject representing a vertical leaflet that performs bending and stretching oscillations.
void position(const Vector< double > &xi, Vector< double > &r) const
Steady version: Get current shape.
virtual ~Leaflet()
Destructor – emtpy.
double T
Period of the oscillations.
double X_0
Origin.
double D_x
Horizontal displacement of tip.
double & d_x()
Amplitude of horizontal tip displacement.
void position(const unsigned &t, const Vector< double > &xi, Vector< double > &r) const
Position vector, r, to the point identified by its 1D Lagrangian coordinate, xi (passed as a 1D Vec...
double Length
Length in terms of Lagrangian coordinates.
double length()
Length of the leaflet.
double d_y()
Amplitude of vertical tip displacement.
unsigned ngeom_data() const
Number of geometric Data in GeomObject: None.
Time * Time_pt
Pointer to the global time object.
Leaflet(const double &length, const double &d_x, const double &d_y, const double &x_0, const double &period, Time *time_pt)
Constructor: Pass length (in Lagrangian coordinates), the amplitude of the horizontal and vertical de...
double x_0()
x-coordinate of leaflet origin
double D_y
Vertical displacement of tip.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...