jeffery_orbit.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// Driver to demonstrate the interaction of a fluid flow with a rigid body.
27// The driver solves a classic problem of the motion of an rigid ellipse in
28// a shear flow. The problem imposes a smooth, but rapid, start-up from
29// no flow to a shear imposed on the boundaries of a rectangular domain
30// and the ellipse soon settles into approximate Jeffery orbits for small
31// enough Reynolds number.
32
33//Generic routines
34#include "generic.h"
35
36// The equations
37#include "navier_stokes.h"
38#include "solid.h"
39#include "constitutive.h"
40#include "rigid_body.h"
41
42// The mesh
43#include "meshes/triangle_mesh.h"
44
45//Thin wrapper to "custom" TaylorHoodElements that overload output functions
47
48#include <algorithm>
49
50using namespace std;
51using namespace oomph;
52
53//==start_of_namespace==============================
54/// Namespace for Problem Parameters
55//==================================================
57 {
58 /// Reynolds number
59 double Re=1.0;
60
61 /// Strouhal number
62 double St = 1.0;
63
64 /// Density ratio (Solid density / Fluid density)
65 double Density_ratio = 1.0;
66
67 /// Initial axis of the elliptical solid in x-direction
68 double A = 0.25;
69
70 /// Initial axis of the elliptical solid in y-direction
71 /// (N.B. 2B = 1 is the reference length scale)
72 double B = 0.5;
73
74 /// Pseudo-solid (mesh) Poisson ratio
75 double Nu=0.3;
76
77 /// Pseudo-solid (mesh) "density"
78 /// Set to zero because we don't want inertia in the node update!
79 double Lambda_sq=0.0;
80
81 /// Constitutive law used to determine the mesh deformation
82 ConstitutiveLaw *Constitutive_law_pt=
83 new GeneralisedHookean(&Problem_Parameter::Nu);
84
85 } // end_of_namespace
86
87
88//=======================================================================
89/// Exact solution for the rotation of an ellipse in unbounded shear flow
90/// as computed by Jeffery (1922)
91//=======================================================================
93{
94 /// Null function
95 double null(const double &t) {return 0.0;}
96
97 /// Angular position as a function of time t
98 double angle(const double &t)
99 {
100 const double a = Problem_Parameter::A;
101 const double b = Problem_Parameter::B;
102
103 return atan((b/a)*tan((a*b*t)/(b*b + a*a)));
104 }
105
106 /// Angular velocity as function of time t
107 double velocity(const double &t)
108 {
109 const double a = Problem_Parameter::A;
110 const double b = Problem_Parameter::B;
111
112 //Get the angle
113 double chi = angle(t);
114
115 //Now return the velocity
116 return (a*a*sin(chi)*sin(chi) + b*b*cos(chi)*cos(chi))/(a*a + b*b);
117 }
118
119 /// Angular acceleration as a function of time t (should always be zero)
120 double acceleration(const double &t)
121 {
122 const double a = Problem_Parameter::A;
123 const double b = Problem_Parameter::A;
124
125 //Get the angle and velocity
126 double chi = angle(t);
127 double chi_dot = velocity(t);
128
129 //Now return the acceleration
130 return 2.0*(a*a - b*b)*(sin(chi)*cos(chi))*chi_dot/(a*a + b*b);
131 }
132}
133
134/// ////////////////////////////////////////////////////////
135/// ////////////////////////////////////////////////////////
136/// ////////////////////////////////////////////////////////
137
138
139//===================start_of_general_ellipse=================================
140/// A geometric object for an ellipse with initial centre of mass at
141/// (centre_x, centre_y) with axis in the x direction given by 2a
142/// and in the y-direction given by 2b. The boundary of the ellipse is
143/// parametrised by its angle.
144//============================================================================
145class GeneralEllipse : public GeomObject
146{
147private:
148
149 //Storage for the centre of mass and semi-major and semi-minor axes
150 double Centre_x, Centre_y, A, B;
151
152public:
153
154 /// Simple Constructor that transfers appropriate geometric
155 /// parameters into internal data
156 GeneralEllipse(const double &centre_x, const double &centre_y,
157 const double &a, const double &b)
158 : GeomObject(1,2), Centre_x(centre_x), Centre_y(centre_y), A(a), B(b)
159 {}
160
161 /// Empty Destructor
163
164 /// Return the position of the ellipse boundary as a function of
165 /// the angle xi[0]
166 void position(const Vector<double> &xi, Vector<double> &r) const
167 {
168 r[0] = Centre_x + A*cos(xi[0]);
169 r[1] = Centre_y + B*sin(xi[0]);
170 }
171
172 //Return the position which is always fixed
173 void position(const unsigned &t,
174 const Vector<double> &xi, Vector<double> &r) const
175 {
176 return position(xi,r);
177 }
178
179};
180//end_of_general_ellipse
181
182
183/// ////////////////////////////////////////////////////////
184/// ////////////////////////////////////////////////////////
185/// ////////////////////////////////////////////////////////
186
187
188//==start_of_problem_class============================================
189/// Unstructured Navier-Stokes ALE Problem for a rigid ellipse
190/// immersed within a viscous fluid
191//====================================================================
192template<class ELEMENT>
194{
195
196public:
197
198 /// Constructor
200
201 /// Destructor
203
204 /// Reset the boundary conditions when timestepping
206 {
207 this->set_boundary_velocity();
208 }
209
210 /// Wipe the meshes of Lagrange multiplier and drag elements
212
213 /// Rebuild the meshes of Lagrange multiplier and drag elements
214 void actions_after_adapt();
215
216 /// Re-apply the no slip condition (imposed indirectly via dependent
217 /// velocities)
219 {
220 // Update mesh -- this applies the auxiliary node update function
221 Fluid_mesh_pt->node_update();
222 }
223
224 /// Set boundary condition, assign auxiliary node update fct.
225 /// Complete the build of all elements, attach power elements that allow
226 /// computation of drag vector
228
229 /// Set the boundary velocity
231
232 /// Function that solves a simplified problem to ensure that
233 /// the positions of the boundary nodes are initially consistent with
234 /// the lagrange multiplier formulation
236
237 /// Doc the solution
238 void doc_solution(const bool& project=false);
239
240 /// Output the exact solution
241 void output_exact_solution(std::ofstream &output_file);
242
243private:
244
245 /// Create elements that enforce prescribed boundary motion
246 /// for the pseudo-solid fluid mesh by Lagrange multipliers
248
249 /// Delete elements that impose the prescribed boundary displacement
250 /// and wipe the associated mesh
252
253 /// Create elements that calculate the drag and torque on
254 /// the boundaries
256
257 /// Delete elements that calculate the drag and torque on the
258 /// boundaries
260
261 /// Pin the degrees of freedom associated with the solid bodies
262 void pin_rigid_body();
263
264 /// Unpin the degrees of freedom associated with the solid bodies
265 void unpin_rigid_body();
266
267 /// Pointers to mesh of Lagrange multiplier elements
269
270 /// Pointer to Fluid_mesh
271 RefineableSolidTriangleMesh<ELEMENT>* Fluid_mesh_pt;
272
273 /// Triangle mesh polygon for outer boundary
274 TriangleMeshPolygon* Outer_boundary_polygon_pt;
275
276 /// Mesh of drag elements
277 Vector<Mesh*> Drag_mesh_pt;
278
279 /// Mesh of the generalised elements for the rigid bodies
281
282 /// Storage for the geom object
283 Vector<GeomObject*> Rigid_body_pt;
284
285 /// Internal DocInfo object
286 DocInfo Doc_info;
287
288 /// File to document the norm of the solution (for validation purposes)
289 ofstream Norm_file;
290
291 /// File to document the motion of the centre of gravity
292 ofstream Cog_file;
293
294 /// File to document the exact motion of the centre of gravity
296
297}; // end_of_problem_class
298
299
300//==start_constructor=====================================================
301/// Constructor: Open output files, construct time steppers, build
302/// fluid mesh, immersed rigid body and combine to form the problem
303//========================================================================
304template<class ELEMENT>
307{
308 // Output directory
309 this->Doc_info.set_directory("RESLT");
310
311 // Open norm file
312 this->Norm_file.open("RESLT/norm.dat");
313
314 // Open file to trace the centre of gravity
315 this->Cog_file.open("RESLT/cog_trace.dat");
316
317 // Open file to document the exact motion of the centre of gravity
318 this->Cog_exact_file.open("RESLT/cog_exact_trace.dat");
319
320 // Allocate the timestepper -- this constructs the Problem's
321 // time object with a sufficient amount of storage to store the
322 // previous timsteps.
323 this->add_time_stepper_pt(new BDF<2>);
324
325 // Allocate a timestepper for the rigid body
326 this->add_time_stepper_pt(new Newmark<2>);
327
328 // Define the boundaries: Polyline with 4 different
329 // boundaries for the outer boundary and 1 internal elliptical hole
330
331 // Build the boundary segments for outer boundary, consisting of
332 //--------------------------------------------------------------
333 // four separate polyline segments
334 //---------------------------------
335 Vector<TriangleMeshCurveSection*> boundary_segment_pt(4);
336
337 //Set the length of the channel
338 double half_length = 5.0;
339 double half_height = 2.5;
340
341 // Initialize boundary segment
342 Vector<Vector<double> > bound_seg(2);
343 for(unsigned i=0;i<2;i++) {bound_seg[i].resize(2);}
344
345 // First boundary segment
346 bound_seg[0][0]=-half_length;
347 bound_seg[0][1]=-half_height;
348 bound_seg[1][0]=-half_length;
349 bound_seg[1][1]=half_height;
350
351 // Specify 1st boundary id
352 unsigned bound_id = 0;
353
354 // Build the 1st boundary segment
355 boundary_segment_pt[0] = new TriangleMeshPolyLine(bound_seg,bound_id);
356
357 // Second boundary segment
358 bound_seg[0][0]=-half_length;
359 bound_seg[0][1]=half_height;
360 bound_seg[1][0]=half_length;
361 bound_seg[1][1]=half_height;
362
363 // Specify 2nd boundary id
364 bound_id = 1;
365
366 // Build the 2nd boundary segment
367 boundary_segment_pt[1] = new TriangleMeshPolyLine(bound_seg,bound_id);
368
369 // Third boundary segment
370 bound_seg[0][0]=half_length;
371 bound_seg[0][1]=half_height;
372 bound_seg[1][0]=half_length;
373 bound_seg[1][1]=-half_height;
374
375 // Specify 3rd boundary id
376 bound_id = 2;
377
378 // Build the 3rd boundary segment
379 boundary_segment_pt[2] = new TriangleMeshPolyLine(bound_seg,bound_id);
380
381 // Fourth boundary segment
382 bound_seg[0][0]=half_length;
383 bound_seg[0][1]=-half_height;
384 bound_seg[1][0]=-half_length;
385 bound_seg[1][1]=-half_height;
386
387 // Specify 4th boundary id
388 bound_id = 3;
389
390 // Build the 4th boundary segment
391 boundary_segment_pt[3] = new TriangleMeshPolyLine(bound_seg,bound_id);
392
393 // Create the triangle mesh polygon for outer boundary using boundary segment
394 Outer_boundary_polygon_pt = new TriangleMeshPolygon(boundary_segment_pt);
395
396 // Now build the moving rigid body
397 //-------------------------------------
398
399 // We have one rigid body
400 Rigid_body_pt.resize(1);
401 Vector<TriangleMeshClosedCurve*> hole_pt(1);
402
403 // Build Rigid Body
404 //-----------------
405 double x_center = 0.0;
406 double y_center = 0.0;
407 double A = Problem_Parameter::A;
408 double B = Problem_Parameter::B;
409 GeomObject* temp_hole_pt = new GeneralEllipse(x_center,y_center,A,B);
410 Rigid_body_pt[0] = new ImmersedRigidBodyElement(temp_hole_pt,
411 this->time_stepper_pt(1));
412
413 // Build the two parts of the curvilinear boundary from the rigid body
414 Vector<TriangleMeshCurveSection*> curvilinear_boundary_pt(2);
415
416 //First section (boundary 4)
417 double zeta_start=0.0;
418 double zeta_end=MathematicalConstants::Pi;
419 unsigned nsegment=8;
420 unsigned boundary_id=4;
421 curvilinear_boundary_pt[0]=new TriangleMeshCurviLine(
422 Rigid_body_pt[0],zeta_start,zeta_end,nsegment,boundary_id);
423
424 //Second section (boundary 5)
425 zeta_start=MathematicalConstants::Pi;
426 zeta_end=2.0*MathematicalConstants::Pi;
427 nsegment=8;
428 boundary_id=5;
429 curvilinear_boundary_pt[1]=new TriangleMeshCurviLine(
430 Rigid_body_pt[0],zeta_start,zeta_end,
431 nsegment,boundary_id);
432
433 // Combine to form a hole in the fluid mesh
434 Vector<double> hole_coords(2);
435 hole_coords[0]=0.0;
436 hole_coords[1]=0.0;
437 Vector<TriangleMeshClosedCurve*> curvilinear_hole_pt(1);
438 hole_pt[0]=
439 new TriangleMeshClosedCurve(
440 curvilinear_boundary_pt,hole_coords);
441
442 // Now build the mesh, based on the boundaries specified by
443 //---------------------------------------------------------
444 // polygons just created
445 //----------------------
446
447 TriangleMeshClosedCurve* closed_curve_pt=Outer_boundary_polygon_pt;
448
449 double uniform_element_area=1.0;
450
451 // Use the TriangleMeshParameters object for gathering all
452 // the necessary arguments for the TriangleMesh object
453 TriangleMeshParameters triangle_mesh_parameters(
454 closed_curve_pt);
455
456 // Define the holes on the domain
457 triangle_mesh_parameters.internal_closed_curve_pt() =
458 hole_pt;
459
460 // Define the maximum element area
461 triangle_mesh_parameters.element_area() =
462 uniform_element_area;
463
464 // Create the mesh
465 Fluid_mesh_pt =
466 new RefineableSolidTriangleMesh<ELEMENT>(
467 triangle_mesh_parameters, this->time_stepper_pt());
468
469 // Set error estimator for bulk mesh
470 Z2ErrorEstimator* error_estimator_pt=new Z2ErrorEstimator;
471 Fluid_mesh_pt->spatial_error_estimator_pt()=error_estimator_pt;
472
473 // Set targets for spatial adaptivity
474 Fluid_mesh_pt->max_permitted_error()=0.005;
475 Fluid_mesh_pt->min_permitted_error()=0.001;
476 Fluid_mesh_pt->max_element_size()=1.0;
477 Fluid_mesh_pt->min_element_size()=0.001;
478
479 // Use coarser mesh during validation
480 if (CommandLineArgs::command_line_flag_has_been_set("--validation"))
481 {
482 Fluid_mesh_pt->min_element_size()=0.01;
483 }
484
485 // Set boundary condition, assign auxiliary node update fct,
486 // complete the build of all elements, attach power elements that allow
487 // computation of drag vector
488 complete_problem_setup();
489
490 //Set the parameters of the rigid body elements
491 ImmersedRigidBodyElement* rigid_element1_pt =
492 dynamic_cast<ImmersedRigidBodyElement*>(Rigid_body_pt[0]);
493 rigid_element1_pt->initial_centre_of_mass(0) = x_center;
494 rigid_element1_pt->initial_centre_of_mass(1) = y_center;
495 rigid_element1_pt->mass_shape() = MathematicalConstants::Pi*A*B;
496 rigid_element1_pt->moment_of_inertia_shape() =
497 0.25*MathematicalConstants::Pi*A*B*(A*A + B*B);
498 rigid_element1_pt->re_pt() = &Problem_Parameter::Re;
499 rigid_element1_pt->st_pt() = &Problem_Parameter::St;
500 rigid_element1_pt->density_ratio_pt() = &Problem_Parameter::Density_ratio;
501
502 //Pin the position of the centre of mass
503 rigid_element1_pt->pin_centre_of_mass_coordinate(0);
504 rigid_element1_pt->pin_centre_of_mass_coordinate(1);
505
506 // Create the mesh for the rigid bodies
507 Rigid_body_mesh_pt = new Mesh;
508 Rigid_body_mesh_pt->add_element_pt(rigid_element1_pt);
509
510 // Create the drag mesh for the rigid bodies
511 Drag_mesh_pt.resize(1);
512 for(unsigned m=0;m<1;m++) {Drag_mesh_pt[m] = new Mesh;}
513 this->create_drag_elements();
514
515 //Add the drag mesh to the appropriate rigid bodies
516 rigid_element1_pt->set_drag_mesh(Drag_mesh_pt[0]);
517
518
519
520 // Create Lagrange multiplier mesh for boundary motion
521 //----------------------------------------------------
522 // Construct the mesh of elements that enforce prescribed boundary motion
523 // of pseudo-solid fluid mesh by Lagrange multipliers
524 Lagrange_multiplier_mesh_pt=new SolidMesh;
525 create_lagrange_multiplier_elements();
526
527
528 // Combine meshes
529 //---------------
530
531 // Add Fluid_mesh_pt sub meshes
532 this->add_sub_mesh(Fluid_mesh_pt);
533
534 // Add Lagrange_multiplier sub meshes
535 this->add_sub_mesh(this->Lagrange_multiplier_mesh_pt);
536
537 this->add_sub_mesh(this->Rigid_body_mesh_pt);
538
539 // Build global mesh
540 this->build_global_mesh();
541
542 // Setup equation numbering scheme
543 cout <<"Number of equations: " << this->assign_eqn_numbers() << std::endl;
544
545} // end_of_constructor
546
547
548//========================================================================
549/// Destructor that cleans up memory and closes files
550//========================================================================
551template<class ELEMENT>
554{
555 // Delete Fluid timestepper
556 delete this->time_stepper_pt(0);
557
558 //Delete solid timestepper
559 delete this->time_stepper_pt(1);
560
561 // Kill data associated with outer boundary
562 unsigned n=Outer_boundary_polygon_pt->npolyline();
563 for (unsigned j=0;j<n;j++)
564 {
565 delete Outer_boundary_polygon_pt->polyline_pt(j);
566 }
567 delete Outer_boundary_polygon_pt;
568
569 //Flush the drag mesh from the rigid body
570 dynamic_cast<ImmersedRigidBodyElement*>(Rigid_body_pt[0])->
571 flush_drag_mesh();
572
573 // Flush Lagrange multiplier mesh
574 delete_lagrange_multiplier_elements();
575 delete Lagrange_multiplier_mesh_pt;
576
577 // Delete error estimator
578 delete Fluid_mesh_pt->spatial_error_estimator_pt();
579
580 // Delete fluid mesh
581 delete Fluid_mesh_pt;
582
583 // Kill const eqn
585
586 // Close norm and trace files
587 this->Cog_exact_file.close();
588 this->Cog_file.close();
589 this->Norm_file.close();
590}
591
592
593//==========start_actions_before_adapt==================================
594/// Actions before adapt: Wipe the mesh of Lagrange multiplier elements
595//=======================================================================
596template<class ELEMENT>
598{
599 //Flush the drag mesh from the rigid body
600 dynamic_cast<ImmersedRigidBodyElement*>(Rigid_body_pt[0])->
601 flush_drag_mesh();
602
603 //Kill the drag element
604 delete_drag_elements();
605
606 // Kill the elements and wipe surface mesh
607 delete_lagrange_multiplier_elements();
608
609 // Rebuild the Problem's global mesh from its various sub-meshes
610 this->rebuild_global_mesh();
611} // end of actions_before_adapt
612
613
614//============start_actions_after_adapt=====================================
615/// Actions after adapt: Rebuild the mesh of Lagrange multiplier elements
616//==========================================================================
617template<class ELEMENT>
619{
620 //Reset the lagrangian coordinates for the solid mechanics
621 //an updated lagrangian approach
622 Fluid_mesh_pt->set_lagrangian_nodal_coordinates();
623
624 // Create the elements that impose the displacement constraint
625 create_lagrange_multiplier_elements();
626
627 // Create the drag elements anew
628 create_drag_elements();
629
630 // Add the drag elements to the rigid body
631 dynamic_cast<ImmersedRigidBodyElement*>(Rigid_body_pt[0])->
632 set_drag_mesh(this->Drag_mesh_pt[0]);
633
634 // Rebuild the Problem's global mesh from its various sub-meshes
635 this->rebuild_global_mesh();
636
637 // Setup the problem again -- remember that fluid mesh has been
638 // completely rebuilt and its element's don't have any
639 // pointers to Re etc. yet
640 complete_problem_setup();
641
642 // Output solution after adaptation/projection
643 bool doc_projection=true;
644 doc_solution(doc_projection);
645}// end of actions_after_adapt
646
647
648//============start_complete_problem_setup=================================
649/// Set boundary condition, assign auxiliary node update fct.
650/// Complete the build of all elements, attach power elements that allow
651/// computation of drag vector
652//=========================================================================
653template<class ELEMENT>
655{
656 // Set the boundary conditions for fluid problem: All nodes are
657 // free by default -- just pin the ones that have Dirichlet conditions
658 // here.
659 unsigned nbound=Fluid_mesh_pt->nboundary();
660 for(unsigned ibound=0;ibound<nbound;ibound++)
661 {
662 unsigned num_nod=Fluid_mesh_pt->nboundary_node(ibound);
663 for (unsigned inod=0;inod<num_nod;inod++)
664 {
665 // Cache pointer to node
666 Node* const nod_pt=Fluid_mesh_pt->boundary_node_pt(ibound,inod);
667
668 //Pin x-velocity unless on inlet (0) and outlet (2) boundaries
669 //of the external rectangular box
670 if((ibound!=0) && (ibound!=2)) {nod_pt->pin(0);}
671 //Pin the y-velocity on all boundaries
672 nod_pt->pin(1);
673
674 // Pin pseudo-solid positions apart from on the
675 // ellipse boundary that is allowed to move
676 // Cache cast pointer to solid node
677 SolidNode* const solid_node_pt = dynamic_cast<SolidNode*>(nod_pt);
678
679 //Pin the solid positions on all external boundaries
680 if(ibound < 4)
681 {
682 solid_node_pt->pin_position(0);
683 solid_node_pt->pin_position(1);
684 }
685 // Unpin the position of all the nodes on hole boundaries:
686 // since they will be moved using Lagrange Multiplier
687 else
688 {
689 solid_node_pt->unpin_position(0);
690 solid_node_pt->unpin_position(1);
691
692 // Assign auxiliary node update fct, which determines the
693 // velocity on the moving boundary using the position history
694 // values
695 // A more accurate version may be obtained by using velocity
696 // based on the actual position of the geometric object,
697 // but this introduces additional dependencies between the
698 // Data of the rigid body and the fluid elements.
699 nod_pt->set_auxiliary_node_update_fct_pt(
700 FSI_functions::apply_no_slip_on_moving_wall);
701 }
702 } //End of loop over boundary nodes
703 } // End loop over boundaries
704
705 // Complete the build of all elements so they are fully functional
706 unsigned n_element = Fluid_mesh_pt->nelement();
707 for(unsigned e=0;e<n_element;e++)
708 {
709 // Upcast from GeneralisedElement to the present element
710 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(Fluid_mesh_pt->element_pt(e));
711
712 // Set the Reynolds number
713 el_pt->re_pt() = &Problem_Parameter::Re;
714
715 // Set the Womersley number (same as Re since St=1)
716 el_pt->re_st_pt() = &Problem_Parameter::Re;
717
718 // Set the constitutive law for pseudo-elastic mesh deformation
719 el_pt->constitutive_law_pt()=Problem_Parameter::Constitutive_law_pt;
720
721 // Set the "density" for pseudo-elastic mesh deformation
722 el_pt->lambda_sq_pt()=&Problem_Parameter::Lambda_sq;
723 }
724
725 // Re-apply Dirichlet boundary conditions for current and history values
726 // (projection ignores boundary conditions!)
727 this->set_boundary_velocity();
728
729} //end_of_complete_problem_setup
730
731
732//=================start_set_boundary_velocity===========================
733/// Set the boundary velocity for current and history values
734//=======================================================================
735template<class ELEMENT>
738{
739 //Loop over top and bottom walls, inlet and outlet
740 for(unsigned ibound=0;ibound<4;ibound++)
741 {
742 //If we are on the inlet or outlet only zero the
743 //y velocity, leave x alone
744 if((ibound==0) || (ibound==2))
745 {
746 unsigned num_nod=this->Fluid_mesh_pt->nboundary_node(ibound);
747 for (unsigned inod=0;inod<num_nod;inod++)
748 {
749 // Get node
750 Node* const nod_pt=this->Fluid_mesh_pt->boundary_node_pt(ibound,inod);
751
752 // Get number of previous (history) values
753 unsigned n_prev=nod_pt->time_stepper_pt()->nprev_values();
754
755 //Zero all current and previous y-velocities
756 for(unsigned t=0;t<=n_prev;t++)
757 {
758 nod_pt->set_value(t,1,0.0);
759 }
760 }
761 }
762 //Otherwise on the top and bottom walls set a smooth x-velocity
763 //and zero y-velocity
764 else
765 {
766 unsigned num_nod=this->Fluid_mesh_pt->nboundary_node(ibound);
767 for (unsigned inod=0;inod<num_nod;inod++)
768 {
769 // Get node
770 Node* nod_pt=this->Fluid_mesh_pt->boundary_node_pt(ibound,inod);
771
772 // Get number of previous (history) values
773 unsigned n_prev=nod_pt->time_stepper_pt()->nprev_values();
774
775 //Now set the boundary velocity
776 double y = nod_pt->x(1);
777 //Get the previous time
778 for(unsigned t=0;t<=n_prev;t++)
779 {
780 //Get the time
781 double time_ = this->time_pt()->time(t);
782
783 //Get the velocity ramp
784 //Initially zero (nothing at all is going on)
785 double ramp = 0.0;
786 double delta = 5.0;
787
788 double e1 = exp(-delta);
789 double a1 = 1.0 - (1.0 + delta + 0.5*delta*delta)*e1;
790 double b1 = (3.0 + 3.0*delta + delta*delta)*e1 - 3.0;
791 double c1 = 3.0 - (3.0 + 2.0*delta + 0.5*delta*delta)*e1;
792 //Smooth start
793 if((time_ >= 0.0) & (time_ <= 1.0))
794 {
795 ramp = a1*time_*time_*time_
796 + b1*time_*time_
797 + c1*time_;
798 }
799 //Coupled to exponential levelling
800 else if (time_ > 1.0)
801 {
802 ramp = 1.0 - exp(-delta*time_);
803 }
804
805 nod_pt->set_value(t,0,-y*ramp);
806 }
807
808 // Zero all current and previous veloc values
809 // for the v-velocity
810 for (unsigned t=0;t<=n_prev;t++)
811 {
812 nod_pt->set_value(t,1,0.0);
813 }
814 }
815 }
816 }
817}
818
819
820//====================start_of_pin_rigid_body=====================
821/// Pin the degrees of freedom associated with the solid bodies
822//================================================================
823template<class ELEMENT>
825{
826 unsigned n_rigid_body = Rigid_body_pt.size();
827 for(unsigned i=0;i<n_rigid_body;++i)
828 {
829 unsigned n_geom_data = Rigid_body_pt[i]->ngeom_data();
830 for(unsigned j=0;j<n_geom_data;j++)
831 {
832 Rigid_body_pt[i]->geom_data_pt(j)->pin_all();
833 }
834 }
835}
836
837
838//=================start_unpin_rigid_body==========================
839/// Unpin the degrees of freedom associated with the solid bodies
840//=================================================================
841template<class ELEMENT>
843{
844 unsigned n_rigid_body = Rigid_body_pt.size();
845 for(unsigned i=0;i<n_rigid_body;++i)
846 {
847 unsigned n_geom_data = Rigid_body_pt[i]->ngeom_data();
848 for(unsigned j=0;j<n_geom_data;j++)
849 {
850 Rigid_body_pt[i]->geom_data_pt(j)->unpin_all();
851 }
852 }
853}
854
855
856//==========start_solve_for_consistent_nodal_positions================
857/// Assemble and solve a simplified problem that ensures that the
858/// positions of the boundary nodes are consistent with the weak
859/// imposition of the displacement boundary conditions on the surface
860/// of the ellipse.
861//===================================================================
862template<class ELEMENT>
865{
866 //First pin all degrees of freedom in the rigid body
867 this->pin_rigid_body();
868
869 //Must reassign equation numbrs
870 this->assign_eqn_numbers();
871
872 //Do a steady solve to map the nodes to the boundary of the ellipse
873 this->steady_newton_solve();
874
875 //Now unpin the rigid body...
876 this->unpin_rigid_body();
877
878 //...and then repin the position of the centre of mass
879 ImmersedRigidBodyElement* rigid_element1_pt =
880 dynamic_cast<ImmersedRigidBodyElement*>(Rigid_body_pt[0]);
881 rigid_element1_pt->pin_centre_of_mass_coordinate(0);
882 rigid_element1_pt->pin_centre_of_mass_coordinate(1);
883
884 //and then reassign equation numbers
885 this->assign_eqn_numbers();
886
887} //end_solve_for_consistent_nodal_positions
888
889
890
891
892//============start_of_create_lagrange_multiplier_elements===============
893/// Create elements that impose the prescribed boundary displacement
894/// for the pseudo-solid fluid mesh
895//=======================================================================
896template<class ELEMENT>
899{
900 // The idea is to apply Lagrange multipliers to the boundaries in
901 // the mesh that have associated geometric objects
902
903 //Find the number of boundaries
904 unsigned n_boundary = Fluid_mesh_pt->nboundary();
905
906 // Loop over the boundaries
907 for(unsigned b=0;b<n_boundary;b++)
908 {
909 //Get the geometric object associated with the boundary
910 GeomObject* boundary_geom_obj_pt =
911 Fluid_mesh_pt->boundary_geom_object_pt(b);
912
913 //Only bother to do anything if there is a geometric object
914 if(boundary_geom_obj_pt!=0)
915 {
916 // How many bulk fluid elements are adjacent to boundary b?
917 unsigned n_element = Fluid_mesh_pt->nboundary_element(b);
918
919 // Loop over the bulk fluid elements adjacent to boundary b?
920 for(unsigned e=0;e<n_element;e++)
921 {
922 // Get pointer to the bulk fluid element that is
923 // adjacent to boundary b
924 ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
925 Fluid_mesh_pt->boundary_element_pt(b,e));
926
927 //Find the index of the face of element e along boundary b
928 int face_index = Fluid_mesh_pt->face_index_at_boundary(b,e);
929
930 // Create new element. Note that we use different Lagrange
931 // multiplier fields for each distinct boundary (here indicated
932 // by b.
933 ImposeDisplacementByLagrangeMultiplierElement<ELEMENT>* el_pt =
934 new ImposeDisplacementByLagrangeMultiplierElement<ELEMENT>(
935 bulk_elem_pt,face_index,b);
936
937 // Add it to the mesh
938 Lagrange_multiplier_mesh_pt->add_element_pt(el_pt);
939
940 // Set the GeomObject that defines the boundary shape and set
941 // which bulk boundary we are attached to (needed to extract
942 // the boundary coordinate from the bulk nodes)
943 el_pt->set_boundary_shape_geom_object_pt(
944 boundary_geom_obj_pt,b);
945
946 // Loop over the nodes to pin Lagrange multiplier
947 unsigned nnod=el_pt->nnode();
948 for(unsigned j=0;j<nnod;j++)
949 {
950 Node* nod_pt = el_pt->node_pt(j);
951
952 // How many nodal values were used by the "bulk" element
953 // that originally created this node?
954 unsigned n_bulk_value=el_pt->nbulk_value(j);
955
956 // Pin two of the four Lagrange multipliers at vertices
957 // This is not totally robust, but will work in this application
958 unsigned nval=nod_pt->nvalue();
959 if (nval==7)
960 {
961 for (unsigned i=0;i<2;i++)
962 {
963 // Pin lagrangian multipliers
964 nod_pt->pin(n_bulk_value+2+i);
965 }
966 }
967 }
968 } // end loop over the element
969 } //End of case if there is a geometric object
970 } //End of loop over boundaries
971}
972// end of create_lagrange_multiplier_elements
973
974
975//===============start_delete_lagrange_multiplier_elements==================
976/// Delete elements that impose the prescribed boundary displacement
977/// and wipe the associated mesh
978//==========================================================================
979template<class ELEMENT>
982{
983 // How many surface elements are in the surface mesh
984 unsigned n_element = Lagrange_multiplier_mesh_pt->nelement();
985
986 // Loop over the surface elements
987 for(unsigned e=0;e<n_element;e++)
988 {
989 // Kill surface element
990 delete Lagrange_multiplier_mesh_pt->element_pt(e);
991 }
992
993 // Wipe the mesh
994 Lagrange_multiplier_mesh_pt->flush_element_and_node_storage();
995
996} // end of delete_lagrange_multiplier_elements
997
998
999
1000
1001//============start_of_create_drag_elements===============
1002/// Create elements that calculate the drag and torque on
1003/// the obstacles in the fluid mesh
1004//=======================================================================
1005template<class ELEMENT>
1007{
1008 //The idea is only to attach drag elements to those
1009 //boundaries associated with each particular rigid body
1010 //The loop is slightly inefficient, but should work in general
1011
1012 // Get the number of rigid bodies
1013 unsigned n_rigid = Rigid_body_pt.size();
1014
1015 // Get the number of boundaries in the mesh
1016 unsigned n_boundary = Fluid_mesh_pt->nboundary();
1017
1018 //Loop over the rigid bodies
1019 for(unsigned r=0;r<n_rigid;r++)
1020 {
1021 //Get the rigid boundary geometric object
1022 ImmersedRigidBodyElement* rigid_el_pt =
1023 dynamic_cast<ImmersedRigidBodyElement*>(this->Rigid_body_pt[r]);
1024
1025 // Loop over all boundaries
1026 for(unsigned b=0;b<n_boundary;b++)
1027 {
1028 //Does the boundary correspond to the current rigid body
1029 if(dynamic_cast<ImmersedRigidBodyElement*>
1030 (Fluid_mesh_pt->boundary_geom_object_pt(b)) == rigid_el_pt)
1031 {
1032 // How many bulk fluid elements are adjacent to boundary b?
1033 unsigned n_element = Fluid_mesh_pt->nboundary_element(b);
1034
1035 // Loop over the bulk fluid elements adjacent to boundary b?
1036 for(unsigned e=0;e<n_element;e++)
1037 {
1038 // Get pointer to the bulk fluid element that is
1039 // adjacent to boundary b
1040 ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
1041 Fluid_mesh_pt->boundary_element_pt(b,e));
1042
1043 //Find the index of the face of element e along boundary b
1044 int face_index = Fluid_mesh_pt->face_index_at_boundary(b,e);
1045
1046 // Create new element. Note that we use different Lagrange
1047 // multiplier fields for each distinct boundary (here indicated
1048 // by b.
1049 NavierStokesSurfaceDragTorqueElement<ELEMENT>* el_pt =
1050 new NavierStokesSurfaceDragTorqueElement<ELEMENT>(
1051 bulk_elem_pt,face_index);
1052
1053 // Add it to the mesh
1054 Drag_mesh_pt[r]->add_element_pt(el_pt);
1055
1056 //Set the original centre of rotation
1057 //from the rigid body in the drag element
1058 for(unsigned i=0;i<2;i++)
1059 {el_pt->centre_of_rotation(i) =
1060 rigid_el_pt->initial_centre_of_mass(i);}
1061
1062 //Set the data to the translation and rotation data
1063 //as well because these data will enter the Jacobian terms
1064 //in the DragTorqueElement
1065 el_pt->set_translation_and_rotation(rigid_el_pt->geom_data_pt(0));
1066 } // end loop over the element
1067 }
1068 } //End of loop over boundaries
1069 } // end loop over the rigid bodies
1070}
1071// end of create_drag_elements
1072
1073
1074//=======================================================================
1075/// Delete elements that calculate the drag and torque on the
1076/// boundaries
1077//=======================================================================
1078template<class ELEMENT>
1080{
1081 unsigned n_bodies = Drag_mesh_pt.size();
1082 for(unsigned n=0;n<n_bodies;n++)
1083 {
1084 // How many surface elements are in the surface mesh
1085 unsigned n_element = Drag_mesh_pt[n]->nelement();
1086
1087 // Loop over the surface elements
1088 for(unsigned e=0;e<n_element;e++)
1089 {
1090 // Kill surface element
1091 delete Drag_mesh_pt[n]->element_pt(e);
1092 }
1093
1094 // Wipe the mesh
1095 Drag_mesh_pt[n]->flush_element_and_node_storage();
1096 }
1097} // end of delete_drag_elements
1098
1099
1100
1101
1102//==start_of_doc_solution=================================================
1103/// Doc the solution
1104//========================================================================
1105template<class ELEMENT>
1107 const bool& project)
1108{
1109
1110 oomph_info << "Docing step: " << this->Doc_info.number()
1111 << std::endl;
1112
1113 ofstream some_file;
1114 char filename[100];
1115
1116 // Number of plot points
1117 unsigned npts;
1118 npts=5;
1119
1120
1121 // Output solution and projection files
1122 if(!project)
1123 {
1124 sprintf(filename,"%s/soln%i.dat",
1125 this->Doc_info.directory().c_str(),
1126 this->Doc_info.number());
1127 }
1128 else
1129 {
1130 sprintf(filename,"%s/proj%i.dat",
1131 this->Doc_info.directory().c_str(),
1132 this->Doc_info.number()-1);
1133 }
1134
1135 // Assemble square of L2 norm
1136 double square_of_l2_norm=0.0;
1137 unsigned nel=Fluid_mesh_pt->nelement();
1138 for (unsigned e=0;e<nel;e++)
1139 {
1140 square_of_l2_norm+=
1141 dynamic_cast<ELEMENT*>(this->Fluid_mesh_pt->element_pt(e))->
1142 square_of_l2_norm();
1143 }
1144
1145 std::cout << "Checking " << sqrt(square_of_l2_norm) << "\n";
1146 this->Norm_file << sqrt(square_of_l2_norm) << "\n";
1147
1148 some_file.open(filename);
1149 some_file << dynamic_cast<ELEMENT*>(this->Fluid_mesh_pt->element_pt(0))
1150 ->variable_identifier();
1151 this->Fluid_mesh_pt->output(some_file,npts);
1152 some_file.close();
1153
1154 // No trace file writing after projection
1155 if(project) return;
1156
1157 //Output the boundary only if not projecting
1158 sprintf(filename,"%s/int%i.dat",
1159 this->Doc_info.directory().c_str(),
1160 this->Doc_info.number());
1161 some_file.open(filename);
1162 this->Lagrange_multiplier_mesh_pt->output(some_file);
1163 some_file.close();
1164
1165 // Increment the doc_info number
1166 this->Doc_info.number()++;
1167
1168 //Output the motion of the centre of gravity
1169 dynamic_cast<ImmersedRigidBodyElement*>(this->Rigid_body_pt[0])->
1170 output_centre_of_gravity(this->Cog_file);
1171
1172 //Output the exact solution
1173 this->output_exact_solution(this->Cog_exact_file);
1174
1175}
1176
1177
1178//=====================================================================
1179/// Output the exact solution
1180//=====================================================================
1181template<class ELEMENT>
1183output_exact_solution(std::ofstream &output_file)
1184 {
1185 //Get the current time
1186 double time = this->time();
1187 //Output the exact solution
1188 output_file << time << " " << Jeffery_Solution::angle(time)
1189 << " " << Jeffery_Solution::velocity(time)
1190 << " " << Jeffery_Solution::acceleration(time) << std::endl;
1191 }
1192
1193
1194
1195//==========start_of_main======================================
1196/// Driver code for immersed ellipse problem
1197//============================================================
1198int main(int argc, char **argv)
1199{
1200 // Store command line arguments
1201 CommandLineArgs::setup(argc,argv);
1202
1203 // Define possible command line arguments and parse the ones that
1204 // were actually specified
1205
1206 // Validation?
1207 CommandLineArgs::specify_command_line_flag("--validation");
1208
1209 // Parse command line
1210 CommandLineArgs::parse_and_assign();
1211
1212 // Doc what has actually been specified on the command line
1213 CommandLineArgs::doc_specified_flags();
1214
1215 // Create problem in initial configuration
1217 ProjectableTaylorHoodElement<MyTaylorHoodElement> > problem;
1218
1219 //Initially ensure that the nodal positions are consistent with
1220 //their weak imposition
1222
1223 // Initialise timestepper
1224 double dt=0.05;
1225 problem.initialise_dt(dt);
1226
1227 // Perform impulsive start
1228 problem.assign_initial_values_impulsive();
1229
1230 // Output initial conditions
1231 problem.doc_solution();
1232
1233 // Solve problem a few times on given mesh
1234 unsigned nstep=3;
1235 for (unsigned i=0;i<nstep;i++)
1236 {
1237 // Solve the problem
1238 problem.unsteady_newton_solve(dt);
1239 problem.doc_solution();
1240 }
1241
1242 // Now do a couple of adaptations
1243 unsigned ncycle=200;
1244 if (CommandLineArgs::command_line_flag_has_been_set("--validation"))
1245 {
1246 ncycle=1;
1247 oomph_info << "Only doing one cycle during validation\n";
1248 }
1249
1250 for (unsigned j=0;j<ncycle;j++)
1251 {
1252 // Adapt the problem
1253 problem.adapt();
1254
1255 //Solve problem a few times
1256 for (unsigned i=0;i<nstep;i++)
1257 {
1258 // Solve the problem
1259 problem.unsteady_newton_solve(dt);
1260 problem.doc_solution();
1261 }
1262 }
1263
1264} //end of main
//////////////////////////////////////////////////////// ////////////////////////////////////////////...
void position(const unsigned &t, const Vector< double > &xi, Vector< double > &r) const
~GeneralEllipse()
Empty Destructor.
GeneralEllipse(const double &centre_x, const double &centre_y, const double &a, const double &b)
Simple Constructor that transfers appropriate geometric parameters into internal data.
void position(const Vector< double > &xi, Vector< double > &r) const
Return the position of the ellipse boundary as a function of the angle xi[0].
//////////////////////////////////////////////////////// ////////////////////////////////////////////...
void pin_rigid_body()
Pin the degrees of freedom associated with the solid bodies.
void doc_solution(const bool &project=false)
Doc the solution.
RefineableSolidTriangleMesh< ELEMENT > * Fluid_mesh_pt
Pointer to Fluid_mesh.
void create_drag_elements()
Create elements that calculate the drag and torque on the boundaries.
DocInfo Doc_info
Internal DocInfo object.
void output_exact_solution(std::ofstream &output_file)
Output the exact solution.
void solve_for_consistent_nodal_positions()
Function that solves a simplified problem to ensure that the positions of the boundary nodes are init...
Mesh * Rigid_body_mesh_pt
Mesh of the generalised elements for the rigid bodies.
void actions_before_implicit_timestep()
Reset the boundary conditions when timestepping.
void create_lagrange_multiplier_elements()
Create elements that enforce prescribed boundary motion for the pseudo-solid fluid mesh by Lagrange m...
ofstream Cog_exact_file
File to document the exact motion of the centre of gravity.
void complete_problem_setup()
Set boundary condition, assign auxiliary node update fct. Complete the build of all elements,...
void unpin_rigid_body()
Unpin the degrees of freedom associated with the solid bodies.
void delete_lagrange_multiplier_elements()
Delete elements that impose the prescribed boundary displacement and wipe the associated mesh.
UnstructuredImmersedEllipseProblem()
Constructor.
ofstream Norm_file
File to document the norm of the solution (for validation purposes)
void delete_drag_elements()
Delete elements that calculate the drag and torque on the boundaries.
Vector< Mesh * > Drag_mesh_pt
Mesh of drag elements.
void actions_before_adapt()
Wipe the meshes of Lagrange multiplier and drag elements.
Vector< GeomObject * > Rigid_body_pt
Storage for the geom object.
void actions_after_adapt()
Rebuild the meshes of Lagrange multiplier and drag elements.
void set_boundary_velocity()
Set the boundary velocity.
ofstream Cog_file
File to document the motion of the centre of gravity.
SolidMesh * Lagrange_multiplier_mesh_pt
Pointers to mesh of Lagrange multiplier elements.
TriangleMeshPolygon * Outer_boundary_polygon_pt
Triangle mesh polygon for outer boundary.
void actions_before_newton_convergence_check()
Re-apply the no slip condition (imposed indirectly via dependent velocities)
int main(int argc, char **argv)
Driver code for immersed ellipse problem.
Exact solution for the rotation of an ellipse in unbounded shear flow as computed by Jeffery (1922)
double velocity(const double &t)
Angular velocity as function of time t.
double angle(const double &t)
Angular position as a function of time t.
double acceleration(const double &t)
Angular acceleration as a function of time t (should always be zero)
Namespace for Problem Parameters.
double Lambda_sq
Pseudo-solid (mesh) "density" Set to zero because we don't want inertia in the node update!
double A
Initial axis of the elliptical solid in x-direction.
double Density_ratio
Density ratio (Solid density / Fluid density)
double St
Strouhal number.
ConstitutiveLaw * Constitutive_law_pt
Constitutive law used to determine the mesh deformation.
double Nu
Pseudo-solid (mesh) Poisson ratio.
double B
Initial axis of the elliptical solid in y-direction (N.B. 2B = 1 is the reference length scale)
double Re
Reynolds number.