spine_two_layer_interface.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 // Driver for two-dimensional two fluid interface problem, where
31 // the mesh is deformed using a spine-based node-update strategy
32 
33 // Generic oomph-lib header
34 #include "generic.h"
35 
36 // Navier-Stokes headers
37 #include "navier_stokes.h"
38 
39 // Interface headers
40 #include "fluid_interface.h"
41 
42 // The mesh
43 #include "meshes/two_layer_spine_mesh.h"
44 
45 using namespace std;
46 
47 using namespace oomph;
48 
49 
50 //==start_of_namespace====================================================
51 /// Namespace for physical parameters
52 //========================================================================
53 namespace Global_Physical_Variables
54 {
55 
56  /// Reynolds number
57  double Re = 5.0;
58 
59  /// Strouhal number
60  double St = 1.0;
61 
62  /// Womersley number (Reynolds x Strouhal, computed automatically)
63  double ReSt;
64 
65  /// Product of Reynolds number and inverse of Froude number
66  double ReInvFr = 5.0; // (Fr = 1)
67 
68  /// \short Ratio of viscosity in upper fluid to viscosity in lower
69  /// fluid. Reynolds number etc. is based on viscosity in lower fluid.
70  double Viscosity_Ratio = 0.1;
71 
72  /// \short Ratio of density in upper fluid to density in lower
73  /// fluid. Reynolds number etc. is based on density in lower fluid.
74  double Density_Ratio = 0.5;
75 
76  /// Capillary number
77  double Ca = 0.01;
78 
79  /// Direction of gravity
80  Vector<double> G(2);
81 
82 } // End of namespace
83 
84 
85 //////////////////////////////////////////////////////////////////////////
86 //////////////////////////////////////////////////////////////////////////
87 //////////////////////////////////////////////////////////////////////////
88 
89 
90 //==start_of_problem_class================================================
91 /// Two fluid interface problem in a rectangular domain which is
92 /// periodic in the x direction
93 //========================================================================
94 template <class ELEMENT, class TIMESTEPPER>
95 class InterfaceProblem : public Problem
96 {
97 
98 public:
99 
100  /// Constructor: Pass the number of elements and the width of the
101  /// domain in the x direction. Also pass the number of elements in
102  /// the y direction of the bottom (fluid 1) and top (fluid 2) layers,
103  /// along with the heights of both layers.
104  InterfaceProblem(const unsigned &n_x, const unsigned &n_y1,
105  const unsigned &n_y2, const double &l_x,
106  const double &h1, const double &h2);
107 
108  /// Destructor (empty)
110 
111  /// Set initial conditions
112  void set_initial_condition();
113 
114  /// Set boundary conditions
115  void set_boundary_conditions();
116 
117  /// Doc the solution
118  void doc_solution(DocInfo &doc_info);
119 
120  /// Do unsteady run up to maximum time t_max with given timestep dt
121  void unsteady_run(const double &t_max, const double &dt);
122 
123 private:
124 
125  /// \short Spine heights/lengths are unknowns in the problem so their
126  /// values get corrected during each Newton step. However, changing
127  /// their value does not automatically change the nodal positions, so
128  /// we need to update all of them here.
130  {
131  Bulk_mesh_pt->node_update();
132  }
133 
134  /// No actions required before solve step
136 
137  /// \short Update after solve can remain empty, because the update
138  /// is performed automatically after every Newton step.
140 
141  /// Deform the mesh/free surface to a prescribed function
142  void deform_free_surface(const double &epsilon, const unsigned &n_periods);
143 
144  /// Fix pressure in element e at pressure dof pdof and set to pvalue
145  void fix_pressure(const unsigned &e,
146  const unsigned &pdof,
147  const double &pvalue)
148  {
149  // Fix the pressure at that element
150  dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e))->
151  fix_pressure(pdof,pvalue);
152  }
153 
154  /// Width of domain
155  double Lx;
156 
157  /// Trace file
158  ofstream Trace_file;
159 
160  /// \short Pointer to the specific bulk mesh
161  TwoLayerSpineMesh<ELEMENT>* Bulk_mesh_pt;
162 
163  /// \short Pointer to the surface mesh
164  Mesh* Surface_mesh_pt;
165 
166 }; // End of problem class
167 
168 
169 
170 //==start_of_constructor==================================================
171 /// Constructor for two fluid interface problem
172 //========================================================================
173 template <class ELEMENT, class TIMESTEPPER>
175 InterfaceProblem(const unsigned &n_x, const unsigned &n_y1,
176  const unsigned &n_y2, const double &l_x,
177  const double& h1, const double &h2) : Lx(l_x)
178 {
179 
180  // Allocate the timestepper (this constructs the time object as well)
181  add_time_stepper_pt(new TIMESTEPPER);
182 
183  // Build and assign mesh (the "true" boolean flag tells the mesh
184  // constructor that the domain is periodic in x)
185  Bulk_mesh_pt = new TwoLayerSpineMesh<ELEMENT>
186  (n_x,n_y1,n_y2,l_x,h1,h2,true,time_stepper_pt());
187 
188  Surface_mesh_pt = new Mesh;
189 
190  //Loop over the horizontal elements
191  for(unsigned i=0;i<n_x;i++)
192  {
193  //Construct a new 1D line element on the face on which the local
194  //coordinate 1 is fixed at its max. value (1) -- Face 2
195  FiniteElement *interface_element_pt = new
196  SpineLineFluidInterfaceElement<ELEMENT>(
197  Bulk_mesh_pt->finite_element_pt(n_x*(n_y1-1)+i),2);
198 
199  this->Surface_mesh_pt->add_element_pt(interface_element_pt);
200  }
201 
202 
203  // --------------------------------------------
204  // Set the boundary conditions for this problem
205  // --------------------------------------------
206 
207  // All nodes are free by default -- just pin the ones that have
208  // Dirichlet conditions here
209 
210  // Loop over mesh boundaries
211  for(unsigned b=0;b<5;b++)
212  {
213  // Determine number of nodes on boundary b
214  const unsigned n_node = Bulk_mesh_pt->nboundary_node(b);
215 
216  // Loop over nodes on boundary b
217  for (unsigned n=0;n<n_node;n++)
218  {
219  // Pin x-component of velocity on all boundaries (no slip/penetration)
220  Bulk_mesh_pt->boundary_node_pt(b,n)->pin(0);
221 
222  // Pin y-component of velocity on both solid boundaries (no penetration)
223  if(b==0 || b==3) { Bulk_mesh_pt->boundary_node_pt(b,n)->pin(1); }
224  }
225  }
226 
227  // ----------------------------------------------------------------
228  // Complete the problem setup to make the elements fully functional
229  // ----------------------------------------------------------------
230 
231  // Determine number of bulk elements in lower/upper fluids
232  const unsigned n_lower = Bulk_mesh_pt->nlower();
233  const unsigned n_upper = Bulk_mesh_pt->nupper();
234 
235  // Loop over bulk elements in lower fluid
236  for(unsigned e=0;e<n_lower;e++)
237  {
238  // Upcast from GeneralisedElement to the present element
239  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->
240  lower_layer_element_pt(e));
241 
242  // Set the Reynolds number
243  el_pt->re_pt() = &Global_Physical_Variables::Re;
244 
245  // Set the Womersley number
246  el_pt->re_st_pt() = &Global_Physical_Variables::ReSt;
247 
248  // Set the product of the Reynolds number and the inverse of the
249  // Froude number
250  el_pt->re_invfr_pt() = &Global_Physical_Variables::ReInvFr;
251 
252  // Set the direction of gravity
253  el_pt->g_pt() = &Global_Physical_Variables::G;
254 
255  } // End of loop over bulk elements in lower fluid
256 
257  // Loop over bulk elements in upper fluid
258  for(unsigned e=0;e<n_upper;e++)
259  {
260  // Upcast from GeneralisedElement to the present element
261  ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->
262  upper_layer_element_pt(e));
263 
264  // Set the Reynolds number
265  el_pt->re_pt() = &Global_Physical_Variables::Re;
266 
267  // Set the Womersley number
268  el_pt->re_st_pt() = &Global_Physical_Variables::ReSt;
269 
270  // Set the product of the Reynolds number and the inverse of the
271  // Froude number
272  el_pt->re_invfr_pt() = &Global_Physical_Variables::ReInvFr;
273 
274  // Set the viscosity ratio
275  el_pt->viscosity_ratio_pt() = &Global_Physical_Variables::Viscosity_Ratio;
276 
277  // Set the density ratio
278  el_pt->density_ratio_pt() = &Global_Physical_Variables::Density_Ratio;
279 
280  // Set the direction of gravity
281  el_pt->g_pt() = &Global_Physical_Variables::G;
282 
283  } // End of loop over bulk elements in upper fluid
284 
285  // Set the pressure in the first element at 'node' 0 to 0.0
286  fix_pressure(0,0,0.0);
287 
288  // Determine number of 1D interface elements in mesh
289  unsigned n_interface_element = Surface_mesh_pt->nelement();
290 
291  // Loop over interface elements
292  for(unsigned e=0;e<n_interface_element;e++)
293  {
294  // Upcast from GeneralisedElement to the present element
295  SpineLineFluidInterfaceElement<ELEMENT>* el_pt =
296  dynamic_cast<SpineLineFluidInterfaceElement<ELEMENT>*>
297  (Surface_mesh_pt->element_pt(e));
298 
299  // Set the Strouhal number
300  el_pt->ca_pt() = &Global_Physical_Variables::St;
301 
302  // Set the Capillary number
303  el_pt->ca_pt() = &Global_Physical_Variables::Ca;
304 
305  } // End of loop over interface elements
306 
307  // Apply the boundary conditions
309 
310  this->add_sub_mesh(Bulk_mesh_pt);
311  this->add_sub_mesh(Surface_mesh_pt);
312 
313  this->build_global_mesh();
314 
315 
316  // Setup equation numbering scheme
317  cout << "Number of equations: " << assign_eqn_numbers() << std::endl;
318 
319 } // End of constructor
320 
321 
322 
323 //==start_of_set_initial_condition========================================
324 /// \short Set initial conditions: Set all nodal velocities to zero and
325 /// initialise the previous velocities and nodal positions to correspond
326 /// to an impulsive start
327 //========================================================================
328 template <class ELEMENT, class TIMESTEPPER>
330 {
331  // Determine number of nodes in mesh
332  const unsigned n_node = Bulk_mesh_pt->nnode();
333 
334  // Loop over all nodes in mesh
335  for(unsigned n=0;n<n_node;n++)
336  {
337  // Loop over the two velocity components
338  for(unsigned i=0;i<2;i++)
339  {
340  // Set velocity component i of node n to zero
341  Bulk_mesh_pt->node_pt(n)->set_value(i,0.0);
342  }
343  }
344 
345  // Initialise the previous velocity values for timestepping
346  // corresponding to an impulsive start
347  assign_initial_values_impulsive();
348 
349 } // End of set_initial_condition
350 
351 
352 
353 //==start_of_set_boundary_conditions======================================
354 /// \short Set boundary conditions: Set both velocity components
355 /// to zero on the top and bottom (solid) walls and the horizontal
356 /// component only to zero on the side (periodic) boundaries
357 //========================================================================
358 template <class ELEMENT, class TIMESTEPPER>
360 {
361 
362  // Loop over mesh boundaries
363  for(unsigned b=0;b<5;b++)
364  {
365  // Determine number of nodes on boundary b
366  const unsigned n_node = Bulk_mesh_pt->nboundary_node(b);
367 
368  // Loop over nodes on boundary b
369  for(unsigned n=0;n<n_node;n++)
370  {
371  // Set x-component of the velocity to zero on all boundaries
372  Bulk_mesh_pt->boundary_node_pt(b,n)->set_value(0,0.0);
373 
374  // Set y-component of the velocity to zero on both solid
375  // boundaries (top and bottom)
376  if(b==0 || b==3)
377  {
378  Bulk_mesh_pt->boundary_node_pt(b,n)->set_value(1,0.0);
379  }
380  } // End of loop over nodes on boundary b
381  } // End of loop over mesh boundaries
382 
383 } // End of set_boundary_conditions
384 
385 
386 
387 //==start_of_deform_free_surface==========================================
388 /// Deform the mesh/free surface to a prescribed function
389 //========================================================================
390 template <class ELEMENT, class TIMESTEPPER>
392 deform_free_surface(const double &epsilon,const unsigned &n_periods)
393 {
394  // Determine number of spines in mesh
395  const unsigned n_spine = Bulk_mesh_pt->nspine();
396 
397  // Loop over spines in mesh
398  for(unsigned i=0;i<n_spine;i++)
399  {
400  // Determine x coordinate of spine
401  double x_value = Bulk_mesh_pt->boundary_node_pt(0,i)->x(0);
402 
403  // Set spine height
404  Bulk_mesh_pt->spine_pt(i)->height() =
405  1.0 + epsilon*(cos(2.0*n_periods*MathematicalConstants::Pi*x_value/Lx));
406  }
407 
408  // Update nodes in bulk mesh
409  Bulk_mesh_pt->node_update();
410 
411 } // End of deform_free_surface
412 
413 
414 
415 //==start_of_doc_solution=================================================
416 /// Doc the solution
417 //========================================================================
418 template <class ELEMENT, class TIMESTEPPER>
420 {
421 
422  // Output the time
423  cout << "Time is now " << time_pt()->time() << std::endl;
424 
425  // Document time and vertical position of left hand side of interface
426  // in trace file
427  Trace_file << time_pt()->time() << " "
428  << Bulk_mesh_pt->spine_pt(0)->height() << " " << std::endl;
429 
430  ofstream some_file;
431  char filename[100];
432 
433  // Set number of plot points (in each coordinate direction)
434  const unsigned npts = 5;
435 
436  // Open solution output file
437  sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
438  doc_info.number());
439  some_file.open(filename);
440 
441  // Output solution to file
442  Bulk_mesh_pt->output(some_file,npts);
443  Surface_mesh_pt->output(some_file,npts);
444 
445  // Close solution output file
446  some_file.close();
447 
448 } // End of doc_solution
449 
450 
451 
452 //==start_of_unsteady_run=================================================
453 /// Perform run up to specified time t_max with given timestep dt
454 //========================================================================
455 template <class ELEMENT, class TIMESTEPPER>
457 unsteady_run(const double &t_max, const double &dt)
458 {
459 
460  // Set value of epsilon
461  const double epsilon = 0.1;
462 
463  // Set number of periods for cosine term
464  const unsigned n_periods = 1;
465 
466  // Deform the mesh/free surface
467  deform_free_surface(epsilon,n_periods);
468 
469  // Initialise DocInfo object
470  DocInfo doc_info;
471 
472  // Set output directory
473  doc_info.set_directory("RESLT");
474 
475  // Initialise counter for solutions
476  doc_info.number()=0;
477 
478  // Open trace file
479  char filename[100];
480  sprintf(filename,"%s/trace.dat",doc_info.directory().c_str());
481  Trace_file.open(filename);
482 
483  // Initialise trace file
484  Trace_file << "time, free surface height" << std::endl;
485 
486  // Initialise timestep
487  initialise_dt(dt);
488 
489  // Set initial conditions
490  set_initial_condition();
491 
492  // Determine number of timesteps
493  const unsigned n_timestep = unsigned(t_max/dt);
494 
495  // Doc initial solution
496  doc_solution(doc_info);
497 
498  // Increment counter for solutions
499  doc_info.number()++;
500 
501  // Timestepping loop
502  for(unsigned t=1;t<=n_timestep;t++)
503  {
504  // Output current timestep to screen
505  cout << "\nTimestep " << t << " of " << n_timestep << std::endl;
506 
507  // Take one fixed timestep
508  unsteady_newton_solve(dt);
509 
510  // Doc solution
511  doc_solution(doc_info);
512 
513  // Increment counter for solutions
514  doc_info.number()++;
515 
516  } // End of timestepping loop
517 
518 } // End of unsteady_run
519 
520 
521 //////////////////////////////////////////////////////////////////////////
522 //////////////////////////////////////////////////////////////////////////
523 //////////////////////////////////////////////////////////////////////////
524 
525 
526 //==start_of_main======================================================
527 /// Driver code for two-dimensional two fluid interface problem
528 //=====================================================================
529 int main(int argc, char* argv[])
530 {
531  // Store command line arguments
532  CommandLineArgs::setup(argc,argv);
533 
534  // Compute the Womersley number
537 
538  /// Maximum time
539  double t_max = 0.6;
540 
541  /// Duration of timestep
542  const double dt = 0.0025;
543 
544  // If we are doing validation run, use smaller number of timesteps
545  if(CommandLineArgs::Argc>1) { t_max = 0.005; }
546 
547  // Number of elements in x direction
548  const unsigned n_x = 12;
549 
550  // Number of elements in y direction in lower fluid (fluid 1)
551  const unsigned n_y1 = 12;
552 
553  // Number of elements in y direction in upper fluid (fluid 2)
554  const unsigned n_y2 = 12;
555 
556  // Width of domain
557  const double l_x = 1.0;
558 
559  // Height of lower fluid layer
560  const double h1 = 1.0;
561 
562  // Height of upper fluid layer
563  const double h2 = 1.0;
564 
565  // Set direction of gravity (vertically downwards)
568 
569  // Set up the spine test problem with QCrouzeixRaviartElements,
570  // using the BDF<2> timestepper
572  problem(n_x,n_y1,n_y2,l_x,h1,h2);
573 
574  // Run the unsteady simulation
575  problem.unsteady_run(t_max,dt);
576 
577 } // End of main
TwoLayerSpineMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the specific bulk mesh.
int main(int argc, char *argv[])
Driver code for two-dimensional two fluid interface problem.
Vector< double > G(2)
Direction of gravity.
Mesh * Surface_mesh_pt
Pointer to the "surface" mesh.
void actions_before_newton_solve()
No actions required before solve step.
double ReSt
Womersley number (Reynolds x Strouhal)
void fix_pressure(const unsigned &e, const unsigned &pdof, const double &pvalue)
Fix pressure in element e at pressure dof pdof and set to pvalue.
~InterfaceProblem()
Destructor (empty)
double Density_Ratio
Ratio of density in upper fluid to density in lower fluid. Reynolds number etc. is based on density i...
void actions_after_newton_solve()
Update after solve can remain empty, because the update is performed automatically after every Newton...
void set_boundary_conditions()
Set boundary conditions.
void actions_before_newton_convergence_check()
Spine heights/lengths are unknowns in the problem so their values get corrected during each Newton st...
void set_initial_condition()
Set initial conditions.
void unsteady_run(const double &t_max, const double &dt)
Do unsteady run up to maximum time t_max with given timestep dt.
double Viscosity_Ratio
Ratio of viscosity in upper fluid to viscosity in lower fluid. Reynolds number etc. is based on viscosity in lower fluid.
void deform_free_surface(const double &epsilon, const unsigned &n_periods)
Deform the mesh/free surface to a prescribed function.
double ReInvFr
Product of Reynolds number and inverse of Froude number.
ElasticRefineableTwoLayerMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the (specific) "bulk" mesh.
void doc_solution(DocInfo &doc_info)
Document the solution.