30 #include "navier_stokes.h"
34 #include "meshes/one_d_lagrangian_mesh.h"
37 #include "meshes/collapsible_channel_mesh.h"
41 using namespace oomph;
110 void position(
const Vector<double>& zeta, Vector<double>& r)
const
121 void position(
const unsigned& t,
const Vector<double>& zeta,
122 Vector<double>& r)
const
137 DenseMatrix<double> &drdzeta,
138 RankThreeTensor<double> &ddrdzeta)
const
186 const Vector<double>& x,
187 const Vector<double>& n,
188 Vector<double>& traction)
208 void load(
const Vector<double>& xi,
const Vector<double>& x,
209 const Vector<double>& N, Vector<double>&
load)
211 for(
unsigned i=0;i<2;i++)
231 template <
class ELEMENT>
240 const unsigned& ncollapsible,
241 const unsigned& ndown,
244 const double& lcollapsible,
251 #ifdef MACRO_ELEMENT_NODE_UPDATE
254 MacroElementNodeUpdateRefineableCollapsibleChannelMesh<ELEMENT>*
bulk_mesh_pt()
259 MacroElementNodeUpdateRefineableCollapsibleChannelMesh<ELEMENT>*
>
271 RefineableAlgebraicCollapsibleChannelMesh<ELEMENT>*
>
287 void actions_before_adapt();
291 void actions_after_adapt();
304 Bulk_mesh_pt->node_update();
317 Mesh*
const &bulk_mesh_pt,
318 Mesh*
const &traction_mesh_pt);
321 void delete_traction_elements(Mesh*
const &traction_mesh_pt);
328 unsigned Ncollapsible;
348 #ifdef MACRO_ELEMENT_NODE_UPDATE
351 MacroElementNodeUpdateRefineableCollapsibleChannelMesh<ELEMENT>*
Bulk_mesh_pt;
362 Mesh* Applied_fluid_traction_mesh_pt;
365 OneDLagrangianMesh<FSIHermiteBeamElement>* Wall_mesh_pt;
384 template <
class ELEMENT>
387 const unsigned& ncollapsible,
388 const unsigned& ndown,
391 const double& lcollapsible,
397 Ncollapsible=ncollapsible;
401 Lcollapsible=lcollapsible;
407 Problem::Max_residuals=1000.0;
410 BDF<2>* fluid_time_stepper_pt=
new BDF<2>;
413 add_time_stepper_pt(fluid_time_stepper_pt);
416 Steady<2>* wall_time_stepper_pt =
new Steady<2>;
419 add_time_stepper_pt(wall_time_stepper_pt);
427 Wall_mesh_pt =
new OneDLagrangianMesh<FSIHermiteBeamElement>
428 (Ncollapsible,Lcollapsible,undeformed_wall_pt,wall_time_stepper_pt);
432 MeshAsGeomObject* wall_geom_object_pt=
433 new MeshAsGeomObject(Wall_mesh_pt);
435 #ifdef MACRO_ELEMENT_NODE_UPDATE
439 new MacroElementNodeUpdateRefineableCollapsibleChannelMesh<ELEMENT>
440 (nup, ncollapsible, ndown, ny,
441 lup, lcollapsible, ldown, ly,
443 fluid_time_stepper_pt);
449 Bulk_mesh_pt->node_update();
455 new RefineableAlgebraicCollapsibleChannelMesh<ELEMENT>
456 (nup, ncollapsible, ndown, ny,
457 lup, lcollapsible, ldown, ly,
460 fluid_time_stepper_pt);
464 if (CommandLineArgs::Argc>1)
467 Bulk_mesh_pt->refine_uniformly();
473 Applied_fluid_traction_mesh_pt =
new Mesh;
477 create_traction_elements(5,Bulk_mesh_pt,Applied_fluid_traction_mesh_pt);
480 add_sub_mesh(Bulk_mesh_pt);
481 add_sub_mesh(Applied_fluid_traction_mesh_pt);
482 add_sub_mesh(Wall_mesh_pt);
488 Z2ErrorEstimator* error_estimator_pt=
new Z2ErrorEstimator;
489 bulk_mesh_pt()->spatial_error_estimator_pt()=error_estimator_pt;
497 unsigned n_element=Bulk_mesh_pt->nelement();
498 for(
unsigned e=0;e<n_element;e++)
501 ELEMENT* el_pt =
dynamic_cast<ELEMENT*
>(Bulk_mesh_pt->element_pt(e));
513 RefineableNavierStokesEquations<2>::
514 pin_redundant_nodal_pressures(Bulk_mesh_pt->element_pt());
524 unsigned num_nod= bulk_mesh_pt()->nboundary_node(ibound);
525 for (
unsigned inod=0;inod<num_nod;inod++)
527 for(
unsigned i=0;i<2;i++)
529 bulk_mesh_pt()->boundary_node_pt(ibound, inod)->pin(i);
534 for(ibound=2;ibound<5;ibound++)
536 num_nod= bulk_mesh_pt()->nboundary_node(ibound);
537 for (
unsigned inod=0;inod<num_nod;inod++)
539 for(
unsigned i=0;i<2;i++)
541 bulk_mesh_pt()->boundary_node_pt(ibound, inod)->pin(i);
548 num_nod= bulk_mesh_pt()->nboundary_node(ibound);
549 for (
unsigned inod=0;inod<num_nod;inod++)
551 bulk_mesh_pt()->boundary_node_pt(ibound, inod)->pin(1);
557 num_nod= bulk_mesh_pt()->nboundary_node(ibound);
558 for (
unsigned inod=0;inod<num_nod;inod++)
561 bulk_mesh_pt()->boundary_node_pt(ibound, inod)->pin(1);
572 unsigned n_el=Applied_fluid_traction_mesh_pt->nelement();
573 for(
unsigned e=0;e<n_el;e++)
576 NavierStokesTractionElement<ELEMENT> *el_pt =
577 dynamic_cast< NavierStokesTractionElement<ELEMENT>*
>(
578 Applied_fluid_traction_mesh_pt->element_pt(e));
592 n_element = wall_mesh_pt()->nelement();
593 for(
unsigned e=0;e<n_element;e++)
596 FSIHermiteBeamElement *elem_pt =
597 dynamic_cast<FSIHermiteBeamElement*
>(wall_mesh_pt()->element_pt(e));
610 elem_pt->undeformed_beam_pt() = undeformed_wall_pt;
615 elem_pt->set_normal_pointing_out_of_fluid();
626 for(
unsigned b=0;b<2;b++)
629 wall_mesh_pt()->boundary_node_pt(b,0)->pin_position(0);
630 wall_mesh_pt()->boundary_node_pt(b,0)->pin_position(1);
640 num_nod= bulk_mesh_pt()->nboundary_node(ibound);
641 unsigned control_nod=num_nod/2;
642 Left_node_pt= bulk_mesh_pt()->boundary_node_pt(ibound, control_nod);
646 num_nod= bulk_mesh_pt()->nboundary_node(ibound);
647 control_nod=num_nod/2;
648 Right_node_pt= bulk_mesh_pt()->boundary_node_pt(ibound, control_nod);
652 num_nod= wall_mesh_pt()->nnode();
653 Wall_node_pt=wall_mesh_pt()->node_pt(Ncollapsible/2);
665 num_nod= bulk_mesh_pt()->nboundary_node(ibound);
666 for (
unsigned inod=0;inod<num_nod;inod++)
668 bulk_mesh_pt()->boundary_node_pt(ibound, inod)->
669 set_auxiliary_node_update_fct_pt(
670 FSI_functions::apply_no_slip_on_moving_wall);
679 FSI_functions::setup_fluid_load_info_for_solid_elements<ELEMENT,2>
680 (
this,3,Bulk_mesh_pt,Wall_mesh_pt);
683 cout <<
"Number of equations: " << assign_eqn_numbers() << std::endl;
694 template <
class ELEMENT>
696 ofstream& trace_file)
699 #ifdef MACRO_ELEMENT_NODE_UPDATE
702 if (CommandLineArgs::Argc>1)
704 FSI_functions::doc_fsi<MacroElementNodeUpdateNode>
705 (Bulk_mesh_pt,Wall_mesh_pt,doc_info);
711 if (CommandLineArgs::Argc>1)
713 FSI_functions::doc_fsi<AlgebraicNode>(Bulk_mesh_pt,Wall_mesh_pt,doc_info);
726 sprintf(filename,
"%s/soln%i.dat",doc_info.directory().c_str(),
728 some_file.open(filename);
729 bulk_mesh_pt()->output(some_file,npts);
733 sprintf(filename,
"%s/beam%i.dat",doc_info.directory().c_str(),
735 some_file.open(filename);
736 wall_mesh_pt()->output(some_file,npts);
742 unsigned nsteps=time_stepper_pt(1)->nprev_values();
743 for (
unsigned t=0;t<=nsteps;t++)
745 sprintf(filename,
"%s/wall%i-%i.dat",doc_info.directory().c_str(),
746 doc_info.number(),t);
747 some_file.open(filename);
748 unsigned n_elem=wall_mesh_pt()->nelement();
749 for (
unsigned ielem=0;ielem<n_elem;ielem++)
751 dynamic_cast<FSIHermiteBeamElement*
>(wall_mesh_pt()->element_pt(ielem))->
752 output(t,some_file,npts);
760 trace_file << time_pt()->time() <<
" "
761 << Wall_node_pt->x(1) <<
" "
762 << Left_node_pt->value(0) <<
" "
763 << Right_node_pt->value(0) <<
" "
776 template <
class ELEMENT>
778 const unsigned &b, Mesh*
const &bulk_mesh_pt, Mesh*
const &traction_mesh_pt)
782 unsigned n_element = bulk_mesh_pt->nboundary_element(b);
785 for(
unsigned e=0;e<n_element;e++)
788 ELEMENT* bulk_elem_pt =
dynamic_cast<ELEMENT*
>
789 (bulk_mesh_pt->boundary_element_pt(b,e));
792 int face_index = bulk_mesh_pt->face_index_at_boundary(b,e);
795 NavierStokesTractionElement<ELEMENT>* flux_element_pt =
796 new NavierStokesTractionElement<ELEMENT>(bulk_elem_pt,face_index);
799 traction_mesh_pt->add_element_pt(flux_element_pt);
810 template<
class ELEMENT>
815 unsigned n_element = surface_mesh_pt->nelement();
818 for(
unsigned e=0;e<n_element;e++)
821 delete surface_mesh_pt->element_pt(e);
825 surface_mesh_pt->flush_element_and_node_storage();
834 template <
class ELEMENT>
838 if (time_stepper_pt()->type()!=
"BDF")
840 std::ostringstream error_stream;
841 error_stream <<
"Timestepper has to be from the BDF family!\n"
842 <<
"You have specified a timestepper from the "
843 << time_stepper_pt()->type() <<
" family" << std::endl;
845 throw OomphLibError(error_stream.str(),
846 OOMPH_CURRENT_FUNCTION,
847 OOMPH_EXCEPTION_LOCATION);
851 bulk_mesh_pt()->node_update();
854 unsigned num_nod = bulk_mesh_pt()->nnode();
855 for (
unsigned n=0;n<num_nod;n++)
859 x[0]=bulk_mesh_pt()->node_pt(n)->x(0);
860 x[1]=bulk_mesh_pt()->node_pt(n)->x(1);
863 bulk_mesh_pt()->node_pt(n)->set_value(0,6.0*(x[1]/Ly)*(1.0-(x[1]/Ly)));
864 bulk_mesh_pt()->node_pt(n)->set_value(1,0.0);
868 bulk_mesh_pt()->assign_initial_values_impulsive();
869 wall_mesh_pt()->assign_initial_values_impulsive();
880 template<
class ELEMENT>
884 delete_traction_elements(Applied_fluid_traction_mesh_pt);
887 rebuild_global_mesh();
896 template<
class ELEMENT>
901 create_traction_elements(5,Bulk_mesh_pt,Applied_fluid_traction_mesh_pt);
904 rebuild_global_mesh();
907 RefineableNavierStokesEquations<2>::
908 unpin_all_pressure_dofs(Bulk_mesh_pt->element_pt());
911 RefineableNavierStokesEquations<2>::
912 pin_redundant_nodal_pressures(Bulk_mesh_pt->element_pt());
916 unsigned n_element=Applied_fluid_traction_mesh_pt->nelement();
917 for(
unsigned e=0;e<n_element;e++)
920 NavierStokesTractionElement<ELEMENT> *el_pt =
921 dynamic_cast<NavierStokesTractionElement<ELEMENT>*
>(
922 Applied_fluid_traction_mesh_pt->element_pt(e));
942 unsigned num_nod= bulk_mesh_pt()->nboundary_node(ibound);
943 for (
unsigned inod=0;inod<num_nod;inod++)
945 bulk_mesh_pt()->boundary_node_pt(ibound, inod)->
946 set_auxiliary_node_update_fct_pt(
947 FSI_functions::apply_no_slip_on_moving_wall);
958 FSI_functions::setup_fluid_load_info_for_solid_elements<ELEMENT,2>
959 (
this,3,Bulk_mesh_pt,Wall_mesh_pt);
971 int main(
int argc,
char* argv[])
975 CommandLineArgs::setup(argc,argv);
978 unsigned coarsening_factor=4;
979 if (CommandLineArgs::Argc>1)
985 unsigned nup=20/coarsening_factor;
986 unsigned ncollapsible=40/coarsening_factor;
987 unsigned ndown=40/coarsening_factor;
988 unsigned ny=16/coarsening_factor;
992 double lcollapsible=10.0;
1004 #ifdef MACRO_ELEMENT_NODE_UPDATE
1010 <MacroElementNodeUpdateElement<RefineableQTaylorHoodElement<2> > >
1011 problem(nup, ncollapsible, ndown, ny,
1012 lup, lcollapsible, ldown, ly);
1018 <MacroElementNodeUpdateElement<RefineableQCrouzeixRaviartElement<2> > >
1019 problem(nup, ncollapsible, ndown, ny,
1020 lup, lcollapsible, ldown, ly);
1030 <AlgebraicElement<RefineableQTaylorHoodElement<2> > >
1031 problem(nup, ncollapsible, ndown, ny,
1032 lup, lcollapsible, ldown, ly);
1038 <AlgebraicElement<RefineableQCrouzeixRaviartElement<2> > >
1039 problem(nup, ncollapsible, ndown, ny,
1040 lup, lcollapsible, ldown, ly);
1058 problem.time_pt()->time()=t_min;
1059 problem.initialise_dt(dt);
1062 problem.set_initial_condition();
1066 doc_info.set_directory(
"RESLT");
1069 ofstream trace_file;
1071 sprintf(filename,
"%s/trace.dat",doc_info.directory().c_str());
1072 trace_file.open(filename);
1075 problem.doc_solution(doc_info, trace_file);
1078 doc_info.number()++;
1081 unsigned nstep = unsigned((t_max-t_min)/dt);
1082 if (CommandLineArgs::Argc>1)
1089 problem.bulk_mesh_pt()->max_permitted_error()=1.0e-3;
1090 problem.bulk_mesh_pt()->min_permitted_error()=1.0e-5;
1093 if (CommandLineArgs::Argc>1)
1096 problem.bulk_mesh_pt()->max_permitted_error()=0.5e-2;
1097 problem.bulk_mesh_pt()->min_permitted_error()=0.5e-4;
1102 unsigned max_adapt=3;
1106 for (
unsigned istep=0;istep<nstep;istep++)
1109 problem.unsteady_newton_solve(dt,max_adapt,first);
1112 problem.doc_solution(doc_info, trace_file);
1115 doc_info.number()++;
MacroElementNodeUpdateRefineableCollapsibleChannelMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the "bulk" mesh.
void actions_after_newton_solve()
Update the problem after solve (empty)
MacroElementNodeUpdateRefineableCollapsibleChannelMesh< ELEMENT > * bulk_mesh_pt()
Access function for the specific bulk (fluid) mesh.
void actions_before_adapt()
Actions before adapt: Wipe the mesh of prescribed traction elements.
RefineableAlgebraicCollapsibleChannelMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the "bulk" mesh.
~FSICollapsibleChannelProblem()
Destructor (empty)
void actions_before_newton_convergence_check()
Update before checking Newton convergence: Update the nodal positions in the fluid mesh in response t...
void actions_before_newton_solve()
Update the problem specs before solve (empty)
void actions_after_adapt()
Actions after adapt: Rebuild the mesh of prescribed traction elements and reset FSI.
OneDLagrangianMesh< FSIHermiteBeamElement > * wall_mesh_pt()
Access function for the wall mesh.
void delete_traction_elements(Mesh *const &traction_mesh_pt)
Delete prescribed traction elements from the surface mesh.
void create_traction_elements(const unsigned &b, Mesh *const &bulk_mesh_pt, Mesh *const &traction_mesh_pt)
Create the prescribed traction elements on boundary b.
RefineableAlgebraicCollapsibleChannelMesh< ELEMENT > * bulk_mesh_pt()
Access function for the specific bulk (fluid) mesh.
void set_initial_condition()
Apply initial conditions.
FSICollapsibleChannelProblem(const unsigned &nup, const unsigned &ncollapsible, const unsigned &ndown, const unsigned &ny, const double &lup, const double &lcollapsible, const double &ldown, const double &ly)
Constructor: The arguments are the number of elements and the lengths of the domain.
void doc_solution(DocInfo &doc_info, ofstream &trace_file)
Doc the solution.
int main(int argc, char *argv[])
Driver code for a collapsible channel problem with FSI. Presence of command line arguments indicates ...
Namespace to define the mapping [0,1] -> [0,1] that re-distributes nodal points across the channel wi...
double squash_fct(const double &s)
Mapping [0,1] -> [0,1] that re-distributes nodal points across the channel width.
double Delta
Boundary layer width.
double Fract_in_BL
Fraction of points in boundary layer.
Namespace for phyical parameters.
double P_ext
External pressure.
double ReSt
Womersley = Reynolds times Strouhal.
void prescribed_traction(const double &t, const Vector< double > &x, const Vector< double > &n, Vector< double > &traction)
Traction applied on the fluid at the left (inflow) boundary.
void load(const Vector< double > &xi, const Vector< double > &x, const Vector< double > &N, Vector< double > &load)
Load function: Apply a constant external pressure to the wall. Note: This is the load without the flu...
double Sigma0
2nd Piola Kirchhoff pre-stress. As in Jensen & Heil (2003) paper.
double Q
Fluid structure interaction parameter: Ratio of stresses used for non-dimensionalisation of fluid to ...
double Re
Reynolds number.
double P_up
Default pressure on the left boundary.
double H
Non-dimensional wall thickness. As in Jensen & Heil (2003) paper.