33 #include "axisym_navier_stokes.h"
34 #include "navier_stokes.h"
37 #include "fluid_interface.h"
40 #include "oomph_crbond_bessel.h"
43 #include "meshes/two_layer_spine_mesh.h"
47 using namespace oomph;
93 template <
class ELEMENT,
class TIMESTEPPER>
104 const unsigned &n_z2,
const double &l_r,
105 const double &h1,
const double &h2);
130 Bulk_mesh_pt->node_update();
145 const unsigned &pdof,
146 const double &pvalue)
149 dynamic_cast<ELEMENT*
>(Bulk_mesh_pt->element_pt(e))->
150 fix_pressure(pdof,pvalue);
164 Mesh* Surface_mesh_pt;
174 template <
class ELEMENT,
class TIMESTEPPER>
177 const unsigned &n_z2,
const double &l_r,
178 const double& h1,
const double &h2) : Lr(l_r)
182 add_time_stepper_pt(
new TIMESTEPPER);
186 (n_r,n_z1,n_z2,l_r,h1,h2,time_stepper_pt());
191 for(
unsigned i=0;i<n_r;i++)
195 FiniteElement *interface_element_pt =
new
196 SpineAxisymmetricFluidInterfaceElement<ELEMENT>(
210 for(
unsigned b=0;b<6;b++)
213 const unsigned n_node =
Bulk_mesh_pt->nboundary_node(b);
216 for (
unsigned n=0;n<n_node;n++)
227 if(b==0 || b==3) {
Bulk_mesh_pt->boundary_node_pt(b,n)->pin(1); }
233 for(
unsigned n=0;n<n_node;n++) {
Bulk_mesh_pt->node_pt(n)->pin(2); }
244 for(
unsigned e=0;e<n_lower;e++)
248 lower_layer_element_pt(e));
266 for(
unsigned e=0;e<n_upper;e++)
270 upper_layer_element_pt(e));
300 for(
unsigned e=0;e<n_interface_element;e++)
303 SpineAxisymmetricFluidInterfaceElement<ELEMENT>* el_pt =
304 dynamic_cast<SpineAxisymmetricFluidInterfaceElement<ELEMENT>*
>
322 this->build_global_mesh();
326 cout <<
"Number of equations: " << assign_eqn_numbers() << std::endl;
337 template <
class ELEMENT,
class TIMESTEPPER>
341 const unsigned n_node = Bulk_mesh_pt->nnode();
344 for(
unsigned n=0;n<n_node;n++)
347 for(
unsigned i=0;i<3;i++)
350 Bulk_mesh_pt->node_pt(n)->set_value(i,0.0);
356 assign_initial_values_impulsive();
367 template <
class ELEMENT,
class TIMESTEPPER>
371 for(
unsigned b=0;b<6;b++)
374 const unsigned n_node = Bulk_mesh_pt->nboundary_node(b);
377 for(
unsigned n=0;n<n_node;n++)
380 Bulk_mesh_pt->boundary_node_pt(b,n)->set_value(0,0.0);
383 Bulk_mesh_pt->boundary_node_pt(b,n)->set_value(2,0.0);
388 Bulk_mesh_pt->boundary_node_pt(b,n)->set_value(1,0.0);
399 template <
class ELEMENT,
class TIMESTEPPER>
404 double j0, j1, y0, y1, j0p, j1p, y0p, y1p;
407 const unsigned n_spine = Bulk_mesh_pt->nspine();
410 for(
unsigned i=0;i<n_spine;i++)
413 const double r_value = Bulk_mesh_pt->boundary_node_pt(0,i)->x(0);
416 CRBond_Bessel::bessjy01a(k*r_value,j0,j1,y0,y1,j0p,j1p,y0p,y1p);
419 Bulk_mesh_pt->spine_pt(i)->height() = 1.0 + epsilon*j0;
423 Bulk_mesh_pt->node_update();
432 template <
class ELEMENT,
class TIMESTEPPER>
437 cout <<
"Time is now " << time_pt()->time() << std::endl;
441 Trace_file << time_pt()->time() <<
" "
442 << Bulk_mesh_pt->spine_pt(0)->height() <<
" " << std::endl;
448 const unsigned npts = 5;
451 sprintf(filename,
"%s/soln%i.dat",doc_info.directory().c_str(),
453 some_file.open(filename);
456 Bulk_mesh_pt->output(some_file,npts);
457 Surface_mesh_pt->output(some_file,npts);
469 template <
class ELEMENT,
class TIMESTEPPER>
475 const double epsilon = 0.1;
478 const double k_bessel = 3.8317;
481 deform_free_surface(epsilon,k_bessel);
487 doc_info.set_directory(
"RESLT");
494 sprintf(filename,
"%s/trace.dat",doc_info.directory().c_str());
495 Trace_file.open(filename);
498 Trace_file <<
"time, interface height" << std::endl;
504 set_initial_condition();
507 const unsigned n_timestep = unsigned(t_max/dt);
510 doc_solution(doc_info);
516 for(
unsigned t=1;t<=n_timestep;t++)
519 cout <<
"\nTimestep " << t <<
" of " << n_timestep << std::endl;
522 unsteady_newton_solve(dt);
525 doc_solution(doc_info);
543 int main(
int argc,
char* argv[])
546 CommandLineArgs::setup(argc,argv);
556 const double dt = 0.005;
559 if(CommandLineArgs::Argc>1) { t_max = 0.01; }
562 const unsigned n_r = 12;
565 const unsigned n_z1 = 12;
568 const unsigned n_z2 = 12;
571 const double l_r = 1.0;
574 const double h1 = 1.0;
577 const double h2 = 1.0;
588 problem(n_r,n_z1,n_z2,l_r,h1,h2);
/////////////////////////////////////////////////////////////////////// /////////////////////////////...
Mesh * Surface_mesh_pt
Pointer to the "surface" mesh.
void set_initial_condition()
Set initial conditions.
TwoLayerSpineMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the specific bulk mesh.
void doc_solution(DocInfo &doc_info)
Doc the solution.
InterfaceProblem()
Constructor.
void set_boundary_conditions()
Set boundary conditions.
~InterfaceProblem()
Destructor (empty)
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.
void deform_free_surface(const double &epsilon, const double &k)
Deform the mesh/free surface to a prescribed function.
void actions_before_newton_convergence_check()
Spine heights/lengths are unknowns in the problem so their values get corrected during each Newton st...
ElasticRefineableTwoLayerMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the (specific) "bulk" mesh.
void actions_before_newton_solve()
No actions required before solve step.
void unsteady_run(const double &t_max, const double &dt)
Do unsteady run up to maximum time t_max with given timestep dt.
void actions_after_newton_solve()
Update after solve can remain empty, because the update is performed automatically after every Newton...
Namespace for physical parameters.
double ReSt
Womersley number (Reynolds x Strouhal)
double St
Strouhal number.
double Density_Ratio
Ratio of density in upper fluid to density in lower fluid. Reynolds number etc. is based on density i...
double Ca
Capillary number.
double ReInvFr
Product of Reynolds number and inverse of Froude number.
double Re
Reynolds number.
double Viscosity_Ratio
Ratio of viscosity in upper fluid to viscosity in lower fluid. Reynolds number etc....
Vector< double > G(3)
Direction of gravity.
int main(int argc, char *argv[])
/////////////////////////////////////////////////////////////////////// /////////////////////////////...