rayleigh_instability_soluble_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: 1104 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-01-09 08:06:50 +0000 (Sat, 09 Jan 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 // Driver for axisymmetric single-layer fluid problem. Plateau-Rayleigh
31 // instability (unstable if H>2*pi*R -> forming drops)
32 // in the presence of a soluble surfactant. The problem is described in
33 // Numerical analaysis of the Rayleigh instability in capillary tubes:
34 // The influence of surfactant solubility by
35 // D. M. Campana & F. A. Saita, Phys. Fluids.
36 // vol 18, 022104 (2006)
37 // The parameters are taken from their Table 3
38 
39 // Generic oomph-lib header
40 #include "generic.h"
41 
42 // Axisymmetric Navier-Stokes headers
43 #include "axisym_navier_stokes.h"
44 
45 // Interface headers
46 #include "fluid_interface.h"
47 //Header for the pesudo-solid mesh updates
48 #include "solid.h"
49 #include "constitutive.h"
50 
51 // The standard rectangular mesh
52 #include "meshes/rectangular_quadmesh.h"
53 
54 // Equations for bulk surfactant transport
56 
57 //Use the oomph and std namespaces
58 using namespace oomph;
59 using namespace std;
60 
61 //==start_of_namespace===================================================
62 /// Namespace for physical parameters
63 /// The parameter values are chosen to be those used in Figures 7, 12
64 /// in Campana & Saita
65 //=======================================================================
66 namespace Global_Physical_Variables
67 {
68 
69  //Film thickness parameter
70  double Film_Thickness = 0.18;
71 
72  /// Reynolds number
73  double Re = 1.0;
74 
75  /// Womersley number
76  double ReSt = Re; // (St = 1)
77 
78  /// Product of Reynolds number and inverse of Froude number
79  double ReInvFr = 0.0; // (Fr = 0)
80 
81  /// Capillary number
82  double Ca = pow(Film_Thickness,3.0);
83 
84  /// External pressure
85  double P_ext = 0.0;
86 
87  /// Direction of gravity
88  Vector<double> G(3);
89 
90  /// Wavelength of the domain
91  double Alpha = 0.8537;
92 
93  /// Free surface cosine deformation parameter
94  double Epsilon = 1.0e-3;
95 
96  /// \short Surface Elasticity number
97  double Beta = 0.1;
98 
99  /// \short Alpha for absorption kinetics
100  double Alpha_absorption = 1.0;
101 
102  /// \short K parameter that
103  /// describes solubility number
104  double K = 0.01;
105 
106  // Bulk Peclet number
107  double Pe = 10000.0;
108  // Bulk Peclet Strouhal number
109  double PeSt = Pe;
110 
111  //We shall adopt the same scaling
112  //as Campana & Saita, so divide
113  //the bulk equation through by
114  //Peclet number to give 1 as the
115  //nominal peclet number
116  double Pe_reference_scale=1.0;
117 
118  //Chose our diffusion in the bulk equations
119  //to be 1/Pe, after dividing through by Pe
120  double Diff = 1.0/Pe;
121 
122  /// \short Surface Peclet number
123  double Peclet_S = 10000.0;
124 
125  /// \short Sufrace Peclet number multiplied by Strouhal number
126  double Peclet_St_S = 1.0;
127 
128  /// \short Pvd file -- a wrapper for all the different
129  /// vtu output files plus information about continuous time
130  /// to facilitate animations in paraview
131  ofstream Pvd_file;
132 
133  /// Pseudo-solid Poisson ratio
134  double Nu = 0.1;
135 
136 } // End of namespace
137 
138 
139 
140 //================================================================
141 ///Interface class to handle the mass transport between bulk
142 ///and surface as well as the surfactant transport along the
143 //interface
144 //================================================================
145 template<class ELEMENT>
148 {
149 private:
150  ///Pointer to adsorption number
151  double *Alpha_pt;
152  ///Pointer to solubility number
153  double *K_pt;
154 
155  ///Storage for the index at
156  ///which the bulk concentration is stored
157  unsigned C_bulk_index;
158 
159 protected:
160 
161  ///Get the bulk surfactant concentration
162  double interpolated_C_bulk(const Vector<double> &s)
163  {
164  //Find number of nodes
165  unsigned n_node = this->nnode();
166 
167  //Get the nodal index at which the unknown is stored
168  const unsigned c_index = C_bulk_index;
169 
170  //Local shape function
171  Shape psi(n_node);
172 
173  //Find values of shape function
174  this->shape(s,psi);
175 
176  //Initialise value of C
177  double C = 0.0;
178 
179  //Loop over the local nodes and sum
180  for(unsigned l=0;l<n_node;l++)
181  {
182  C += this->nodal_value(l,c_index)*psi(l);
183  }
184 
185  return(C);
186  }
187 
188  /// The time derivative of the bulk surface concentration
189  double dc_bulk_dt_surface(const unsigned &l) const
190  {
191  // Get the data's timestepper
192  TimeStepper* time_stepper_pt= this->node_pt(l)->time_stepper_pt();
193 
194  //Initialise dudt
195  double dcdt=0.0;
196  //Loop over the timesteps, if there is a non Steady timestepper
197  if (time_stepper_pt->type()!="Steady")
198  {
199  //Find the index at which the variable is stored
200  const unsigned c_index = C_bulk_index;
201 
202  // Number of timsteps (past & present)
203  const unsigned n_time = time_stepper_pt->ntstorage();
204 
205  for(unsigned t=0;t<n_time;t++)
206  {
207  dcdt += time_stepper_pt->weight(1,t)*this->nodal_value(t,l,c_index);
208  }
209  }
210  return dcdt;
211  }
212 
213  /// \short Overload the Helper function to calculate the residuals and
214  /// jacobian entries. This particular function ensures that the
215  /// additional entries are calculated inside the integration loop
217  Vector<double> &residuals, DenseMatrix<double> &jacobian,
218  const unsigned &flag,const Shape &psif, const DShape &dpsifds,
219  const DShape &dpsifdS, const DShape &dpsifdS_div,
220  const Vector<double> &s,
221  const Vector<double> &interpolated_x, const Vector<double> &interpolated_n,
222  const double &W,const double &J)
223  {
224  //Call the standard transport equations
227  residuals, jacobian, flag,psif,dpsifds,dpsifdS,dpsifdS_div,
228  s, interpolated_x, interpolated_n, W, J);
229 
230  //Add the additional mass transfer terms
231  const double k = this->k();
232  const double alpha = this->alpha();
233  //Find out how many nodes there are
234  unsigned n_node = this->nnode();
235 
236  //Storage for the local equation numbers and unknowns
237  int local_eqn = 0, local_unknown = 0;
238 
239  //Surface advection-diffusion equation
240 
241  //Find the index at which the bulk concentration is stored
242  unsigned c_bulk_index = this->C_bulk_index;
243 
244  //Now calculate the bulk and surface concentrations at this point
245  //Assuming the same shape functions are used (which they are)
246  double interpolated_C = 0.0;
247  double interpolated_C_bulk = 0.0;
248 
249  //Loop over the shape functions
250  for(unsigned l=0;l<n_node;l++)
251  {
252  const double psi = psif(l);
253  const double C_ = this->nodal_value(l,this->C_index[l]);
254 
255  interpolated_C += C_*psi;
256  interpolated_C_bulk += this->nodal_value(l,c_bulk_index)*psi;
257  }
258 
259  //Pre compute the flux between the surface and bulk
260  double flux = alpha*(interpolated_C_bulk - interpolated_C);
261 
262  //Now we add the flux term to the appropriate residuals
263  for(unsigned l=0;l<n_node;l++)
264  {
265  //BULK FLUX
266 
267  //Read out the appropriate local equation
268  local_eqn = this->nodal_local_eqn(l,c_bulk_index);
269 
270  //If not a boundary condition
271  if(local_eqn >= 0)
272  {
273  //Add flux to the bulk
274  residuals[local_eqn] -= k*flux*psif(l)*W*J;
275 
276  if(flag)
277  {
278  for(unsigned l2=0;l2<n_node;l2++)
279  {
280  local_unknown = this->nodal_local_eqn(l2,this->C_index[l2]);
281  if(local_unknown >= 0)
282  {
283  jacobian(local_eqn,local_unknown) += k*alpha*psif(l2)*psif(l)*W*J;
284  }
285 
286  local_unknown = this->nodal_local_eqn(l2,c_bulk_index);
287  if(local_unknown >= 0)
288  {
289  jacobian(local_eqn,local_unknown) -= k*alpha*psif(l2)*psif(l)*W*J;
290  }
291  }
292  }
293  } //End of contribution to bulk
294 
295 
296  //SURFACE FLUX
297  //Read out the appropriate local equation
298  local_eqn = this->nodal_local_eqn(l,this->C_index[l]);
299 
300  //If not a boundary condition
301  if(local_eqn >= 0)
302  {
303  //Add flux to the surface
304  residuals[local_eqn] -= flux*psif(l)*W*J;
305 
306  if(flag)
307  {
308  for(unsigned l2=0;l2<n_node;l2++)
309  {
310  local_unknown = this->nodal_local_eqn(l2,this->C_index[l2]);
311  if(local_unknown >= 0)
312  {
313  jacobian(local_eqn,local_unknown) += alpha*psif(l2)*psif(l)*W*J;
314  }
315 
316  local_unknown = this->nodal_local_eqn(l2,c_bulk_index);
317  if(local_unknown >= 0)
318  {
319  jacobian(local_eqn,local_unknown) -= alpha*psif(l2)*psif(l)*W*J;
320  }
321  }
322  } //End of contribution to surface
323 
324  }
325  }
326 
327  }
328 
329 public:
330  ///Constructor that passes the bulk element and face index down
331  ///to the underlying
333  FiniteElement* const &element_pt, const int &face_index) :
335  (element_pt,face_index)
336  {
337  //Initialise the values
338  Alpha_pt = &this->Default_Physical_Constant_Value;
339  K_pt = &this->Default_Physical_Constant_Value;
340 
341  //Hack because I know the bulk index in this case
342  C_bulk_index = 3;
343  }
344 
345 
346  ///Return the adsorption number
347  double alpha() {return *Alpha_pt;}
348 
349  ///Return the solubility nubmer
350  double k() {return *K_pt;}
351 
352  ///Access function for pointer to adsorption number
353  double* &alpha_pt() {return Alpha_pt;}
354 
355  ///Access function for pointer to solubility number
356  double* &k_pt() {return K_pt;}
357 
358  ///Overload the output function
359  void output(std::ostream &outfile, const unsigned &n_plot)
360  {
361  outfile.precision(16);
362 
363  //Set output Vector
364  Vector<double> s(1);
365 
366  //Tecplot header info
367  outfile << "ZONE I=" << n_plot << std::endl;
368 
369  const unsigned n_node = this->nnode();
370  const unsigned dim = this->dim();
371 
372  Shape psi(n_node);
373  DShape dpsi(n_node,dim);
374 
375  //Loop over plot points
376  for(unsigned l=0;l<n_plot;l++)
377  {
378  s[0] = -1.0 + l*2.0/(n_plot-1);
379 
380  this->dshape_local(s,psi,dpsi);
381  Vector<double> interpolated_tangent(2,0.0);
382  for(unsigned l=0;l<n_node;l++)
383  {
384  const double dpsi_ = dpsi(l,0);
385  for(unsigned i=0;i<2;i++)
386  {
387  interpolated_tangent[i] += this->nodal_position(l,i)*dpsi_;
388  }
389  }
390 
391  //Tangent
392  double t_vel = (this->interpolated_u(s,0)*interpolated_tangent[0] + this->interpolated_u(s,1)*interpolated_tangent[1])/
393  sqrt(interpolated_tangent[0]*interpolated_tangent[0] + interpolated_tangent[1]*interpolated_tangent[1]);
394 
395 
396  //Output the x,y,u,v
397  for(unsigned i=0;i<2;i++) outfile << this->interpolated_x(s,i) << " ";
398  for(unsigned i=0;i<2;i++) outfile << this->interpolated_u(s,i) << " ";
399  //Output a dummy pressure
400  outfile << 0.0 << " ";
401  //Output the concentrations
402  outfile << this->interpolated_C(s) << " ";
403  outfile << interpolated_C_bulk(s) << " ";
404  //Output the interfacial tension
405  outfile << this->sigma(s) << " " << t_vel << std::endl;
406  }
407  outfile << std::endl;
408  }
409 
410 };
411 
412 //==start_of_problem_class===============================================
413 /// Single axisymmetric fluid interface problem including the
414 /// transport of an soluble surfactant.
415 //=======================================================================
416 template<class ELEMENT, class TIMESTEPPER>
417 class InterfaceProblem : public Problem
418 {
419 
420 public:
421 
422  /// Constructor: Pass the number of elements in radial and axial directions
423  /// and the length of the domain in the z direction)
424  InterfaceProblem(const unsigned &n_r, const unsigned &n_z,
425  const double &l_z);
426 
427  /// Destructor (empty)
429 
430  /// Set initial conditions: Set all nodal velocities to zero and
431  /// initialise the previous velocities to correspond to an impulsive
432  /// start
434  {
435  // Determine number of nodes in mesh
436  const unsigned n_node = Bulk_mesh_pt->nnode();
437 
438  // Loop over all nodes in mesh
439  for(unsigned n=0;n<n_node;n++)
440  {
441  // Loop over the three velocity components
442  for(unsigned i=0;i<3;i++)
443  {
444  // Set velocity component i of node n to zero
445  Bulk_mesh_pt->node_pt(n)->set_value(i,0.0);
446  }
447  //Set the bulk concentration to be 1
448  Bulk_mesh_pt->node_pt(n)->set_value(3,1.0);
449  }
450 
451  // Initialise the previous velocity values for timestepping
452  // corresponding to an impulsive start
453  assign_initial_values_impulsive();
454 
455  } // End of set_initial_condition
456 
457 
458  /// The global temporal error norm, based on the movement of the nodes
459  /// in the radial direction only (because that's the only direction
460  /// in which they move!)
462  {
463  //Temp
464  double global_error = 0.0;
465 
466  //Find out how many nodes there are in the problem
467  const unsigned n_node = Bulk_mesh_pt->nnode();
468 
469  //Loop over the nodes and calculate the errors in the positions
470  for(unsigned n=0;n<n_node;n++)
471  {
472  //Set the dimensions to be restricted to the radial direction only
473  const unsigned n_dim = 1;
474  //Set the position error to zero
475  double node_position_error = 0.0;
476  //Loop over the dimensions
477  for(unsigned i=0;i<n_dim;i++)
478  {
479  //Get position error
480  double error =
481  Bulk_mesh_pt->node_pt(n)->position_time_stepper_pt()->
482  temporal_error_in_position(Bulk_mesh_pt->node_pt(n),i);
483 
484  //Add the square of the individual error to the position error
485  node_position_error += error*error;
486  }
487 
488  //Divide the position error by the number of dimensions
489  node_position_error /= n_dim;
490  //Now add to the global error
491  global_error += node_position_error;
492  }
493 
494  //Now the global error must be divided by the number of nodes
495  global_error /= n_node;
496 
497  //Return the square root of the errr
498  return sqrt(global_error);
499  }
500 
501 
502  /// \short Access function for the specific mesh
503  ElasticRectangularQuadMesh<ELEMENT>* Bulk_mesh_pt;
504 
505  /// \short Mesh for the free surface (interface) elements
506  Mesh* Interface_mesh_pt;
507 
508  /// \short Pointer to the constitutive law
509  ConstitutiveLaw* Constitutive_law_pt;
510 
511  /// \short Node for documentatin
512  Node* Document_node_pt;
513 
514  /// Doc the solution
515  void doc_solution(DocInfo &doc_info);
516 
517  /// Do unsteady run up to maximum time t_max with given timestep dt
518  void unsteady_run(const double &t_max, const double &dt);
519 
520  /// Compute the total mass of the insoluble surfactant
522  {
523  //Initialise to zero
524  double mass = 0.0;
525 
526  // Determine number of 1D interface elements in mesh
527  const unsigned n_interface_element = Interface_mesh_pt->nelement();
528 
529  // Loop over the interface elements
530  for(unsigned e=0;e<n_interface_element;e++)
531  {
532  // Upcast from GeneralisedElement to the present element
535  (Interface_mesh_pt->element_pt(e));
536  //Add contribution from each element
537  mass += el_pt->integrate_c();
538  }
539  return mass;
540  } // End of compute_total_mass
541 
542 
543 private:
544 
545  /// Deform the mesh/free surface to a prescribed function
546  void deform_free_surface(const double &epsilon)
547  {
548  //Loop over all nodes in the mesh
549  const unsigned n_node = Bulk_mesh_pt->nnode();
550  for(unsigned n=0;n<n_node;n++)
551  {
552  //Find out the r value
553  double r = Bulk_mesh_pt->node_pt(n)->x(0);
554  // Determine z coordinate of spine
555  double z_value = Bulk_mesh_pt->node_pt(n)->x(1);
556 
557  //Set the thickess of the filme
558  double thickness =
560  + epsilon*cos(Global_Physical_Variables::Alpha*z_value));
561 
562  //Now scale the position accordingly
563  Bulk_mesh_pt->node_pt(n)->x(0) = 1.0 - (1.0-r)*thickness;
564  }
565 
566  //Reset the lagrangian coordinates
567  Bulk_mesh_pt->set_lagrangian_nodal_coordinates();
568  } // End of deform_free_surface
569 
570  /// Trace file
571 ofstream Trace_file;
572 
573 }; // End of problem class
574 
575 
576 
577 //==start_of_constructor==================================================
578 /// Constructor for single fluid interface problem
579 //========================================================================
580 template<class ELEMENT, class TIMESTEPPER>
582 InterfaceProblem(const unsigned &n_r,
583  const unsigned &n_z,
584  const double &l_z)
585 
586 {
587  // Allocate the timestepper (this constructs the time object as well)
588  add_time_stepper_pt(new TIMESTEPPER(true));
589 
590  // Build and assign mesh (the "false" boolean flag tells the mesh
591  // constructor that the domain is not periodic in r)
592  Bulk_mesh_pt = new ElasticRectangularQuadMesh<ELEMENT>(n_r,n_z,1.0,l_z,time_stepper_pt());
593 
594  //Create "surface mesh" that will only contain the interface elements
595  Interface_mesh_pt = new Mesh;
596  {
597  // How many bulk elements are adjacent to boundary b?
598  // Boundary 1 is the outer boundary
599  // The boundaries are labelled
600  // 2
601  // 3 1
602  // 0
603 
604  unsigned n_element = Bulk_mesh_pt->nboundary_element(3);
605 
606  // Loop over the bulk elements adjacent to boundary b?
607  for(unsigned e=0;e<n_element;e++)
608  {
609  // Get pointer to the bulk element that is adjacent to boundary b
610  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
611  Bulk_mesh_pt->boundary_element_pt(3,e));
612 
613  // Find the index of the face of element e along boundary b
614  int face_index = Bulk_mesh_pt->face_index_at_boundary(3,e);
615 
616  // Build the corresponding free surface element
619 
620  //Add the prescribed-flux element to the surface mesh
621  Interface_mesh_pt->add_element_pt(interface_element_pt);
622 
623  } //end of loop over bulk elements adjacent to boundary b
624  }
625 
626 
627  //Use the first node as our document
628  Document_node_pt = Bulk_mesh_pt->node_pt(0);
629 
630  // Add the two sub meshes to the problem
631  add_sub_mesh(Bulk_mesh_pt);
632  add_sub_mesh(Interface_mesh_pt);
633 
634  // Combine all submeshes into a single Mesh
635  build_global_mesh();
636 
637 
638  // --------------------------------------------
639  // Set the boundary conditions for this problem
640  // --------------------------------------------
641 
642  //Pin all azimuthal velocities
643  //and all vertical positions so the
644  //nodes can only move horizontally
645  {
646  unsigned n_node = this->Bulk_mesh_pt->nnode();
647  for(unsigned n=0;n<n_node;++n)
648  {
649  this->Bulk_mesh_pt->node_pt(n)->pin(2);
650  this->Bulk_mesh_pt->node_pt(n)->pin_position(1);
651  }
652  }
653 
654  // All nodes are free by default -- just pin the ones that have
655  // Dirichlet conditions here
656  unsigned ibound = 3;
657  unsigned n_node = Bulk_mesh_pt->nboundary_node(ibound);
658  for(unsigned n=0;n<n_node;n++)
659  {
660  Bulk_mesh_pt->boundary_node_pt(ibound,n)->set_value(4,1.0);
661  }
662 
663  // Determine number of mesh boundaries
664  const unsigned n_boundary = Bulk_mesh_pt->nboundary();
665 
666  // Loop over mesh boundaries
667  for(unsigned b=0;b<n_boundary;b++)
668  {
669  // Determine number of nodes on boundary b
670  const unsigned n_node = Bulk_mesh_pt->nboundary_node(b);
671 
672  // Loop over nodes on boundary b
673  for(unsigned n=0;n<n_node;n++)
674  {
675  // Pin azimuthal velocity on bounds
676  Bulk_mesh_pt->boundary_node_pt(b,n)->pin(2);
677 
678  // Pin velocity on the outer wall
679  if(b==1)
680  {
681  Bulk_mesh_pt->boundary_node_pt(b,n)->pin(0);
682  Bulk_mesh_pt->boundary_node_pt(b,n)->pin(1);
683  Bulk_mesh_pt->boundary_node_pt(b,n)->pin_position(0);
684  }
685  // Pin axial velocity on bottom and top walls (no penetration)
686  if(b==0 || b==2)
687  {
688  Bulk_mesh_pt->boundary_node_pt(b,n)->pin(1);
689  }
690  } // End of loop over nodes on boundary b
691  } // End of loop over mesh boundaries
692 
693  // ----------------------------------------------------------------
694  // Complete the problem setup to make the elements fully functional
695  // ----------------------------------------------------------------
696  Constitutive_law_pt = new GeneralisedHookean(&Global_Physical_Variables::Nu);
697 
698  // Determine number of bulk elements in mesh
699  const unsigned n_bulk = Bulk_mesh_pt->nelement();
700 
701  // Loop over the bulk elements
702  for(unsigned e=0;e<n_bulk;e++)
703  {
704  // Upcast from GeneralisedElement to the present element
705  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
706 
707  // Set the Reynolds number
708  el_pt->re_pt() = &Global_Physical_Variables::Re;
709 
710  // Set the Womersley number
711  el_pt->re_st_pt() = &Global_Physical_Variables::ReSt;
712 
713  // Set the product of the Reynolds number and the inverse of the
714  // Froude number
715  el_pt->re_invfr_pt() = &Global_Physical_Variables::ReInvFr;
716 
717  // Set the direction of gravity
718  el_pt->g_pt() = &Global_Physical_Variables::G;
719 
720  // Set the constitutive law
721  el_pt->constitutive_law_pt() = Constitutive_law_pt;
722 
723  //Set the peclet numbers
725 
727  //Set the diffusion scale
728  el_pt->d_pt() = &Global_Physical_Variables::Diff;
729  } // End of loop over bulk elements
730 
731  // Create a Data object whose single value stores the external pressure
732  Data* external_pressure_data_pt = new Data(1);
733 
734  // Pin and set the external pressure to some arbitrary value
735  double p_ext = Global_Physical_Variables::P_ext;
736 
737  external_pressure_data_pt->pin(0);
738  external_pressure_data_pt->set_value(0,p_ext);
739 
740  // Determine number of 1D interface elements in mesh
741  const unsigned n_interface_element = Interface_mesh_pt->nelement();
742 
743  // Loop over the interface elements
744  for(unsigned e=0;e<n_interface_element;e++)
745  {
746  // Upcast from GeneralisedElement to the present element
749  (Interface_mesh_pt->element_pt(e));
750 
751  // Set the Capillary number
753 
754  // Pass the Data item that contains the single external pressure value
755  el_pt->set_external_pressure_data(external_pressure_data_pt);
756 
757  // Set the surface elasticity number
759 
760  // Set the sorption
762 
763  // Set the K parameter
765 
766  // Set the surface peclect number
768 
769  // Set the surface peclect number multiplied by strouhal number
771 
772  } // End of loop over interface elements
773 
774  // Setup equation numbering scheme
775  cout << "Number of equations: " << assign_eqn_numbers() << std::endl;
776 
777 } // End of constructor
778 
779 
780 
781 //==start_of_doc_solution=================================================
782 /// Doc the solution
783 //========================================================================
784 template<class ELEMENT, class TIMESTEPPER>
786 doc_solution(DocInfo &doc_info)
787 {
788 
789  // Output the time
790  double t= time_pt()->time();
791  cout << "Time is now " << t << std::endl;
792 
793  // Document in trace file
794  Trace_file << time_pt()->time() << " "
795  << Document_node_pt->x(0)
796  << " " << this->compute_total_mass() << std::endl;
797 
798  ofstream some_file;
799  char filename[100];
800 
801  // Set number of plot points (in each coordinate direction)
802  const unsigned npts = 5;
803 
804  // Open solution output file
805 // sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
806 // doc_info.number());
807 // some_file.open(filename);
808 
809  // Output solution to file
810 // Bulk_mesh_pt->output(some_file,npts);
811 // some_file.close();
812  //Put interface in separate file
813  sprintf(filename,"%s/int%i.dat",doc_info.directory().c_str(),
814  doc_info.number());
815  some_file.open(filename);
816  Interface_mesh_pt->output(some_file,npts);
817  some_file.close();
818 
819  // Write file as a tecplot text object...
820 // some_file << "TEXT X=2.5,Y=93.6,F=HELV,HU=POINT,C=BLUE,H=26,T=\"time = "
821  // << time_pt()->time() << "\"";
822  // ...and draw a horizontal line whose length is proportional
823  // to the elapsed time
824  //some_file << "GEOMETRY X=2.5,Y=98,T=LINE,C=BLUE,LT=0.4" << std::endl;
825  //some_file << "1" << std::endl;
826  //some_file << "2" << std::endl;
827  //some_file << " 0 0" << std::endl;
828  //some_file << time_pt()->time()*20.0 << " 0" << std::endl;
829 
830  // Close solution output file
831 // some_file.close();
832 
833  // Output solution to file in paraview format
834  //sprintf(filename,"%s/soln%i.vtu",doc_info.directory().c_str(),
835  // doc_info.number());
836  //some_file.open(filename);
837  //Bulk_mesh_pt->output_paraview(some_file,npts);
838  //some_file.close();
839 
840  // Write pvd information
841  //string file_name="soln"+StringConversion::to_string(doc_info.number())
842  // +".vtu";
843  //ParaviewHelper::write_pvd_information(Global_Physical_Variables::Pvd_file,
844  // file_name,t);
845 
846 } // End of doc_solution
847 
848 
849 
850 //==start_of_unsteady_run=================================================
851 /// Perform run up to specified time t_max with given timestep dt
852 //========================================================================
853 template<class ELEMENT, class TIMESTEPPER>
855 unsteady_run(const double &t_max, const double &dt)
856 {
857 
858  // Set value of epsilon
859  double epsilon = Global_Physical_Variables::Epsilon;
860 
861  // Deform the mesh/free surface
862  deform_free_surface(epsilon);
863 
864  // Initialise DocInfo object
865  DocInfo doc_info;
866 
867  // Set output directory
868  doc_info.set_directory("RESLT");
869 
870  // Initialise counter for solutions
871  doc_info.number()=0;
872 
873  // Open trace file
874  char filename[100];
875  sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
876  Trace_file.open(filename);
877 
878  // Initialise trace file
879  Trace_file << "time" << ", "
880  << "edge spine height" << ", "
881  << "mass " << ", " << std::endl;
882 
883  // Initialise timestep
884  initialise_dt(dt);
885 
886  // Set initial conditions
887  set_initial_condition();
888 
889  // Determine number of timesteps
890  const unsigned n_timestep = unsigned(t_max/dt);
891 
892  // Open pvd file -- a wrapper for all the different
893  // vtu output files plus information about continuous time
894  // to facilitate animations in paraview
895  //sprintf(filename,"%s/soln.pvd",doc_info.directory().c_str());
896  //Global_Physical_Variables::Pvd_file.open(filename);
897  //ParaviewHelper::write_pvd_header(Global_Physical_Variables::Pvd_file);
898 
899  // Doc initial solution
900  doc_solution(doc_info);
901 
902  // Increment counter for solutions
903  doc_info.number()++;
904 
905  //double dt_desired = dt;
906 
907  // Timestepping loop
908  for(unsigned t=1;t<=n_timestep;t++)
909  {
910  // Output current timestep to screen
911  cout << "\nTimestep " << t << " of " << n_timestep << std::endl;
912 
913  // Take one fixed timestep
914  unsteady_newton_solve(dt);
915 
916  //double dt_actual =
917  // adaptive_unsteady_newton_solve(dt_desired,1.0e-6);
918  //dt_desired = dt_actual;
919 
920 
921  // Doc solution
922  doc_solution(doc_info);
923 
924  // Increment counter for solutions
925  doc_info.number()++;
926 
927  //Reset the lagrangian coordinates
928  Bulk_mesh_pt->set_lagrangian_nodal_coordinates();
929 
930  } // End of timestepping loop
931 
932  // write footer and close pvd file
933  //ParaviewHelper::write_pvd_footer(Global_Physical_Variables::Pvd_file);
934  //Global_Physical_Variables::Pvd_file.close();
935 
936 } // End of unsteady_run
937 
938 
939 ///////////////////////////////////////////////////////////////////////
940 ///////////////////////////////////////////////////////////////////////
941 ///////////////////////////////////////////////////////////////////////
942 
943 
944 //==start_of_main=========================================================
945 /// Driver code for single fluid axisymmetric horizontal interface problem
946 //========================================================================
947 int main(int argc, char* argv[])
948 {
949 
950  // Store command line arguments
951  CommandLineArgs::setup(argc,argv);
952 
953  /// Maximum time
954  double t_max = 1000.0;
955 
956  /// Duration of timestep
957  const double dt = 0.1;
958 
959  // If we are doing validation run, use smaller number of timesteps
960  if(CommandLineArgs::Argc>1)
961  {
962  t_max = 0.5;
963  }
964 
965  // Number of elements in radial (r) direction
966  const unsigned n_r = 10;
967 
968  // Number of elements in axial (z) direction
969  const unsigned n_z = 80;
970 
971  // Height of domain
972  const double l_z = MathematicalConstants::Pi/Global_Physical_Variables::Alpha;
973 
974  // Set direction of gravity (vertically downwards)
978 
979  // Set up the spine test problem with AxisymmetricQCrouzeixRaviartElements,
980  // using the BDF<2> timestepper
982  problem(n_r,n_z,l_z);
983 
984  // Run the unsteady simulation
985  problem.unsteady_run(t_max,dt);
986 } // End of main
987 
double K
K parameter that describes solubility number.
double Peclet_S
Surface Peclet number.
ElasticAxisymmetricSolubleSurfactantTransportInterfaceElement(FiniteElement *const &element_pt, const int &face_index)
int main(int argc, char *argv[])
Driver code for single fluid axisymmetric horizontal interface problem.
double Peclet_St_S
Sufrace Peclet number multiplied by Strouhal number.
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.
ofstream Pvd_file
Pvd file – a wrapper for all the different vtu output files plus information about continuous time to...
double dc_bulk_dt_surface(const unsigned &l) const
The time derivative of the bulk surface concentration.
void add_additional_residual_contributions_interface(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag, const Shape &psif, const DShape &dpsifds, const DShape &dpsifdS, const DShape &dpsifdS_div, const Vector< double > &s, const Vector< double > &interpolated_x, const Vector< double > &interpolated_n, const double &W, const double &J)
Helper function to calculate the additional contributions to be added at each integration point...
double ReSt
Womersley = Reynolds times Strouhal.
double *& alpha_pt()
Access function for pointer to adsorption number.
double integrate_c()
Compute the concentration intergated over the surface area.
void output(std::ostream &outfile, const unsigned &n_plot)
Overload the output function.
double *& k_pt()
Access function for pointer to solubility number.
ElasticRectangularQuadMesh< ELEMENT > * Bulk_mesh_pt
Access function for the specific mesh.
double Epsilon
Free surface cosine deformation parameter.
double *& beta_pt()
Access function for pointer to the Elasticity number.
double ReInvFr
Product of Reynolds and Froude number.
Vector< double > G(3)
Direction of gravity.
void unsteady_run(const unsigned &nstep)
Run an unsteady simulation with specified number of steps.
double Beta
Surface Elasticity number (weak case)
void deform_free_surface(const double &epsilon)
Deform the mesh/free surface to a prescribed function.
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.
double interpolated_C_bulk(const Vector< double > &s)
Get the bulk surfactant concentration.
void add_additional_residual_contributions_interface(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag, const Shape &psif, const DShape &dpsifds, const DShape &dpsifdS, const DShape &dpsifdS_div, const Vector< double > &s, const Vector< double > &interpolated_x, const Vector< double > &interpolated_n, const double &W, const double &J)
Overload the Helper function to calculate the residuals and jacobian entries. This particular functio...
ConstitutiveLaw * Constitutive_law_pt
Pointer to the constitutive law.
double Alpha_absorption
Alpha for absorption kinetics.
double Alpha
Wavelength of the domain.
double compute_total_mass()
Compute the total mass of the insoluble surfactant.
double *& ca_pt()
Pointer to the Capillary number.