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