surfactant_transport_elements.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: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
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 //Non-inline functions for fluid free surface elements
31 
32 //OOMPH-LIB headers
34 
35 namespace oomph
36 {
37 
38 
39 //Define the default physical value to be one
41 
42  //=====================================================================
43  ///Get the surfactant concentration
44  //=====================================================================
46  {
47  //Find number of nodes
48  unsigned n_node = this->nnode();
49 
50  //Local shape function
51  Shape psi(n_node);
52 
53  //Find values of shape function
54  this->shape(s,psi);
55 
56  //Initialise value of C
57  double C = 0.0;
58 
59  //Loop over the local nodes and sum
60  for(unsigned l=0;l<n_node;l++)
61  {
62  C += this->nodal_value(l,this->C_index[l])*psi(l);
63  }
64 
65  return(C);
66  }
67 
68  //=====================================================================
69  /// The time derivative of the surface concentration
70  //=====================================================================
71  double SurfactantTransportInterfaceElement::dcdt_surface(const unsigned &l) const
72  {
73  // Get the data's timestepper
74  TimeStepper* time_stepper_pt= this->node_pt(l)->time_stepper_pt();
75 
76  //Initialise dudt
77  double dcdt=0.0;
78  //Loop over the timesteps, if there is a non Steady timestepper
79  if (time_stepper_pt->type()!="Steady")
80  {
81  // Number of timsteps (past & present)
82  const unsigned n_time = time_stepper_pt->ntstorage();
83 
84  for(unsigned t=0;t<n_time;t++)
85  {
86  dcdt += time_stepper_pt->weight(1,t)*this->nodal_value(t,l,this->C_index[l]);
87  }
88  }
89  return dcdt;
90  }
91 
92  //=====================================================================
93  /// The surface tension function is linear in the
94  /// concentration with constant of proportionality equal
95  /// to the elasticity number.
96  //=====================================================================
97  double SurfactantTransportInterfaceElement::sigma(const Vector<double> &s)
98  {
99  //Find the number of shape functions
100  const unsigned n_node = this->nnode();
101  //Now get the shape fuctions at the local coordinate
102  Shape psi(n_node);
103  this->shape(s,psi);
104 
105  //Now interpolate the temperature and surfactant concentration
106  double C=0.0;
107  for(unsigned l=0;l<n_node;l++)
108  {
109  C += this->nodal_value(l,this->C_index[l])*psi(l);
110  }
111 
112  //Get the Elasticity numbers
113  double Beta = this->beta();
114  //Return the variable surface tension
115  return (1.0 - Beta*(C-1.0));
116  }
117 
118  //=======================================================================
119  /// \short Overload the Helper function to calculate the residuals and
120  /// jacobian entries. This particular function ensures that the
121  /// additional entries are calculated inside the integration loop
123  Vector<double> &residuals, DenseMatrix<double> &jacobian,
124  const unsigned &flag,const Shape &psif, const DShape &dpsifds,
125  const DShape &dpsifdS, const DShape &dpsifdS_div,
126  const Vector<double> &s,
127  const Vector<double> &interpolated_x, const Vector<double> &interpolated_n,
128  const double &W,const double &J)
129  {
130  //Find out how many nodes there are
131  unsigned n_node = this->nnode();
132 
133  //Storage for the local equation numbers and unknowns
134  int local_eqn = 0, local_unknown = 0;
135 
136  //Surface advection-diffusion equation
137 
138  //Find the index at which the concentration is stored
139  Vector<unsigned> u_index = this->U_index_interface;
140 
141  //Read out the surface peclect number
142  const double Pe_s = this->peclet_s();
143  //const double PeSt_s = this->peclet_strouhal_s();
144 
145  //Now calculate the concentration and derivatives at this point
146  //Assuming the same shape functions are used (which they are)
147  double interpolated_C = 0.0;
148  double dCdt = 0.0;
149  //The tangent vectors and velocity
150  const unsigned n_dim = this->node_pt(0)->ndim();
151  Vector<double> interpolated_u(n_dim,0.0);
152  Vector<double> mesh_velocity(n_dim,0.0);
153  Vector<double> interpolated_grad_C(n_dim,0.0);
154  double interpolated_div_u = 0.0;
155 
156  //Loop over the shape functions
157  for(unsigned l=0;l<n_node;l++)
158  {
159  const double psi = psif(l);
160  const double C_ = this->nodal_value(l,this->C_index[l]);
161 
162  interpolated_C += C_*psi;
163  dCdt += dcdt_surface(l)*psi;
164  //Velocity and Mesh Velocity
165  for(unsigned j=0;j<n_dim;j++)
166  {
167  const double u_ = this->nodal_value(l,u_index[j]);
168  interpolated_u[j] += u_*psi;
169  mesh_velocity[j] += this->dnodal_position_dt(l,j)*psi;
170  interpolated_grad_C[j] += C_*dpsifdS(l,j);
171  interpolated_div_u += u_*dpsifdS_div(l,j);
172  }
173  }
174 
175  //Pre-compute advection term
176  double interpolated_advection_term = interpolated_C*interpolated_div_u;
177  for(unsigned i=0;i<n_dim;i++)
178  {
179  interpolated_advection_term += (interpolated_u[i] - mesh_velocity[i])*interpolated_grad_C[i];
180  }
181 
182 
183  //Now we add the flux term to the appropriate residuals
184  for(unsigned l=0;l<n_node;l++)
185  {
186  //Read out the apprporiate local equation
187  local_eqn = this->nodal_local_eqn(l,this->C_index[l]);
188 
189  //If not a boundary condition
190  if(local_eqn >= 0)
191  {
192  //Time derivative term
193  residuals[local_eqn] += dCdt*psif(l)*W*J;
194 
195  //First Advection term
196  residuals[local_eqn] += interpolated_advection_term*psif(l)*W*J;
197 
198  //Diffusion term
199  double diffusion_term = 0.0;
200  for(unsigned i=0;i<n_dim;i++)
201  {
202  diffusion_term += interpolated_grad_C[i]*dpsifdS(l,i);
203  }
204  residuals[local_eqn] += (1.0/Pe_s)*diffusion_term*W*J;
205 
206  //We also need to worry about the jacobian terms
207  if(flag)
208  {
209  //Loop over the nodes again
210  for(unsigned l2=0;l2<n_node;l2++)
211  {
212  //Get the time stepper
213  TimeStepper* time_stepper_pt=this->node_pt(l2)->time_stepper_pt();
214 
215  //Get the unknown c_index
216  local_unknown =this->nodal_local_eqn(l2,this->C_index[l2]);
217 
218  if(local_unknown >=0)
219  {
220  jacobian(local_eqn,local_unknown) +=
221  time_stepper_pt->weight(1,0)*psif(l2)*psif(l)*W*J;
222 
223  jacobian(local_eqn,local_unknown) +=
224  psif(l2)*interpolated_div_u*psif(l)*W*J;
225 
226  for(unsigned i=0;i<n_dim;i++)
227  {
228  jacobian(local_eqn,local_unknown) +=
229  (interpolated_u[i] - mesh_velocity[i])*dpsifdS(l2,i)*psif(l)*W*J;
230  }
231 
232  for(unsigned i=0;i<n_dim;i++)
233  {
234  jacobian(local_eqn,local_unknown) += (1.0/Pe_s)*dpsifdS(l2,i)*dpsifdS(l,i)*W*J;
235  }
236  }
237 
238 
239  //Loop over the velocity components
240  for(unsigned i2=0;i2<n_dim;i2++)
241  {
242 
243  //Get the unknown
244  local_unknown = this->nodal_local_eqn(l2,u_index[i2]);
245 
246 
247  //If not a boundary condition
248  if(local_unknown >= 0)
249  {
250  //Bits from the advection term
251  jacobian(local_eqn,local_unknown) += (interpolated_C*dpsifdS_div(l2,i2) +
252  psif(l2)*interpolated_grad_C[i2])
253  *psif(l)*W*J;
254  }
255  }
256  }
257  }
258  }
259 
260  if(flag)
261  {
262  const double dsigma = this->dsigma_dC(s);
263  const double Ca = this->ca();
264  for(unsigned l2=0;l2<n_node;l2++)
265  {
266  local_unknown = this->nodal_local_eqn(l2,this->C_index[l2]);
267  if(local_unknown >= 0)
268  {
269  const double psi_ = psif(l2);
270  for(unsigned i=0;i<n_dim;i++)
271  {
272  //Add the Jacobian contribution from the surface tension
273  local_eqn = this->nodal_local_eqn(l,u_index[i]);
274  if(local_eqn >= 0)
275  {
276  jacobian(local_eqn,local_unknown) -= (dsigma/Ca)*psi_*dpsifdS_div(l,i)*J*W;
277  }
278  }
279  }
280  }
281  }
282 
283  } //End of loop over the nodes
284 
285  }
286 
287  //=======================================================
288  ///Overload the output function
289  //=======================================================
291  std::ostream &outfile, const unsigned &n_plot)
292  {
293  outfile.precision(16);
294 
295  const unsigned el_dim = this->dim();
296  const unsigned n_dim = this->nodal_dimension();
297  const unsigned n_velocity = this->U_index_interface.size();
298 
299  //Set output Vector
300  Vector<double> s(el_dim);
301  Vector<double> n(n_dim);
302  Vector<double> u(n_velocity);
303 
304 
305  outfile << this->tecplot_zone_string(n_plot);
306 
307  // Loop over plot points
308  unsigned num_plot_points=this->nplot_points(n_plot);
309  for (unsigned iplot=0;iplot<num_plot_points;iplot++)
310  {
311  // Get local coordinates of plot point
312  this->get_s_plot(iplot,n_plot,s);
313  //Get the outer unit normal
314  this->outer_unit_normal(s,n);
315 
316  double u_n = 0.0;
317  for(unsigned i=0;i<n_velocity;i++)
318  {
319  u[i] = this->interpolated_u(s,i);
320  }
321 
322  //Not the same as above for axisymmetric case
323  for(unsigned i=0;i<n_dim;i++)
324  {
325  u_n += u[i]*n[i];
326  }
327 
328  //Output the x,y,u,v
329  for(unsigned i=0;i<n_dim;i++) outfile << this->interpolated_x(s,i) << " ";
330  for(unsigned i=0;i<n_dim;i++) outfile << u[i] << " ";
331  //Output a dummy pressure
332  outfile << 0.0 << " ";
333  //Output the concentration
334  outfile << interpolated_C(s) << " ";
335  //Output the interfacial tension
336  outfile << sigma(s) << " ";
337  for(unsigned i=0;i<n_dim;i++)
338  {outfile << u[i] - u_n*n[i] << " ";}
339  outfile << std::endl;
340  }
341  outfile << std::endl;
342  }
343 
344  //=======================================================================
345  /// \short Compute the concentration intergated over the surface area
346  //=======================================================================
348  {
349  //Find out how many nodes there are
350  unsigned n_node = this->nnode();
351 
352  const unsigned el_dim = this->dim();
353  const unsigned n_dim = this->nodal_dimension();
354 
355  //Set up memeory for the shape functions
356  Shape psif(n_node);
357  DShape dpsifds(n_node,el_dim);
358  DShape dpsifdS(n_node,n_dim);
359  DShape dpsifdS_div(n_node,n_dim);
360 
361  //Set the value of n_intpt
362  unsigned n_intpt = this->integral_pt()->nweight();
363 
364  //Storage for the local coordinate
365  Vector<double> s(el_dim);
366 
367  //Storage for the total mass
368  double mass = 0.0;
369 
370  //Loop over the integration points
371  for(unsigned ipt=0;ipt<n_intpt;ipt++)
372  {
373  //Get the local coordinate at the integration point
374  for(unsigned i=0;i<el_dim;i++) {s[i] = this->integral_pt()->knot(ipt,i);}
375 
376  //Get the integral weight
377  double W = this->integral_pt()->weight(ipt);
378 
379  //Call the derivatives of the shape function
380  this->dshape_local_at_knot(ipt,psif,dpsifds);
381 
382  Vector<double> interpolated_x(n_dim,0.0);
383  DenseMatrix<double> interpolated_t(el_dim,n_dim,0.0);
384  double interpolated_c = 0.0;
385 
386  //Loop over the shape functions
387  for(unsigned l=0;l<n_node;l++)
388  {
389  const double psi_ = psif(l);
390  interpolated_c += this->nodal_value(l,this->C_index[l])*psi_;
391  //Loop over directional components
392  for(unsigned i=0;i<n_dim;i++)
393  {
394  // Coordinate
395  interpolated_x[i] += this->nodal_position(l,i)*psi_;
396 
397  //Calculate the tangent vectors
398  for(unsigned j=0;j<el_dim;j++)
399  {
400  interpolated_t(j,i) += this->nodal_position(l,i)*dpsifds(l,j);
401  }
402  }
403  }
404 
405  //Calculate the surface gradient and divergence
406  double J =
407  this->compute_surface_derivatives(psif,dpsifds,
408  interpolated_t,
409  interpolated_x,
410  dpsifdS,
411  dpsifdS_div);
412 
413  mass += interpolated_c*W*J;
414  }
415  return mass;
416  }
417 
418 
419 }
double dcdt_surface(const unsigned &l) const
The time derivative of the surface concentration.
double peclet_s()
Return the surface peclect number.
double u(const unsigned &j, const unsigned &i)
Return the i-th velocity component at local node j.
const double & ca() const
The value of the Capillary number.
double integrate_c()
Compute the concentration intergated over the surface area.
double interpolated_C(const Vector< double > &s)
Get the surfactant concentration.
void output(std::ostream &outfile)
Overload the output function.
static double Default_Physical_Constant_Value
Default value of the physical constants.
double interpolated_u(const Vector< double > &s, const unsigned &i)
Calculate the i-th velocity component at the local coordinate s.
double Beta
Surface Elasticity number (weak case)
virtual double compute_surface_derivatives(const Shape &psi, const DShape &dpsids, const DenseMatrix< double > &interpolated_t, const Vector< double > &interpolated_x, DShape &dpsidS, DShape &dpsidS_div)=0
Compute the surface gradient and surface divergence operators given the shape functions, derivatives, tangent vectors and position. All derivatives and tangent vectors should be formed with respect to the local coordinates.
Vector< unsigned > U_index_interface
Nodal index at which the i-th velocity component is stored.
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...