rayleigh_instability_insoluble_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 an insoluble surfactant. The problem is described in
33 // A 2-D model of Rayleigh instability in capillary tubes --- surfactant
34 // effects by D. Campana, J. Di Paolo & F. A. Saita, Int. J. Multiphase Flow,
35 // vol 30, pp 431--454, (2004).
36 // The initial version of this code was developed by Aman Rajvardhan.
37 
38 // Generic oomph-lib header
39 #include "generic.h"
40 
41 // Axisymmetric Navier-Stokes headers
42 #include "axisym_navier_stokes.h"
43 
44 // Interface headers
45 #include "fluid_interface.h"
46 
47 // The mesh, including horizontal spines
48 #include "meshes/horizontal_single_layer_spine_mesh.h"
49 
50 //Use the oomph and std namespaces
51 using namespace oomph;
52 using namespace std;
53 
54 //==start_of_namespace===================================================
55 /// Namespace for physical parameters
56 /// The parameter values are chosen to be those used in Figures 8, 9
57 /// in Campana et al.
58 //=======================================================================
59 namespace Global_Physical_Variables
60 {
61 
62  //Film thickness parameter
63  double Film_Thickness = 0.2;
64 
65  /// Reynolds number
66  double Re = 40.0;
67 
68  /// Womersley number
69  double ReSt = Re; // (St = 1)
70 
71  /// Product of Reynolds number and inverse of Froude number
72  double ReInvFr = 0.0; // (Fr = 0)
73 
74  /// Capillary number
75  double Ca = pow(Film_Thickness,3.0);
76 
77  /// External pressure
78  double P_ext = 0.0;
79 
80  /// Direction of gravity
81  Vector<double> G(3);
82 
83  /// Wavelength of the domain
84  double Alpha = 1.047;
85 
86  /// Free surface cosine deformation parameter
87  double Epsilon = 1.0e-3;
88 
89  /// \short Surface Elasticity number (weak case)
90  double Beta = 3.6e-3;
91 
92  /// \short Surface Peclet number
93  double Peclet_S = 4032.0;
94 
95  /// \short Sufrace Peclet number multiplied by Strouhal number
96  double Peclet_St_S = 1.0;
97 
98  /// \short Pvd file -- a wrapper for all the different
99  /// vtu output files plus information about continuous time
100  /// to facilitate animations in paraview
101  ofstream Pvd_file;
102 
103 } // End of namespace
104 
105 
106 
107 namespace oomph
108 {
109 
110 //==================================================================
111 ///Spine-based Marangoni surface tension elements that add
112 ///a linear dependence on concentration
113 ///of a surface chemical to the surface tension,
114 ///which decreases with increasing concentration.
115 ///The non-dimensionalisation is the same as Campana et al (2004)
116 ///but we may wish to revisit this.
117 //=================================================================
118 template<class ELEMENT>
121 {
122 private:
123  /// Pointer to an Elasticity number
124  double *Beta_pt;
125 
126  /// Pointer to Surface Peclet number
127  double *Peclet_S_pt;
128 
129  /// Pointer to the surface Peclect Strouhal number
131 
132  /// Index at which the surfactant concentration is stored at the
133  /// nodes
134  unsigned C_index;
135 
136  /// Default value of the physical constants
138 
139 protected:
140 
141  ///Get the surfactant concentration
142  double interpolated_C(const Vector<double> &s)
143  {
144  //Find number of nodes
145  unsigned n_node = this->nnode();
146 
147  //Get the nodal index at which the unknown is stored
148  const unsigned c_index = C_index;
149 
150  //Local shape function
151  Shape psi(n_node);
152 
153  //Find values of shape function
154  this->shape(s,psi);
155 
156  //Initialise value of C
157  double C = 0.0;
158 
159  //Loop over the local nodes and sum
160  for(unsigned l=0;l<n_node;l++)
161  {
162  C += this->nodal_value(l,c_index)*psi(l);
163  }
164 
165  return(C);
166  }
167 
168  /// The time derivative of the surface concentration
169  double dcdt_surface(const unsigned &l) const
170  {
171  // Get the data's timestepper
172  TimeStepper* time_stepper_pt= this->node_pt(l)->time_stepper_pt();
173 
174  //Initialise dudt
175  double dcdt=0.0;
176  //Loop over the timesteps, if there is a non Steady timestepper
177  if (time_stepper_pt->type()!="Steady")
178  {
179  //Find the index at which the variable is stored
180  const unsigned c_index = C_index;
181 
182  // Number of timsteps (past & present)
183  const unsigned n_time = time_stepper_pt->ntstorage();
184 
185  for(unsigned t=0;t<n_time;t++)
186  {
187  dcdt += time_stepper_pt->weight(1,t)*this->nodal_value(t,l,c_index);
188  }
189  }
190  return dcdt;
191  }
192 
193  /// The surface tension function is linear in the
194  /// concentration with constant of proportionality equal
195  /// to the elasticity number.
196  double sigma(const Vector<double> &s)
197  {
198  //Find the number of shape functions
199  const unsigned n_node = this->nnode();
200  //Now get the shape fuctions at the local coordinate
201  Shape psi(n_node);
202  this->shape(s,psi);
203 
204  //Now interpolate the surfactant concentration
205  double C=0.0;
206  for(unsigned l=0;l<n_node;l++)
207  {
208  C += this->nodal_value(l,C_index)*psi(l);
209  }
210 
211  //Get the Elasticity number
212  double Beta = this->beta();
213  //Return the variable surface tension
214  return (1.0 - Beta*(C-1.0));
215  } // End of sigma
216 
217  /// \short Fill in the contribution to the residuals
218  /// Calculate the contribution to the jacobian
219  void fill_in_contribution_to_jacobian(Vector<double> &residuals,DenseMatrix<double> &jacobian)
220  {
221  //Call the generic routine with the flag set to 1
222  this->fill_in_generic_residual_contribution_interface(residuals,jacobian,1);
223  {
224  //Use finite differences to handle concentration variations in the
225  //bulk equations
226  const unsigned n_node = this->nnode();
227  //Find the number of dofs in the element
228  const unsigned n_dof = this->ndof();
229  //Create newres vector
230  Vector<double> newres(n_dof);
231 
232  //Integer storage for local unknown
233  int local_unknown=0;
234 
235  //Use the default finite difference step
236  const double fd_step = this->Default_fd_jacobian_step;
237 
238  //Loop over the nodes again
239  for(unsigned n=0;n<n_node;n++)
240  {
241  //Get the number of values stored at the node
242  unsigned c_index = this->C_index;
243 
244  //Get the local equation number
245  local_unknown = this->nodal_local_eqn(n,c_index);
246  //If it's not pinned
247  if(local_unknown >= 0)
248  {
249  //Store a pointer to the nodal data value
250  double *value_pt = this->node_pt(n)->value_pt(c_index);
251 
252  //Save the old value of the Nodal data
253  double old_var = *value_pt;
254 
255  //Increment the value of the Nodal data
256  *value_pt += fd_step;
257 
258  //Calculate the new residuals
259  this->get_residuals(newres);
260 
261  //Do finite differences
262  for(unsigned m=0;m<n_dof;m++)
263  {
264  double sum = (newres[m] - residuals[m])/fd_step;
265  //Stick the entry into the Jacobian matrix
266  jacobian(m,local_unknown) = sum;
267  }
268 
269  //Reset the Nodal data
270  *value_pt = old_var;
271  }
272  }
273  }
274 
275  //Call the generic routine to handle the spine variables
276  SpineElement<FaceGeometry<ELEMENT> >::fill_in_jacobian_from_geometric_data(jacobian);
277  }
278 
279 
280  /// \short Overload the Helper function to calculate the residuals and
281  /// jacobian entries. This particular function ensures that the
282  /// additional entries are calculated inside the integration loop
283  void add_additional_residual_contributions_interface(Vector<double> &residuals, DenseMatrix<double> &jacobian,
284  const unsigned &flag,const Shape &psif, const DShape &dpsifds,
285  const DShape &dpsifdS, const DShape &dpsifdS_div,
286  const Vector<double> &s,
287  const Vector<double> &interpolated_x, const Vector<double> &interpolated_n,
288  const double &W,const double &J)
289  {
290  //Flag to control whether the Campana formulation (false)
291  // or our own (true) is used
292  bool Integrated_curvature = true;
293 
294  //Find out how many nodes there are
295  unsigned n_node = this->nnode();
296 
297  //Storage for the local equation numbers and unknowns
298  int local_eqn = 0, local_unknown = 0;
299 
300  //Surface advection-diffusion equation
301 
302  //Find the index at which the concentration is stored
303  unsigned c_index = this->C_index;
304  Vector<unsigned> u_index = this->U_index_interface;
305 
306  //Read out the surface peclect number
307  const double Pe_s = this->peclet_s();
308  //const double PeSt_s = this->peclet_strouhal_s();
309 
310  //Read out the radial position
311  const double r = interpolated_x[0];
312 
313  //Rescale the jacobian to be the "raw" version
314  const double J_raw = J/r;
315 
316  //Now calculate the concentration at this point
317  //Assuming the same shape functions are used (which they are)
318  double interpolated_C = 0.0;
319  double interpolated_dCds = 0.0;
320  double dCdt = 0.0;
321  //The tangent vector
322  const unsigned ndim = this->node_pt(0)->ndim();
323  Vector<double> interpolated_tangent(ndim,0.0);
324  Vector<double> interpolated_u(ndim,0.0);
325  Vector<double> mesh_velocity(ndim,0.0);
326  Vector<double> interpolated_duds(ndim,0.0);
327  if(ndim+1 != u_index.size())
328  {
329  throw OomphLibError("Dimension Incompatibility",
330  OOMPH_CURRENT_FUNCTION,
331  OOMPH_EXCEPTION_LOCATION);
332  }
333 
334  for(unsigned l=0;l<n_node;l++)
335  {
336  const double psi = psif(l);
337  const double dpsi = dpsifds(l,0);
338  interpolated_C += this->nodal_value(l,c_index)*psi;
339  interpolated_dCds += this->nodal_value(l,c_index)*dpsi;
340  dCdt += dcdt_surface(l)*psi;
341  for(unsigned i=0;i<ndim;i++)
342  {
343  interpolated_tangent[i] += this->nodal_position(l,i)*dpsi;
344  interpolated_u[i] += this->nodal_value(l,u_index[i])*psi;
345  interpolated_duds[i] += this->nodal_value(l,u_index[i])*dpsi;
346  }
347  //Mesh Velocity
348  for(unsigned j=0;j<ndim;j++)
349  {
350  mesh_velocity[j] += this->dnodal_position_dt(l,j)*psi;
351  }
352  }
353 
354 
355  double u_tangent = 0.0, t_length = 0.0;
356  for(unsigned i=0;i<ndim;i++)
357  {
358  u_tangent += interpolated_u[i]*interpolated_tangent[i];
359  t_length += interpolated_tangent[i]*interpolated_tangent[i];
360  }
361 
362  //Work out the second derivative of position
363  Vector<double> d2xds(2,0.0);
364  for(unsigned i=0;i<2;i++)
365  {
366  d2xds[i] = this->nodal_position(0,i) +
367  this->nodal_position(2,i) - 2.0*this->nodal_position(1,i);
368  }
369 
370  //Let's do the first component of the curvature
371  double k1 = 0.0;
372  for(unsigned i=0;i<2;i++)
373  {
374  k1 += (d2xds[i]/(J_raw*J_raw) - interpolated_tangent[i]*(
375  interpolated_tangent[0]*d2xds[0]
376  + interpolated_tangent[1]*d2xds[1])/(J_raw*J_raw*J_raw*J_raw))*interpolated_n[i];
377  }
378 
379  //Second component of the curvature
380  double k2 = - (interpolated_n[0] / r);
381 
382  //Now we add the flux term to the appropriate residuals
383  for(unsigned l=0;l<n_node;l++)
384  {
385  //Read out the apprporiate local equation
386  local_eqn = this->nodal_local_eqn(l,c_index);
387 
388  //If not a boundary condition
389  if(local_eqn >= 0)
390  {
391  //Time derivative
392  residuals[local_eqn] += dCdt*psif(l)*r*W*J_raw;
393 
394  double tmp = 0.0;
395  //Compute our version in which second derivatives are not required
396  if(Integrated_curvature)
397  {
398  for(unsigned i=0;i<2;i++)
399  {
400  tmp +=
401  (interpolated_u[i] - mesh_velocity[i])*interpolated_tangent[i];
402  }
403  //First Advection term
404  residuals[local_eqn] += tmp*interpolated_dCds*psif(l)*r*W/J_raw;
405  //Additional term from axisymmetric formulation
406  residuals[local_eqn] += interpolated_C*interpolated_u[0]*psif(l)*W*J_raw;
407  //Second Advection term
408  residuals[local_eqn] += interpolated_C*
409  (interpolated_duds[0]*interpolated_tangent[0]
410  + interpolated_duds[1]*interpolated_tangent[1])*r*W*psif(l)/J_raw;
411  }
412  //This is the Campana formulation
413  else
414  {
415  //ALE term
416  for(unsigned i=0;i<2;i++)
417  {
418  tmp += -mesh_velocity[i]*interpolated_tangent[i];
419  }
420  residuals[local_eqn] += tmp*interpolated_dCds*psif(l)*r*W/J_raw;
421  // Curvature (normal velocity) term
422  residuals[local_eqn] -=
423  interpolated_C*(k1 + k2)*psif(l)*r*W*J_raw*(
424  interpolated_u[0]*interpolated_n[0]
425  + interpolated_u[1]*interpolated_n[1]);
426  //Integrated by parts tangential advection term
427  residuals[local_eqn] -= interpolated_C*(
428  interpolated_tangent[0]*interpolated_u[0] +
429  interpolated_tangent[1]*interpolated_u[1])*dpsifds(l,0)*r*W/J_raw;
430  }
431 
432  //Diffusion term
433  residuals[local_eqn] += (1.0/Pe_s)*interpolated_dCds*dpsifds(l,0)*r*W/J_raw;
434 
435  //We also need to worry about the jacobian terms
436  if(flag)
437  {
438  //Loop over the nodes again
439  for(unsigned l2=0;l2<n_node;l2++)
440  {
441  //Get the time stepper
442  TimeStepper* time_stepper_pt=this->node_pt(l2)->time_stepper_pt();
443 
444  //Get the unknown c_index
445  local_unknown =this->nodal_local_eqn(l2,c_index);
446 
447  if(local_unknown >=0)
448  {
449  jacobian(local_eqn,local_unknown) +=
450  time_stepper_pt->weight(1,0)*psif(l2)*psif(l)*r*W*J_raw;
451 
452  if(Integrated_curvature)
453  {
454  jacobian(local_eqn,local_unknown) += ((interpolated_u[0] - mesh_velocity[0])*interpolated_tangent[0]
455  + (interpolated_u[1] - mesh_velocity[1])*interpolated_tangent[1])*dpsifds(l2,0)
456  *psif(l)*r*W/J_raw;
457 
458 
459  jacobian(local_eqn,local_unknown) += psif(l2)*interpolated_u[0]*psif(l)*W*J_raw;
460 
461  jacobian(local_eqn,local_unknown) += psif(l2)*(interpolated_tangent[0]*interpolated_duds[0]
462  + interpolated_tangent[1]*interpolated_duds[1])*psif(l)*r*W/J_raw;
463  }
464  else
465  {
466  jacobian(local_eqn,local_unknown) -=
467  (mesh_velocity[0]*interpolated_tangent[0] + mesh_velocity[1]*interpolated_tangent[1])*dpsifds(l2,0)*psif(l)*r*W/J_raw;
468 
469  jacobian(local_eqn,local_unknown) -= psif(l2)*(k1 + k2)*psif(l)*r*W*J_raw*(interpolated_u[0]*interpolated_n[0]
470  + interpolated_u[1]*interpolated_n[1]);
471 
472  jacobian(local_eqn,local_unknown) -= psif(l2)*(interpolated_tangent[0]*interpolated_u[0] + interpolated_tangent[1]*
473  interpolated_u[1])*dpsifds(l,0)*r*W/J_raw;
474  }
475 
476  jacobian(local_eqn,local_unknown) += (1.0/Pe_s)*dpsifds(l2,0)*dpsifds(l,0)*r*W/J_raw;
477 
478  }
479 
480 
481  //Loop over the velocity components
482  for(unsigned i2=0;i2<ndim;i2++)
483  {
484 
485  //Get the unknown
486  local_unknown = this->nodal_local_eqn(l2,u_index[i2]);
487 
488 
489  //If not a boundary condition
490  if(local_unknown >= 0)
491  {
492 
493  // jacobian(local_eqn,local_unknown) += time_stepper_pt->weight(1,0)*psif(l2)*PeSt_s*psif(l)*r*W*J_raw;
494  if(Integrated_curvature)
495  {
496  jacobian(local_eqn,local_unknown) +=
497  interpolated_dCds*psif(l2)*interpolated_tangent[i2]*psif(l)*r*W/J_raw;
498 
499  if(i2==0)
500  {
501  jacobian(local_eqn,local_unknown) += interpolated_C*psif(l2)*W*J_raw;
502  }
503  }
504  else
505  {
506  jacobian(local_eqn,local_unknown) -= interpolated_C*(k1 + k2)*psif(l)*r*W*J_raw*psif(l2)*interpolated_n[i2];
507 
508  jacobian(local_eqn,local_unknown) -= interpolated_C*interpolated_tangent[i2]*psif(l2)*dpsifds(l,0)*r*W/J_raw;
509  }
510 
511  jacobian(local_eqn,local_unknown) +=
512  interpolated_C*interpolated_tangent[i2]*dpsifds(l2,0)*psif(l)*r*W/J_raw;
513  }
514  }
515  }
516  }
517  }
518  } //End of loop over the nodes
519 
520  }
521 
522 
523  /// Add the element's contribution to its residuals vector,
524  /// jacobian matrix and mass matrix
525  void fill_in_contribution_to_jacobian_and_mass_matrix( Vector<double> &residuals, DenseMatrix<double> &jacobian,
526  DenseMatrix<double> &mass_matrix)
527  {
528  //Add the contribution to the jacobian
529  this->fill_in_contribution_to_jacobian(residuals,jacobian);
530  //No mass matrix terms, but should probably do kinematic bit here
531  }
532 
533 public:
534  ///Constructor that passes the bulk element and face index down
535  ///to the underlying
536  SpineAxisymmetricMarangoniSurfactantFluidInterfaceElement(FiniteElement* const &element_pt, const int &face_index) : SpineAxisymmetricFluidInterfaceElement<ELEMENT>
537  (element_pt,face_index)
538  {
539  //Initialise the values
540  Beta_pt = &Default_Physical_Constant_Value;
541  Peclet_S_pt = &Default_Physical_Constant_Value;
542  Peclet_Strouhal_S_pt = &Default_Physical_Constant_Value;
543 
544  //Add the additional surfactant terms to these surface elements
545 
546  //Read out the number of nodes on the face
547  //For some reason I need to specify the this pointer here(!)
548  unsigned n_node_face = this->nnode();
549  //Set the additional data values in the face
550  //There is one additional values at each node --- the lagrange multiplier
551  Vector<unsigned> additional_data_values(n_node_face);
552  for(unsigned i=0;i<n_node_face;i++) additional_data_values[i] = 1;
553  //Resize the data arrays accordingly
554  this->resize_nodes(additional_data_values);
555 
556  //The C_index is the new final value
557  //Minor HACK HERE
558  C_index = this->node_pt(0)->nvalue()-1;
559 }
560 
561  ///Return the Elasticity number
562  double beta() {return *Beta_pt;}
563 
564  ///Return the surface peclect number
565  double peclet_s() {return *Peclet_S_pt;}
566 
567  ///Return the surface peclect strouhal number
568  double peclet_strouhal_s() {return *Peclet_Strouhal_S_pt;}
569 
570  ///Access function for pointer to the Elasticity number
571  double* &beta_pt() {return Beta_pt;}
572 
573  ///Access function for pointer to the surface Peclet number
574  double* &peclet_s_pt() {return Peclet_S_pt;}
575 
576  ///Access function for pointer to the surface Peclet x Strouhal number
577  double* &peclet_strouhal_s_pt() {return Peclet_Strouhal_S_pt;}
578 
579  void output(std::ostream &outfile, const unsigned &n_plot)
580  {
581  outfile.precision(16);
582 
583  //Set output Vector
584  Vector<double> s(1);
585 
586  //Tecplot header info
587  outfile << "ZONE I=" << n_plot << std::endl;
588 
589  const unsigned n_node = this->nnode();
590  const unsigned dim = this->dim();
591 
592  Shape psi(n_node);
593  DShape dpsi(n_node,dim);
594 
595  //Loop over plot points
596  for(unsigned l=0;l<n_plot;l++)
597  {
598  s[0] = -1.0 + l*2.0/(n_plot-1);
599 
600  this->dshape_local(s,psi,dpsi);
601  Vector<double> interpolated_tangent(2,0.0);
602  for(unsigned l=0;l<n_node;l++)
603  {
604  const double dpsi_ = dpsi(l,0);
605  for(unsigned i=0;i<2;i++)
606  {
607  interpolated_tangent[i] += this->nodal_position(l,i)*dpsi_;
608  }
609  }
610 
611  //Tangent
612  double t_vel = (this->interpolated_u(s,0)*interpolated_tangent[0] + this->interpolated_u(s,1)*interpolated_tangent[1])/
613  sqrt(interpolated_tangent[0]*interpolated_tangent[0] + interpolated_tangent[1]*interpolated_tangent[1]);
614 
615 
616  //Output the x,y,u,v
617  for(unsigned i=0;i<2;i++) outfile << this->interpolated_x(s,i) << " ";
618  for(unsigned i=0;i<2;i++) outfile << this->interpolated_u(s,i) << " ";
619  //Output a dummy pressure
620  outfile << 0.0 << " ";
621  //Output the concentration
622  outfile << interpolated_C(s) << " ";
623  //Output the interfacial tension
624  outfile << sigma(s) << " " << t_vel << std::endl;
625  }
626  outfile << std::endl;
627  }
628 
629 
630  /// \short Compute the concentration intergated over the area
631  double integrate_c() const
632  {
633  //Find out how many nodes there are
634  unsigned n_node = this->nnode();
635 
636  //Set up memeory for the shape functions
637  Shape psif(n_node);
638  DShape dpsifds(n_node,1);
639 
640  //Set the value of n_intpt
641  unsigned n_intpt = this->integral_pt()->nweight();
642 
643  //Storage for the local coordinate
644  Vector<double> s(1);
645 
646  //Storage for the total mass
647  double mass = 0.0;
648 
649  //Loop over the integration points
650  for(unsigned ipt=0;ipt<n_intpt;ipt++)
651  {
652  //Get the local coordinate at the integration point
653  s[0] = this->integral_pt()->knot(ipt,0);
654 
655  //Get the integral weight
656  double W = this->integral_pt()->weight(ipt);
657 
658  //Call the derivatives of the shape function
659  this->dshape_local_at_knot(ipt,psif,dpsifds);
660 
661  //Define and zero the tangent Vectors and local velocities
662  Vector<double> interpolated_x(2,0.0);
663  Vector<double> interpolated_t(2,0.0);
664  double interpolated_c = 0.0;
665 
666 
667  //Loop over the shape functions to compute concentration and tangent
668  for(unsigned l=0;l<n_node;l++)
669  {
670  interpolated_c += this->nodal_value(l,C_index)*psif(l);
671  //Calculate the tangent vector
672  for(unsigned i=0;i<2;i++)
673  {
674  interpolated_x[i] += this->nodal_position(l,i)*psif(l);
675  interpolated_t[i] += this->nodal_position(l,i)*dpsifds(l,0);
676  }
677  }
678 
679  //The first positional coordinate is the radial coordinate
680  double r = interpolated_x[0];
681 
682  //Calculate the length of the tangent Vector
683  double tlength = interpolated_t[0]*interpolated_t[0] +
684  interpolated_t[1]*interpolated_t[1];
685 
686  //Set the Jacobian of the line element
687  double J = sqrt(tlength);
688 
689  mass += interpolated_c*r*W*J;
690  }
691  return mass;
692  }
693 
694 };
695 
696 
697 //Define the default physical value to be one
698 template<class ELEMENT>
700 
701 }
702 
703 
704 namespace oomph
705 {
706 
707 //======================================================================
708 /// Inherit from the standard Horizontal single-layer SpineMesh
709 /// and modify the spine_node_update() function so that it is appropriate
710 /// for an annular film, rather than a fluid fibre.
711 //======================================================================
712 template <class ELEMENT>
714  public HorizontalSingleLayerSpineMesh<ELEMENT>
715 {
716 
717 public:
718 
719  /// \short Constructor: Pass number of elements in x-direction, number of
720  /// elements in y-direction, radial extent, axial length , and pointer
721  /// to timestepper (defaults to Steady timestepper)
722  MyHorizontalSingleLayerSpineMesh(const unsigned &nx,
723  const unsigned &ny,
724  const double &lx,
725  const double &ly,
726  TimeStepper* time_stepper_pt=
727  &Mesh::Default_TimeStepper) :
728  HorizontalSingleLayerSpineMesh<ELEMENT>(nx,ny,lx,ly,time_stepper_pt) {}
729 
730  /// \short Node update function assumed spines rooted at the wall
731  /// fixed to be at r=1 and directed inwards to r=0.
732  virtual void spine_node_update(SpineNode* spine_node_pt)
733  {
734  //Get fraction along the spine
735  double W = spine_node_pt->fraction();
736 
737  //Get spine height
738  double H = spine_node_pt->h();
739 
740  //Set the value of y
741  spine_node_pt->x(0) = 1.0-(1.0-W)*H;
742  }
743 
744 };
745 
746 }
747 
748 
749 //==start_of_problem_class===============================================
750 /// Single axisymmetric fluid interface problem including the
751 /// transport of an insoluble surfactant.
752 //=======================================================================
753 template<class ELEMENT, class TIMESTEPPER>
754 class InterfaceProblem : public Problem
755 {
756 
757 public:
758 
759  /// Constructor: Pass the number of elements in radial and axial directions
760  /// and the length of the domain in the z direction)
761  InterfaceProblem(const unsigned &n_r, const unsigned &n_z,
762  const double &l_z);
763 
764  /// Destructor (empty)
766 
767  /// Spine heights/lengths are unknowns in the problem so their values get
768  /// corrected during each Newton step. However, changing their value does
769  /// not automatically change the nodal positions, so we need to update all
770  /// of them here.
772  {
773  Bulk_mesh_pt->node_update();
774  }
775 
776  /// Set initial conditions: Set all nodal velocities to zero and
777  /// initialise the previous velocities to correspond to an impulsive
778  /// start
780  {
781  // Determine number of nodes in mesh
782  const unsigned n_node = Bulk_mesh_pt->nnode();
783 
784  // Loop over all nodes in mesh
785  for(unsigned n=0;n<n_node;n++)
786  {
787  // Loop over the three velocity components
788  for(unsigned i=0;i<3;i++)
789  {
790  // Set velocity component i of node n to zero
791  Bulk_mesh_pt->node_pt(n)->set_value(i,0.0);
792  }
793  }
794 
795  // Initialise the previous velocity values for timestepping
796  // corresponding to an impulsive start
797  assign_initial_values_impulsive();
798 
799  } // End of set_initial_condition
800 
801 
802  /// The global temporal error norm, based on the movement of the nodes
803  /// in the radial direction only (because that's the only direction
804  /// in which they move!)
806  {
807  //Temp
808  double global_error = 0.0;
809 
810  //Find out how many nodes there are in the problem
811  const unsigned n_node = Bulk_mesh_pt->nnode();
812 
813  //Loop over the nodes and calculate the errors in the positions
814  for(unsigned n=0;n<n_node;n++)
815  {
816  //Set the dimensions to be restricted to the radial direction only
817  const unsigned n_dim = 1;
818  //Set the position error to zero
819  double node_position_error = 0.0;
820  //Loop over the dimensions
821  for(unsigned i=0;i<n_dim;i++)
822  {
823  //Get position error
824  double error =
825  Bulk_mesh_pt->node_pt(n)->position_time_stepper_pt()->
826  temporal_error_in_position(Bulk_mesh_pt->node_pt(n),i);
827 
828  //Add the square of the individual error to the position error
829  node_position_error += error*error;
830  }
831 
832  //Divide the position error by the number of dimensions
833  node_position_error /= n_dim;
834  //Now add to the global error
835  global_error += node_position_error;
836  }
837 
838  //Now the global error must be divided by the number of nodes
839  global_error /= n_node;
840 
841  //Return the square root of the errr
842  return sqrt(global_error);
843  }
844 
845 
846  /// \short Access function for the specific mesh
848 
849  /// \short Mesh for the free surface (interface) elements
851 
852  /// Doc the solution
853  void doc_solution(DocInfo &doc_info);
854 
855  /// Do unsteady run up to maximum time t_max with given timestep dt
856  void unsteady_run(const double &t_max, const double &dt);
857 
858  /// Compute the total mass of the insoluble surfactant
860  {
861  //Initialise to zero
862  double mass = 0.0;
863 
864  // Determine number of 1D interface elements in mesh
865  const unsigned n_interface_element = Interface_mesh_pt->nelement();
866 
867  // Loop over the interface elements
868  for(unsigned e=0;e<n_interface_element;e++)
869  {
870  // Upcast from GeneralisedElement to the present element
873  (Interface_mesh_pt->element_pt(e));
874  //Add contribution from each element
875  mass += el_pt->integrate_c();
876  }
877  return mass;
878  } // End of compute_total_mass
879 
880 
881 private:
882 
883  /// Deform the mesh/free surface to a prescribed function
884  void deform_free_surface(const double &epsilon)
885  {
886  // Determine number of spines in mesh
887  const unsigned n_spine = Bulk_mesh_pt->nspine();
888 
889  // Loop over spines in mesh
890  for(unsigned i=0;i<n_spine;i++)
891  {
892 
893  // Determine z coordinate of spine
894  double z_value = Bulk_mesh_pt->boundary_node_pt(1,i)->x(1);
895 
896  // Set spine height
897  Bulk_mesh_pt->spine_pt(i)->height() =
899  + epsilon*cos(Global_Physical_Variables::Alpha*z_value));
900 
901  } // End of loop over spines
902 
903  // Update nodes in bulk mesh
904  Bulk_mesh_pt->node_update();
905 
906  } // End of deform_free_surface
907 
908  /// Trace file
909  ofstream Trace_file;
910 
911 }; // End of problem class
912 
913 
914 
915 //==start_of_constructor==================================================
916 /// Constructor for single fluid interface problem
917 //========================================================================
918 template<class ELEMENT, class TIMESTEPPER>
920 InterfaceProblem(const unsigned &n_r,
921  const unsigned &n_z,
922  const double &l_z)
923 
924 {
925  // Allocate the timestepper (this constructs the time object as well)
926  add_time_stepper_pt(new TIMESTEPPER(true));
927 
928  // Build and assign mesh (the "false" boolean flag tells the mesh
929  // constructor that the domain is not periodic in r)
930  Bulk_mesh_pt = new MyHorizontalSingleLayerSpineMesh<ELEMENT>(n_r,n_z,1.0,l_z,time_stepper_pt());
931 
932  //Create "surface mesh" that will only contain the interface elements
933  Interface_mesh_pt = new Mesh;
934  {
935  // How many bulk elements are adjacent to boundary b?
936  // Boundary 1 is the outer boundary
937  // The boundaries are labelled
938  // 2
939  // 3 1
940  // 0
941 
942  unsigned n_element = Bulk_mesh_pt->nboundary_element(3);
943 
944  // Loop over the bulk elements adjacent to boundary b?
945  for(unsigned e=0;e<n_element;e++)
946  {
947  // Get pointer to the bulk element that is adjacent to boundary b
948  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
949  Bulk_mesh_pt->boundary_element_pt(3,e));
950 
951  // Find the index of the face of element e along boundary b
952  int face_index = Bulk_mesh_pt->face_index_at_boundary(3,e);
953 
954  // Build the corresponding free surface element
957 
958  //Add the prescribed-flux element to the surface mesh
959  Interface_mesh_pt->add_element_pt(interface_element_pt);
960 
961  } //end of loop over bulk elements adjacent to boundary b
962  }
963 
964 
965  // Add the two sub meshes to the problem
966  add_sub_mesh(Bulk_mesh_pt);
967  add_sub_mesh(Interface_mesh_pt);
968 
969  // Combine all submeshes into a single Mesh
970  build_global_mesh();
971 
972 
973  // --------------------------------------------
974  // Set the boundary conditions for this problem
975  // --------------------------------------------
976 
977  //Pin all azimuthal velocities
978  {
979  unsigned n_node = this->Bulk_mesh_pt->nnode();
980  for(unsigned n=0;n<n_node;++n)
981  {
982  this->Bulk_mesh_pt->node_pt(n)->pin(2);
983  }
984  }
985 
986  // All nodes are free by default -- just pin the ones that have
987  // Dirichlet conditions here
988  unsigned ibound = 3;
989  unsigned n_node = Bulk_mesh_pt->nboundary_node(ibound);
990  for(unsigned n=0;n<n_node;n++)
991  {
992  Bulk_mesh_pt->boundary_node_pt(ibound,n)->set_value(3,1.0);
993  }
994 
995  // Determine number of mesh boundaries
996  const unsigned n_boundary = Bulk_mesh_pt->nboundary();
997 
998  // Loop over mesh boundaries
999  for(unsigned b=0;b<n_boundary;b++)
1000  {
1001  // Determine number of nodes on boundary b
1002  const unsigned n_node = Bulk_mesh_pt->nboundary_node(b);
1003 
1004  // Loop over nodes on boundary b
1005  for(unsigned n=0;n<n_node;n++)
1006  {
1007  // Pin azimuthal velocity on bounds
1008  Bulk_mesh_pt->boundary_node_pt(b,n)->pin(2);
1009 
1010  // Pin velocity on the outer wall
1011  if(b==1)
1012  {
1013  Bulk_mesh_pt->boundary_node_pt(b,n)->pin(0);
1014  Bulk_mesh_pt->boundary_node_pt(b,n)->pin(1);
1015  }
1016  // Pin axial velocity on bottom and top walls (no penetration)
1017  if(b==0 || b==2)
1018  {
1019  Bulk_mesh_pt->boundary_node_pt(b,n)->pin(1);
1020  }
1021  } // End of loop over nodes on boundary b
1022  } // End of loop over mesh boundaries
1023 
1024  // ----------------------------------------------------------------
1025  // Complete the problem setup to make the elements fully functional
1026  // ----------------------------------------------------------------
1027 
1028  // Determine number of bulk elements in mesh
1029  const unsigned n_bulk = Bulk_mesh_pt->nelement();
1030 
1031  // Loop over the bulk elements
1032  for(unsigned e=0;e<n_bulk;e++)
1033  {
1034  // Upcast from GeneralisedElement to the present element
1035  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));
1036 
1037  // Set the Reynolds number
1038  el_pt->re_pt() = &Global_Physical_Variables::Re;
1039 
1040  // Set the Womersley number
1041  el_pt->re_st_pt() = &Global_Physical_Variables::ReSt;
1042 
1043  // Set the product of the Reynolds number and the inverse of the
1044  // Froude number
1045  el_pt->re_invfr_pt() = &Global_Physical_Variables::ReInvFr;
1046 
1047  // Set the direction of gravity
1048  el_pt->g_pt() = &Global_Physical_Variables::G;
1049 
1050  } // End of loop over bulk elements
1051 
1052  // Create a Data object whose single value stores the external pressure
1053  Data* external_pressure_data_pt = new Data(1);
1054 
1055  // Pin and set the external pressure to some arbitrary value
1056  double p_ext = Global_Physical_Variables::P_ext;
1057 
1058  external_pressure_data_pt->pin(0);
1059  external_pressure_data_pt->set_value(0,p_ext);
1060 
1061  // Determine number of 1D interface elements in mesh
1062  const unsigned n_interface_element = Interface_mesh_pt->nelement();
1063 
1064  // Loop over the interface elements
1065  for(unsigned e=0;e<n_interface_element;e++)
1066  {
1067  // Upcast from GeneralisedElement to the present element
1070  (Interface_mesh_pt->element_pt(e));
1071 
1072  // Set the Capillary number
1074 
1075  // Pass the Data item that contains the single external pressure value
1076  el_pt->set_external_pressure_data(external_pressure_data_pt);
1077 
1078  // Set the surface elasticity number
1080 
1081  // Set the surface peclect number
1083 
1084  // Set the surface peclect number multiplied by strouhal number
1086 
1087  } // End of loop over interface elements
1088 
1089  // Setup equation numbering scheme
1090  cout << "Number of equations: " << assign_eqn_numbers() << std::endl;
1091 
1092 } // End of constructor
1093 
1094 
1095 
1096 //==start_of_doc_solution=================================================
1097 /// Doc the solution
1098 //========================================================================
1099 template<class ELEMENT, class TIMESTEPPER>
1101 doc_solution(DocInfo &doc_info)
1102 {
1103 
1104  // Output the time
1105  double t= time_pt()->time();
1106  cout << "Time is now " << t << std::endl;
1107 
1108  // Document in trace file
1109  Trace_file << time_pt()->time() << " "
1110  << Bulk_mesh_pt->spine_pt(0)->height()
1111  << " " << this->compute_total_mass() << std::endl;
1112 
1113  ofstream some_file;
1114  char filename[100];
1115 
1116  // Set number of plot points (in each coordinate direction)
1117  const unsigned npts = 5;
1118 
1119  // Open solution output file
1120  //sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
1121  // doc_info.number());
1122  //some_file.open(filename);
1123 
1124  // Output solution to file
1125  //Bulk_mesh_pt->output(some_file,npts);
1126  //some_file.close();
1127  //Put interface in separate file
1128  sprintf(filename,"%s/int%i.dat",doc_info.directory().c_str(),
1129  doc_info.number());
1130  some_file.open(filename);
1131  Interface_mesh_pt->output(some_file,npts);
1132  some_file.close();
1133 
1134  // Write file as a tecplot text object...
1135 // some_file << "TEXT X=2.5,Y=93.6,F=HELV,HU=POINT,C=BLUE,H=26,T=\"time = "
1136  // << time_pt()->time() << "\"";
1137  // ...and draw a horizontal line whose length is proportional
1138  // to the elapsed time
1139  //some_file << "GEOMETRY X=2.5,Y=98,T=LINE,C=BLUE,LT=0.4" << std::endl;
1140  //some_file << "1" << std::endl;
1141  //some_file << "2" << std::endl;
1142  //some_file << " 0 0" << std::endl;
1143  //some_file << time_pt()->time()*20.0 << " 0" << std::endl;
1144 
1145  // Close solution output file
1146 // some_file.close();
1147 
1148  // Output solution to file in paraview format
1149  //sprintf(filename,"%s/soln%i.vtu",doc_info.directory().c_str(),
1150  // doc_info.number());
1151  //some_file.open(filename);
1152  //Bulk_mesh_pt->output_paraview(some_file,npts);
1153  //some_file.close();
1154 
1155  // Write pvd information
1156  string file_name="soln"+StringConversion::to_string(doc_info.number())
1157  +".vtu";
1158  ParaviewHelper::write_pvd_information(Global_Physical_Variables::Pvd_file,
1159  file_name,t);
1160 
1161 } // End of doc_solution
1162 
1163 
1164 
1165 //==start_of_unsteady_run=================================================
1166 /// Perform run up to specified time t_max with given timestep dt
1167 //========================================================================
1168 template<class ELEMENT, class TIMESTEPPER>
1170 unsteady_run(const double &t_max, const double &dt)
1171 {
1172 
1173  // Set value of epsilon
1174  double epsilon = Global_Physical_Variables::Epsilon;
1175 
1176  // Deform the mesh/free surface
1177  deform_free_surface(epsilon);
1178 
1179  // Initialise DocInfo object
1180  DocInfo doc_info;
1181 
1182  // Set output directory
1183  doc_info.set_directory("RESLT");
1184 
1185  // Initialise counter for solutions
1186  doc_info.number()=0;
1187 
1188  // Open trace file
1189  char filename[100];
1190  sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
1191  Trace_file.open(filename);
1192 
1193  // Initialise trace file
1194  Trace_file << "time" << ", "
1195  << "edge spine height" << ", "
1196  << "contact angle left" << ", "
1197  << "contact angle right" << ", " << std::endl;
1198 
1199  // Initialise timestep
1200  initialise_dt(dt);
1201 
1202  // Set initial conditions
1203  set_initial_condition();
1204 
1205  // Determine number of timesteps
1206  const unsigned n_timestep = unsigned(t_max/dt);
1207 
1208  // Open pvd file -- a wrapper for all the different
1209  // vtu output files plus information about continuous time
1210  // to facilitate animations in paraview
1211  sprintf(filename,"%s/soln.pvd",doc_info.directory().c_str());
1212  Global_Physical_Variables::Pvd_file.open(filename);
1213  ParaviewHelper::write_pvd_header(Global_Physical_Variables::Pvd_file);
1214 
1215  // Doc initial solution
1216  doc_solution(doc_info);
1217 
1218  // Increment counter for solutions
1219  doc_info.number()++;
1220 
1221  //double dt_desired = dt;
1222 
1223  // Timestepping loop
1224  for(unsigned t=1;t<=n_timestep;t++)
1225  {
1226  // Output current timestep to screen
1227  cout << "\nTimestep " << t << " of " << n_timestep << std::endl;
1228 
1229  // Take one fixed timestep
1230  unsteady_newton_solve(dt);
1231 
1232  //double dt_actual =
1233  // adaptive_unsteady_newton_solve(dt_desired,1.0e-6);
1234  //dt_desired = dt_actual;
1235 
1236 
1237  // Doc solution
1238  doc_solution(doc_info);
1239 
1240  // Increment counter for solutions
1241  doc_info.number()++;
1242  } // End of timestepping loop
1243 
1244  // write footer and close pvd file
1245  ParaviewHelper::write_pvd_footer(Global_Physical_Variables::Pvd_file);
1247 
1248 } // End of unsteady_run
1249 
1250 
1251 ///////////////////////////////////////////////////////////////////////
1252 ///////////////////////////////////////////////////////////////////////
1253 ///////////////////////////////////////////////////////////////////////
1254 
1255 
1256 //==start_of_main=========================================================
1257 /// Driver code for single fluid axisymmetric horizontal interface problem
1258 //========================================================================
1259 int main(int argc, char* argv[])
1260 {
1261 
1262  // Store command line arguments
1263  CommandLineArgs::setup(argc,argv);
1264 
1265  /// Maximum time
1266  double t_max = 1000.0;
1267 
1268  /// Duration of timestep
1269  const double dt = 0.1;
1270 
1271  // If we are doing validation run, use smaller number of timesteps
1272  if(CommandLineArgs::Argc>1)
1273  {
1274  t_max = 0.5;
1275  }
1276 
1277  // Number of elements in radial (r) direction
1278  const unsigned n_r = 10;
1279 
1280  // Number of elements in axial (z) direction
1281  const unsigned n_z = 80;
1282 
1283  // Height of domain
1284  const double l_z = MathematicalConstants::Pi/Global_Physical_Variables::Alpha;
1285 
1286  // Set direction of gravity (vertically downwards)
1290 
1291  // Set up the spine test problem with AxisymmetricQCrouzeixRaviartElements,
1292  // using the BDF<2> timestepper
1294  problem(n_r,n_z,l_z);
1295 
1296  // Run the unsteady simulation
1297  problem.unsteady_run(t_max,dt);
1298 } // End of main
1299 
double Peclet_S
Surface Peclet number.
virtual void spine_node_update(SpineNode *spine_node_pt)
Node update function assumed spines rooted at the wall fixed to be at r=1 and directed inwards to r=0...
void fill_in_contribution_to_jacobian(Vector< double > &residuals, DenseMatrix< double > &jacobian)
Fill in the contribution to the residuals Calculate the contribution to the jacobian.
double Peclet_St_S
Sufrace Peclet number multiplied by Strouhal number.
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...
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...
static double Default_Physical_Constant_Value
Default value of the physical constants.
ofstream Pvd_file
Pvd file – a wrapper for all the different vtu output files plus information about continuous time to...
double ReSt
Womersley = Reynolds times Strouhal.
double interpolated_C(const Vector< double > &s)
Get the surfactant concentration.
double Epsilon
Free surface cosine deformation parameter.
double integrate_c() const
Compute the concentration intergated over the area.
double ReInvFr
Product of Reynolds and Froude number.
int main(int argc, char *argv[])
Driver code for single fluid axisymmetric horizontal interface problem.
Vector< double > G(3)
Direction of gravity.
MyHorizontalSingleLayerSpineMesh(const unsigned &nx, const unsigned &ny, const double &lx, const double &ly, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass number of elements in x-direction, number of elements in y-direction, radial extent, axial length , and pointer to timestepper (defaults to Steady timestepper)
void fill_in_contribution_to_jacobian_and_mass_matrix(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix)
MyHorizontalSingleLayerSpineMesh< ELEMENT > * Bulk_mesh_pt
Access function for the specific mesh.
double *& peclet_s_pt()
Access function for pointer to the surface Peclet number.
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.
double *& beta_pt()
Access function for pointer to the Elasticity number.
void doc_solution(DocInfo &doc_info)
Doc the solution.
Mesh * Interface_mesh_pt
Mesh for the free surface (interface) elements.
double dcdt_surface(const unsigned &l) const
The time derivative of the surface concentration.
SpineAxisymmetricMarangoniSurfactantFluidInterfaceElement(FiniteElement *const &element_pt, const int &face_index)
double Alpha
Wavelength of the domain.
double *& peclet_strouhal_s_pt()
Access function for pointer to the surface Peclet x Strouhal number.
double compute_total_mass()
Compute the total mass of the insoluble surfactant.
double *& ca_pt()
Pointer to the Capillary number.