3d_rayleigh_instability_surfactant.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// Version 1.0; svn revision $LastChangedRevision: 1207 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-06-24 15:10:32 +0100 (Fri, 24 Jun 2016) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //Three-dimensional free-surface test case including insoluble
31 //surfactant transport. This is a 3D implementation of the
32 //axisymmetric Rayleigh--Plateau problem solved in
33 //rayleigh_instability_insoluble_surfactant.cc, which
34 //should be a hard test
35 //because it's implemented in Cartesian coordinates.
36 //The main aim is to test the implementation of the surface transport
37 //equations in 3D (2D surface).
38 
39 // The oomphlib headers
40 #include "generic.h"
41 #include "navier_stokes.h"
42 #include "fluid_interface.h"
43 
44 // The basic mesh
45 #include "meshes/single_layer_cubic_spine_mesh.h"
46 
47 using namespace std;
48 using namespace oomph;
49 
50 
51 //==start_of_namespace==================================================
52 /// Namepspace for global parameters, chosen from Campana et al. as
53 /// in the axisymmetric problem.
54 //======================================================================
55 namespace Global_Physical_Variables
56 {
57  //Film thickness parameter
58  double Film_Thickness = 0.2;
59 
60  /// Reynolds number
61  double Re = 40.0;
62 
63  /// Womersley = Reynolds times Strouhal
64  double ReSt = Re;
65 
66  /// Product of Reynolds and Froude number
67  double ReInvFr = 0.0;
68 
69  /// Capillary number
70  double Ca = pow(Film_Thickness,3.0);
71 
72  /// Direction of gravity
73  Vector<double> G(3);
74 
75 // Wavelength of the domain
76  double Alpha = 1.047;
77 
78  /// Free surface cosine deformation parameter
79  double Epsilon = 1.0e-3;
80 
81  /// \short Surface Elasticity number (weak case)
82  double Beta = 3.6e-3;
83 
84  /// \short Surface Peclet number
85  double Peclet_S = 4032.0;
86 
87  /// \short Sufrace Peclet number multiplied by Strouhal number
88  double Peclet_St_S = 1.0;
89 
90 } // End of namespace
91 
92 
93 
94 //=======================================================================
95 /// Function-type-object to perform comparison of elements in y-direction
96 //=======================================================================
98 {
99 public:
100 
101  /// Comparison. Are the values identical or not?
102  bool operator()(GeneralisedElement* const &x, GeneralisedElement* const &y)
103  const
104  {
105  FiniteElement* cast_x = dynamic_cast<FiniteElement*>(x);
106  FiniteElement* cast_y = dynamic_cast<FiniteElement*>(y);
107 
108  if((cast_x ==0) || (cast_y==0)) {return 0;}
109  else
110  {return (cast_x->node_pt(0)->x(1)< cast_y->node_pt(0)->x(1));}
111 
112  }
113 };
114 
115 
116 namespace oomph
117 {
118 //===================================================================
119 /// Deform the existing cubic spine mesh into a annular section
120 /// with spines directed radially inwards from the wall
121 //===================================================================
122 
123 template<class ELEMENT>
124 class AnnularSpineMesh : public SingleLayerCubicSpineMesh<ELEMENT>
125 {
126 
127 public:
128  AnnularSpineMesh(const unsigned &n_r, const unsigned &n_y,
129  const unsigned &n_theta,
130  const double &r_min, const double &r_max,
131  const double &l_y,
132  const double &theta_min, const double &theta_max,
133  TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper):
134  //This will make a cubic mesh with n_r in the x-direction
135  //n_y in the y-direction and n_theta in the z-direction
136  //The coordinates will run from 0 to r_max, 0 to l_y and 0 to theta_max
137  SingleLayerCubicSpineMesh<ELEMENT>(n_theta,n_y,n_r,r_max,l_y,theta_max,
138  time_stepper_pt)
139  {
140 
141  //Find out how many nodes there are
142  unsigned n_node = this->nnode();
143 
144  //Loop over all nodes
145  for (unsigned n=0;n<n_node;n++)
146  {
147  //pointer to node
148  Node* nod_pt=this->node_pt(n);
149  SpineNode* spine_node_pt=dynamic_cast<SpineNode*>(nod_pt);
150  //Get x/y/z coordinates
151  double x_old=nod_pt->x(0);
152  double y_old=nod_pt->x(1);
153  double z_old=nod_pt->x(2);
154 
155  //Mapping
156 
157  double r = r_min + (r_max-r_min)*z_old;
158  double theta = (theta_min + (theta_max-theta_min)*x_old);
159  double y = y_old;
160 
161 
162  if(spine_node_pt->spine_pt()->ngeom_parameter()==0)
163  {
164  spine_node_pt->h() =
166  spine_node_pt->spine_pt()->add_geom_parameter(theta);
167  }
168 
169  //cout << spine_node_pt->spine_pt()->ngeom_parameter() << std::endl;
170  //Set new nodal coordinates
171  nod_pt->x(0)=r*cos(theta);
172  nod_pt->x(2)=r*sin(theta);
173  nod_pt->x(1)=y;
174  //cout << "one" << theta << std::endl;
175  }
176 
177  }
178 
179 
180  virtual void spine_node_update(SpineNode* spine_node_pt)
181  {
182  //Get fraction along the spine
183  double W = spine_node_pt->fraction();
184  //Get theta along the spine
185  double theta = spine_node_pt->spine_pt()->geom_parameter(0);
186  //Get spine height
187  double H = spine_node_pt->h();
188  //Set the value of y
189  spine_node_pt->x(0) = (1.0-W*H)*cos(theta);
190  spine_node_pt->x(2) = (1.0-W*H)*sin(theta);
191  }
192 };
193 }
194 
195 //======================================================================
196 /// Single fluid interface problem including transport of an
197 /// insoluble surfactant.
198 //======================================================================
199 template<class ELEMENT, class TIMESTEPPER>
200 class InterfaceProblem : public Problem
201 {
202 
203 public:
204 
205  /// Constructor: Pass number of elements in x and y directions. Also lengths
206  /// of the domain in x- and y-directions and the height of the layer
207 
208  InterfaceProblem(const unsigned &n_r, const unsigned &n_y,
209  const unsigned &n_theta, const double &r_min,
210  const double &r_max, const double &l_y,
211  const double &theta_max);
212 
213  /// Spine heights/lengths are unknowns in the problem so their
214  /// values get corrected during each Newton step. However,
215  /// changing their value does not automatically change the
216  /// nodal positions, so we need to update all of them
217  void actions_before_newton_convergence_check(){Bulk_mesh_pt->node_update();}
218 
219  /// Run an unsteady simulation with specified number of steps
220  void unsteady_run(const unsigned& nstep);
221 
222  /// Doc the solution
223  void doc_solution(DocInfo& doc_info);
224 
225  /// Compute the total mass
227  {
228  double mass = 0.0;
229 
230  // Determine number of 1D interface elements in mesh
231  const unsigned n_interface_element = Surface_mesh_pt->nelement();
232 
233  // Loop over the interface elements
234  for(unsigned e=0;e<n_interface_element;e++)
235  {
236  // Upcast from GeneralisedElement to the present element
239  (Surface_mesh_pt->element_pt(e));
240 
241  mass += el_pt->integrate_c();
242  }
243  return mass;
244  }
245 
246 
247 
248  private:
249 
250  /// Trace file
251  ofstream Trace_file;
252 
253  /// Axial lengths of domain
254  double R_max;
255 
256  double L_y;
257 
258  /// Pointer to bulk mesh
260 
261  /// Pointer to the surface mes
263 
264  /// Pointer to a node for documentation purposes
266 
267 };
268 
269 
270 //====================================================================
271 /// Problem constructor
272 //====================================================================
273 template<class ELEMENT, class TIMESTEPPER>
275 (const unsigned &n_r, const unsigned &n_y,const unsigned &n_theta,
276  const double &r_min,
277  const double &r_max, const double &l_y, const double &theta_max)
278  : R_max(r_max), L_y(l_y)
279 {
280  //this->linear_solver_pt() = new HSL_MA42;
281 
282  //static_cast<HSL_MA42*>(this->linear_solver_pt())->lenbuf_factor0() = 3.0;
283  //static_cast<HSL_MA42*>(this->linear_solver_pt())->lenbuf_factor1() = 3.0;
284  //static_cast<HSL_MA42*>(this->linear_solver_pt())->lenbuf_factor2() = 3.0;
285 
286  //Allocate the timestepper
287  add_time_stepper_pt(new TIMESTEPPER);
288 
289  //Now create the bulk mesh -- this should be your new annular mesh
290  Bulk_mesh_pt = new AnnularSpineMesh<ELEMENT>
291  (n_r,n_y,n_theta,r_min,r_max,l_y,0.0,theta_max,time_stepper_pt());
292 
293  //Set the documented node
294  Document_node_pt = Bulk_mesh_pt->boundary_node_pt(5,0);
295 
296 
297  //Create the surface mesh that will contain the interface elements
298  //First create storage, but with no elements or nodes
299  Surface_mesh_pt = new Mesh;
300 
301  // Loop over those elements adjacent to the free surface,
302  // which we shall choose to be the upper surface
303  for(unsigned e1=0;e1<n_y;e1++)
304  {
305  for(unsigned e2=0;e2<n_theta;e2++)
306  {
307  // Set a pointer to the bulk element we wish to our interface
308  // element to
309  FiniteElement* bulk_element_pt =
310  Bulk_mesh_pt->finite_element_pt(n_theta*n_y*(n_r-1) + e2 + e1*n_theta);
311 
312  // Create the interface element (on face 3 of the bulk element)
313  FiniteElement* interface_element_pt =
315 
316  // Add the interface element to the surface mesh
317  this->Surface_mesh_pt->add_element_pt(interface_element_pt);
318  }
319  }
320 
321  // Add the two sub-meshes to the problem
322  add_sub_mesh(Bulk_mesh_pt);
323  add_sub_mesh(Surface_mesh_pt);
324 
325  // Combine all sub-meshes into a single mesh
326  build_global_mesh();
327 
328  //Pin all nodes on the bottom
329  unsigned long n_boundary_node = Bulk_mesh_pt->nboundary_node(0);
330  for(unsigned long n=0;n<n_boundary_node;n++)
331  {
332  for(unsigned i=0;i<3;i++)
333  {
334  Bulk_mesh_pt->boundary_node_pt(0,n)->pin(i);
335  }
336  }
337 
338  //On the front and back (y=const) pin in y-direction
339  for(unsigned b=1;b<4;b+=2)
340  {
341  n_boundary_node = Bulk_mesh_pt->nboundary_node(b);
342  for(unsigned long n=0;n<n_boundary_node;n++)
343  {
344  Bulk_mesh_pt->boundary_node_pt(b,n)->pin(1);
345  }
346  }
347 
348  //On sides pin in z-direction
349  for(unsigned b=4;b<5;b+=2)
350  {
351  n_boundary_node = Bulk_mesh_pt->nboundary_node(b);
352  for(unsigned long n=0;n<n_boundary_node;n++)
353  {
354  Bulk_mesh_pt->boundary_node_pt(b,n)->pin(2);
355  }
356  }
357 
358  // pinning the top surface
359  for(unsigned b=2;b<3;b++)
360  {
361  n_boundary_node = Bulk_mesh_pt->nboundary_node(b);
362  for(unsigned long n=0;n<n_boundary_node;n++)
363  {
364  Bulk_mesh_pt->boundary_node_pt(b,n)->pin(0);
365  }
366  }
367 
368  //Loop over the upper surface and set the surface concentration
369  {
370  unsigned b=5;
371  n_boundary_node = Bulk_mesh_pt->nboundary_node(b);
372  for(unsigned long n=0;n<n_boundary_node;n++)
373  {
374  Bulk_mesh_pt->boundary_node_pt(b,n)->set_value(3,1.0);;
375  }
376  }
377 
378 
379  //Create a Data object whose single value stores the
380  //external pressure
381  Data* external_pressure_data_pt = new Data(1);
382 
383  // Set and pin the external pressure to some random value
384  external_pressure_data_pt->set_value(0,1.31);
385  external_pressure_data_pt->pin(0);
386 
387  //Complete the problem setup to make the elements fully functional
388 
389  //Loop over the elements in the layer
390  unsigned n_bulk=Bulk_mesh_pt->nelement();
391  for(unsigned e=0;e<n_bulk;e++)
392  {
393  //Cast to a fluid element
394  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
395 
396  //Set the Reynolds number, etc
397  el_pt->re_pt() = &Global_Physical_Variables::Re;
398  el_pt->re_st_pt() = &Global_Physical_Variables::ReSt;
399  el_pt->re_invfr_pt() = &Global_Physical_Variables::ReInvFr;
400  el_pt->g_pt() = &Global_Physical_Variables::G;
401  }
402 
403  //Loop over 2D interface elements and set capillary number and
404  //the external pressure
405  unsigned long interface_element_pt_range =
406  Surface_mesh_pt->nelement();
407  for(unsigned e=0;e<interface_element_pt_range;e++)
408  {
409  //Cast to a interface element
412  (Surface_mesh_pt->element_pt(e));
413 
414  //Set the Capillary number
416 
417  //Pass the Data item that contains the single external pressure value
418  el_pt->set_external_pressure_data(external_pressure_data_pt);
419 
420  // Set the surface elasticity number
422 
423  // Set the surface peclect number
425 
426  // Set the surface peclect number multiplied by strouhal number
428 
429  }
430 
431  //Do equation numbering
432  cout << assign_eqn_numbers() << std::endl;
433 
434  //Now sort the elements to minimise frontwidth
435  std::sort(mesh_pt()->element_pt().begin(),
436  mesh_pt()->element_pt().end(),
437  ElementCmp());
438 
439 }
440 
441 
442 
443 //========================================================================
444 /// Doc the solution
445 //========================================================================
446 template<class ELEMENT, class TIMESTEPPER>
448 {
449 
450  ofstream some_file;
451  char filename[100];
452 
453  // Number of plot points
454  unsigned npts=5;
455 
456  //Output the time
457  cout << "Time is now " << time_pt()->time() << std::endl;
458 
459  // Doc
460  Trace_file << time_pt()->time();
461  Trace_file << " " << Document_node_pt->x(0) << " "
462  << this->compute_total_mass()
463  << std::endl;
464 
465 
466  // Output solution, bulk elements followed by surface elements
467  /*sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
468  doc_info.number());
469  some_file.open(filename);
470  Bulk_mesh_pt->output(some_file,npts);
471  some_file.close();*/
472  //Change the file name and output surface in separate file
473  sprintf(filename,"%s/surface%i.dat",doc_info.directory().c_str(),
474  doc_info.number());
475  some_file.open(filename);
476  Surface_mesh_pt->output(some_file,npts);
477  some_file.close();
478 
479 }
480 
481 
482 
483 
484  //=============================================================================
485 /// Unsteady run with specified number of steps
486 //=============================================================================
487 template<class ELEMENT, class TIMESTEPPER>
489 {
490 
491  // Increase maximum residual
492  Problem::Max_residuals=500.0;
493 
494  //Distort the mesh
495  const double epsilon=Global_Physical_Variables::Epsilon;
496  unsigned Nspine = Bulk_mesh_pt->nspine();
497  for(unsigned i=0;i<Nspine;i++)
498  {
499  double y_value = Bulk_mesh_pt->boundary_node_pt(0,i)->x(1);
500 
501  Bulk_mesh_pt->spine_pt(i)->height() =
503  + epsilon*cos(Global_Physical_Variables::Alpha*y_value));
504  }
505 
506  //Make sure to update
507  Bulk_mesh_pt->node_update();
508 
509  // Doc info object
510  DocInfo doc_info;
511 
512  // Set output directory
513  doc_info.set_directory("RESLT");
514  doc_info.number()=0;
515 
516  // Open trace file
517  char filename[100];
518  sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
519  Trace_file.open(filename);
520 
521  Trace_file << "VARIABLES=\"time\",";
522  Trace_file << "\"h<sub>left</sub>\",\"h<sub>right</sub>\"";
523  Trace_file << "\n";
524 
525  //Set value of dt
526  double dt = 0.1;
527 
528  //Initialise all time values
529  assign_initial_values_impulsive(dt);
530 
531  //Doc initial solution
532  doc_solution(doc_info);
533 
534 //Loop over the timesteps
535  for(unsigned t=1;t<=nstep;t++)
536  {
537  cout << "TIMESTEP " << t << std::endl;
538 
539  //Take one fixed timestep
540  unsteady_newton_solve(dt);
541 
542  // Doc solution
543  doc_info.number()++;
544  doc_solution(doc_info);
545  }
546 
547 }
548 
549 
550 //==start_of_main========================================================
551 /// Driver code for unsteady two-layer fluid problem. If there are
552 /// any command line arguments, we regard this as a validation run
553 /// and perform only two steps.
554 
555 // In my version we will change nsteps in the programs
556 //======================================================================
557 int main(int argc, char *argv[])
558 {
559 
560  // Set physical parameters:
561 
562  //Set direction of gravity: Vertically downwards
566 
567  //Set the film thickness
569 
570  //Axial lngth of domain
571  double r_max = 1.0;
572  double r_min = r_max - h;
573 
574  const double pi = MathematicalConstants::Pi;
575 
576  double alpha = Global_Physical_Variables::Alpha;
577 
578  double l_y = pi/alpha;
579 
580  double theta_max = 0.5*pi;
581 
582  // Number of elements in r and theta direction
583  unsigned n_r = 4;
584 
585  unsigned n_theta = 4;
586 
587  // Number of elements in axial direction
588  unsigned n_y = 20;
589 
590 
591  {
592  //Set up the spine test problem
593  //The minimum values are all assumed to be zero.
595  problem(n_r,n_y,n_theta,r_min,r_max,l_y,theta_max);
596 
597  // Number of steps:
598  unsigned nstep;
599  if(argc > 1)
600  {
601  // Validation run: Just five steps
602  nstep=5;
603  }
604  else
605  {
606  // Full run otherwise
607  nstep=1000;
608  }
609 
610  // Run the unsteady simulation
611  problem.unsteady_run(nstep);
612  }
613 } // End of main
614 
615 
bool operator()(GeneralisedElement *const &x, GeneralisedElement *const &y) const
Comparison. Are the values identical or not?
double Peclet_S
Surface Peclet number.
double Peclet_St_S
Sufrace Peclet number multiplied by Strouhal number.
virtual void spine_node_update(SpineNode *spine_node_pt)
InterfaceProblem(const unsigned &n_r, const unsigned &n_y, const unsigned &n_theta, const double &r_min, const double &r_max, const double &l_y, const double &theta_max)
Problem constructor.
void set_external_pressure_data(Data *external_pressure_data_pt)
Set the Data that contains the single pressure value that specifies the "external pressure" for the i...
double *& peclet_s_pt()
Access function for pointer to the surface Peclet number.
double ReSt
Womersley = Reynolds times Strouhal.
double integrate_c()
Compute the concentration intergated over the surface area.
AnnularSpineMesh< ELEMENT > * Bulk_mesh_pt
Pointer to bulk mesh.
double Epsilon
Free surface cosine deformation parameter.
double R_max
Axial lengths of domain.
double *& beta_pt()
Access function for pointer to the Elasticity number.
Mesh * Surface_mesh_pt
Pointer to the surface mes.
double ReInvFr
Product of Reynolds and Froude number.
int main(int argc, char *argv[])
Vector< double > G(3)
Direction of gravity.
Node * Document_node_pt
Pointer to a node for documentation purposes.
void unsteady_run(const unsigned &nstep)
Run an unsteady simulation with specified number of steps.
double Beta
Surface Elasticity number (weak case)
AnnularSpineMesh(const unsigned &n_r, const unsigned &n_y, const unsigned &n_theta, const double &r_min, const double &r_max, const double &l_y, const double &theta_min, const double &theta_max, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
void doc_solution(DocInfo &doc_info)
Doc the solution.
double *& peclet_strouhal_s_pt()
Access function for pointer to the surface Peclet x Strouhal number.
Function-type-object to perform comparison of elements in y-direction.
double Alpha
Wavelength of the domain.
double compute_total_mass()
Compute the total mass.
double *& ca_pt()
Pointer to the Capillary number.