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-2022 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
31namespace 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
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 //=====================================================================
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
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
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
TimeStepper *& time_stepper_pt()
Return the pointer to the timestepper.
Definition: nodes.h:238
void outer_unit_normal(const Vector< double > &s, Vector< double > &unit_normal) const
Compute outer unit normal at the specified local coordinate.
Definition: elements.cc:6006
double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s. Overloaded to get information from bulk...
Definition: elements.h:4528
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
virtual void dshape_local_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsids) const
Return the geometric shape function and its derivative w.r.t. the local coordinates at the ipt-th int...
Definition: elements.cc:3239
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
Definition: elements.h:2593
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3161
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node.
Definition: elements.h:1432
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
Definition: elements.h:2611
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:3148
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3186
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
Definition: elements.h:2317
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
unsigned nodal_dimension() const
Return the required Eulerian dimension of the nodes in this element.
Definition: elements.h:2484
double dnodal_position_dt(const unsigned &n, const unsigned &i) const
Return the i-th component of nodal velocity: dx/dt at local node n.
Definition: elements.h:2333
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,...
const double & ca() const
The value of the Capillary number.
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.
TimeStepper *& time_stepper_pt()
Access function for pointer to time stepper: Null if object is not time-dependent.
Definition: geom_objects.h:192
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:1054
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
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.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
unsigned ntstorage() const
Return the number of doubles required to represent history (one for steady)
Definition: timesteppers.h:601
virtual double weight(const unsigned &i, const unsigned &j) const
Access function for j-th weight for the i-th derivative.
Definition: timesteppers.h:594
std::string type() const
Return string that indicates the type of the timestepper (e.g. "BDF", "Newmark", etc....
Definition: timesteppers.h:490
//////////////////////////////////////////////////////////////////// ////////////////////////////////...