29 #include "navier_stokes.h"
32 #include "constitutive.h"
35 #include "meshes/one_d_lagrangian_mesh.h"
38 #include "meshes/collapsible_channel_mesh.h"
42 using namespace oomph;
48 template <
class ELEMENT>
50 public virtual RefineableCollapsibleChannelMesh<ELEMENT>,
51 public virtual SolidMesh
60 const unsigned& ncollapsible,
61 const unsigned& ndown,
64 const double& lcollapsible,
68 TimeStepper* time_stepper_pt=
69 &Mesh::Default_TimeStepper) :
70 CollapsibleChannelMesh<ELEMENT>(nup, ncollapsible, ndown, ny,
71 lup, lcollapsible, ldown, ly,
74 RefineableCollapsibleChannelMesh<ELEMENT>(nup, ncollapsible, ndown, ny,
75 lup, lcollapsible, ldown, ly,
83 set_lagrangian_nodal_coordinates();
155 void position(
const Vector<double>& zeta, Vector<double>& r)
const
166 void position(
const unsigned& t,
const Vector<double>& zeta,
167 Vector<double>& r)
const
182 DenseMatrix<double> &drdzeta,
183 RankThreeTensor<double> &ddrdzeta)
const
237 const Vector<double>& x,
238 const Vector<double>& n,
239 Vector<double>& traction)
259 void load(
const Vector<double>& xi,
const Vector<double>& x,
260 const Vector<double>& N, Vector<double>&
load)
262 for(
unsigned i=0;i<2;i++)
282 template <
class ELEMENT>
291 const unsigned& ncollapsible,
292 const unsigned& ndown,
295 const double& lcollapsible,
340 unsigned num_nod=bulk_mesh_pt()->nboundary_node(ibound);
341 for (
unsigned inod=0;inod<num_nod;inod++)
344 Node* node_pt=bulk_mesh_pt()->boundary_node_pt(ibound,inod);
347 FSI_functions::apply_no_slip_on_moving_wall(node_pt);
362 Mesh*
const &bulk_mesh_pt,
363 Mesh*
const &traction_mesh_pt);
381 unsigned Ncollapsible;
406 Mesh* Applied_fluid_traction_mesh_pt;
409 SolidMesh* Lagrange_multiplier_mesh_pt;
412 OneDLagrangianMesh<FSIHermiteBeamElement>* Wall_mesh_pt;
415 MeshAsGeomObject* Wall_geom_object_pt;
427 ConstitutiveLaw *Constitutive_law_pt;
437 template <
class ELEMENT>
440 const unsigned& ncollapsible,
441 const unsigned& ndown,
444 const double& lcollapsible,
450 Ncollapsible=ncollapsible;
454 Lcollapsible=lcollapsible;
460 Problem::Max_residuals=1000.0;
463 BDF<2>* fluid_time_stepper_pt=
new BDF<2>;
466 add_time_stepper_pt(fluid_time_stepper_pt);
469 Steady<2>* wall_time_stepper_pt =
new Steady<2>;
472 add_time_stepper_pt(wall_time_stepper_pt);
480 Wall_mesh_pt =
new OneDLagrangianMesh<FSIHermiteBeamElement>
481 (Ncollapsible,Lcollapsible,undeformed_wall_pt,wall_time_stepper_pt);
486 new MeshAsGeomObject(Wall_mesh_pt);
493 (nup, ncollapsible, ndown, ny,
494 lup, lcollapsible, ldown, ly,
496 fluid_time_stepper_pt);
502 bool update_all_solid_nodes=
true;
503 Bulk_mesh_pt->node_update(update_all_solid_nodes);
504 Bulk_mesh_pt->set_lagrangian_nodal_coordinates();
506 if (CommandLineArgs::Argc>1)
509 Bulk_mesh_pt->refine_uniformly();
515 Applied_fluid_traction_mesh_pt =
new Mesh;
519 create_traction_elements(5,Bulk_mesh_pt,Applied_fluid_traction_mesh_pt);
523 Lagrange_multiplier_mesh_pt=
new SolidMesh;
524 create_lagrange_multiplier_elements();
527 add_sub_mesh(Bulk_mesh_pt);
528 add_sub_mesh(Applied_fluid_traction_mesh_pt);
529 add_sub_mesh(Wall_mesh_pt);
530 add_sub_mesh(Lagrange_multiplier_mesh_pt);
536 Z2ErrorEstimator* error_estimator_pt=
new Z2ErrorEstimator;
537 bulk_mesh_pt()->spatial_error_estimator_pt()=error_estimator_pt;
548 unsigned n_element=Bulk_mesh_pt->nelement();
549 for(
unsigned e=0;e<n_element;e++)
553 ELEMENT* el_pt =
dynamic_cast<ELEMENT*
>(Bulk_mesh_pt->element_pt(e));
562 el_pt->constitutive_law_pt() = Constitutive_law_pt;
571 RefineableNavierStokesEquations<2>::
572 pin_redundant_nodal_pressures(Bulk_mesh_pt->element_pt());
582 unsigned num_nod= bulk_mesh_pt()->nboundary_node(ibound);
583 for (
unsigned inod=0;inod<num_nod;inod++)
585 for(
unsigned i=0;i<2;i++)
587 bulk_mesh_pt()->boundary_node_pt(ibound, inod)->pin(i);
592 for(ibound=2;ibound<5;ibound++)
594 num_nod= bulk_mesh_pt()->nboundary_node(ibound);
595 for (
unsigned inod=0;inod<num_nod;inod++)
597 for(
unsigned i=0;i<2;i++)
599 bulk_mesh_pt()->boundary_node_pt(ibound, inod)->pin(i);
606 num_nod= bulk_mesh_pt()->nboundary_node(ibound);
607 for (
unsigned inod=0;inod<num_nod;inod++)
609 bulk_mesh_pt()->boundary_node_pt(ibound, inod)->pin(1);
615 num_nod= bulk_mesh_pt()->nboundary_node(ibound);
616 for (
unsigned inod=0;inod<num_nod;inod++)
619 bulk_mesh_pt()->boundary_node_pt(ibound, inod)->pin(1);
625 for (
unsigned ibound=0;ibound<6;ibound++)
629 unsigned num_nod= bulk_mesh_pt()->nboundary_node(ibound);
630 for (
unsigned inod=0;inod<num_nod;inod++)
632 for(
unsigned i=0;i<2;i++)
634 dynamic_cast<SolidNode*
>(bulk_mesh_pt()->
635 boundary_node_pt(ibound, inod))
643 unsigned nnod=bulk_mesh_pt()->nnode();
644 for (
unsigned j=0;j<nnod;j++)
646 SolidNode* nod_pt=
dynamic_cast<SolidNode*
>(bulk_mesh_pt()->node_pt(j));
647 if ((nod_pt->x(0)<=lup)||(nod_pt->x(0)>=(lup+lcollapsible)))
649 for(
unsigned i=0;i<2;i++)
651 nod_pt->pin_position(i);
662 unsigned n_el=Applied_fluid_traction_mesh_pt->nelement();
663 for(
unsigned e=0;e<n_el;e++)
666 NavierStokesTractionElement<ELEMENT> *el_pt =
667 dynamic_cast< NavierStokesTractionElement<ELEMENT>*
>(
668 Applied_fluid_traction_mesh_pt->element_pt(e));
682 n_element = wall_mesh_pt()->nelement();
683 for(
unsigned e=0;e<n_element;e++)
686 FSIHermiteBeamElement *elem_pt =
687 dynamic_cast<FSIHermiteBeamElement*
>(wall_mesh_pt()->element_pt(e));
700 elem_pt->undeformed_beam_pt() = undeformed_wall_pt;
705 elem_pt->set_normal_pointing_out_of_fluid();
716 for(
unsigned b=0;b<2;b++)
719 wall_mesh_pt()->boundary_node_pt(b,0)->pin_position(0);
720 wall_mesh_pt()->boundary_node_pt(b,0)->pin_position(1);
730 num_nod= bulk_mesh_pt()->nboundary_node(ibound);
731 unsigned control_nod=num_nod/2;
732 Left_node_pt= bulk_mesh_pt()->boundary_node_pt(ibound, control_nod);
736 num_nod= bulk_mesh_pt()->nboundary_node(ibound);
737 control_nod=num_nod/2;
738 Right_node_pt= bulk_mesh_pt()->boundary_node_pt(ibound, control_nod);
742 num_nod= wall_mesh_pt()->nnode();
743 Wall_node_pt=wall_mesh_pt()->node_pt(Ncollapsible/2);
756 num_nod= bulk_mesh_pt()->nboundary_node(ibound);
757 for (
unsigned inod=0;inod<num_nod;inod++)
759 bulk_mesh_pt()->boundary_node_pt(ibound, inod)->
760 set_auxiliary_node_update_fct_pt(
761 FSI_functions::apply_no_slip_on_moving_wall);
768 FSI_functions::setup_fluid_load_info_for_solid_elements<ELEMENT,2>
769 (
this,3,Bulk_mesh_pt,Wall_mesh_pt);
772 cout <<
"Number of equations: " << assign_eqn_numbers() << std::endl;
783 template <
class ELEMENT>
785 ofstream& trace_file)
795 sprintf(filename,
"%s/soln%i.dat",doc_info.directory().c_str(),
797 some_file.open(filename);
798 bulk_mesh_pt()->output(some_file,npts);
802 sprintf(filename,
"%s/beam%i.dat",doc_info.directory().c_str(),
804 some_file.open(filename);
805 wall_mesh_pt()->output(some_file,npts);
811 unsigned nsteps=time_stepper_pt(1)->nprev_values();
812 for (
unsigned t=0;t<=nsteps;t++)
814 sprintf(filename,
"%s/wall%i-%i.dat",doc_info.directory().c_str(),
815 doc_info.number(),t);
816 some_file.open(filename);
817 unsigned n_elem=wall_mesh_pt()->nelement();
818 for (
unsigned ielem=0;ielem<n_elem;ielem++)
820 dynamic_cast<FSIHermiteBeamElement*
>(wall_mesh_pt()->element_pt(ielem))->
821 output(t,some_file,npts);
829 trace_file << time_pt()->time() <<
" "
830 << Wall_node_pt->x(1) <<
" "
831 << Left_node_pt->value(0) <<
" "
832 << Right_node_pt->value(0) <<
" "
845 template <
class ELEMENT>
847 const unsigned &b, Mesh*
const &bulk_mesh_pt, Mesh*
const &traction_mesh_pt)
851 unsigned n_element = bulk_mesh_pt->nboundary_element(b);
854 for(
unsigned e=0;e<n_element;e++)
857 ELEMENT* bulk_elem_pt =
dynamic_cast<ELEMENT*
>
858 (bulk_mesh_pt->boundary_element_pt(b,e));
861 int face_index = bulk_mesh_pt->face_index_at_boundary(b,e);
864 NavierStokesTractionElement<ELEMENT>* flux_element_pt =
865 new NavierStokesTractionElement<ELEMENT>(bulk_elem_pt,face_index);
868 traction_mesh_pt->add_element_pt(flux_element_pt);
879 template<
class ELEMENT>
884 unsigned n_element = surface_mesh_pt->nelement();
887 for(
unsigned e=0;e<n_element;e++)
890 delete surface_mesh_pt->element_pt(e);
894 surface_mesh_pt->flush_element_and_node_storage();
902 template<
class ELEMENT>
910 unsigned n_element = bulk_mesh_pt()->nboundary_element(b);
913 for(
unsigned e=0;e<n_element;e++)
916 ELEMENT* bulk_elem_pt =
dynamic_cast<ELEMENT*
>(
917 bulk_mesh_pt()->boundary_element_pt(b,e));
920 int face_index = bulk_mesh_pt()->face_index_at_boundary(b,e);
923 Lagrange_multiplier_mesh_pt->add_element_pt(
924 new ImposeDisplacementByLagrangeMultiplierElement<ELEMENT>(
925 bulk_elem_pt,face_index));
931 n_element=Lagrange_multiplier_mesh_pt->nelement();
932 for(
unsigned i=0;i<n_element;i++)
935 ImposeDisplacementByLagrangeMultiplierElement<ELEMENT> *el_pt =
936 dynamic_cast<ImposeDisplacementByLagrangeMultiplierElement<ELEMENT>*
>
937 (Lagrange_multiplier_mesh_pt->element_pt(i));
942 el_pt->set_boundary_shape_geom_object_pt(Wall_geom_object_pt,b);
945 unsigned nnod=el_pt->nnode();
946 for (
unsigned j=0;j<nnod;j++)
948 Node* nod_pt = el_pt->node_pt(j);
951 if ((nod_pt->is_on_boundary(2))||(nod_pt->is_on_boundary(4)))
955 unsigned n_bulk_value=el_pt->nbulk_value(j);
958 unsigned nval=nod_pt->nvalue();
959 for (
unsigned j=n_bulk_value;j<nval;j++)
976 template<
class ELEMENT>
981 unsigned n_element = Lagrange_multiplier_mesh_pt->nelement();
984 for(
unsigned e=0;e<n_element;e++)
987 delete Lagrange_multiplier_mesh_pt->element_pt(e);
991 Lagrange_multiplier_mesh_pt->flush_element_and_node_storage();
1000 template <
class ELEMENT>
1004 if (time_stepper_pt()->type()!=
"BDF")
1006 std::ostringstream error_stream;
1007 error_stream <<
"Timestepper has to be from the BDF family!\n"
1008 <<
"You have specified a timestepper from the "
1009 << time_stepper_pt()->type() <<
" family" << std::endl;
1011 throw OomphLibError(error_stream.str(),
1012 OOMPH_CURRENT_FUNCTION,
1013 OOMPH_EXCEPTION_LOCATION);
1017 unsigned num_nod = bulk_mesh_pt()->nnode();
1018 for (
unsigned n=0;n<num_nod;n++)
1021 Vector<double> x(2);
1022 x[0]=bulk_mesh_pt()->node_pt(n)->x(0);
1023 x[1]=bulk_mesh_pt()->node_pt(n)->x(1);
1026 bulk_mesh_pt()->node_pt(n)->set_value(0,6.0*(x[1]/Ly)*(1.0-(x[1]/Ly)));
1027 bulk_mesh_pt()->node_pt(n)->set_value(1,0.0);
1031 bulk_mesh_pt()->assign_initial_values_impulsive();
1032 wall_mesh_pt()->assign_initial_values_impulsive();
1043 template<
class ELEMENT>
1047 delete_traction_elements(Applied_fluid_traction_mesh_pt);
1050 delete_lagrange_multiplier_elements();
1053 rebuild_global_mesh();
1062 template<
class ELEMENT>
1067 create_traction_elements(5,Bulk_mesh_pt,Applied_fluid_traction_mesh_pt);
1072 create_lagrange_multiplier_elements();
1075 rebuild_global_mesh();
1078 RefineableNavierStokesEquations<2>::
1079 unpin_all_pressure_dofs(Bulk_mesh_pt->element_pt());
1082 RefineableNavierStokesEquations<2>::
1083 pin_redundant_nodal_pressures(Bulk_mesh_pt->element_pt());
1087 unsigned n_element=Applied_fluid_traction_mesh_pt->nelement();
1088 for(
unsigned e=0;e<n_element;e++)
1091 NavierStokesTractionElement<ELEMENT> *el_pt =
1092 dynamic_cast<NavierStokesTractionElement<ELEMENT>*
>(
1093 Applied_fluid_traction_mesh_pt->element_pt(e));
1113 unsigned num_nod=bulk_mesh_pt()->nboundary_node(ibound);
1114 for (
unsigned inod=0;inod<num_nod;inod++)
1116 bulk_mesh_pt()->boundary_node_pt(ibound, inod)->
1117 set_auxiliary_node_update_fct_pt(
1118 FSI_functions::apply_no_slip_on_moving_wall);
1129 FSI_functions::setup_fluid_load_info_for_solid_elements<ELEMENT,2>
1130 (
this,3,Bulk_mesh_pt,Wall_mesh_pt);
1144 CommandLineArgs::setup(argc,argv);
1147 unsigned coarsening_factor=4;
1148 if (CommandLineArgs::Argc>1)
1150 coarsening_factor=4;
1154 unsigned nup=20/coarsening_factor;
1155 unsigned ncollapsible=40/coarsening_factor;
1156 unsigned ndown=40/coarsening_factor;
1157 unsigned ny=16/coarsening_factor;
1161 double lcollapsible=10.0;
1177 <RefineablePseudoSolidNodeUpdateElement<
1178 RefineableQTaylorHoodElement<2>,
1179 RefineableQPVDElement<2,3> > >
1180 problem(nup, ncollapsible, ndown, ny,
1181 lup, lcollapsible, ldown, ly);
1187 <RefineablePseudoSolidNodeUpdateElement<
1188 RefineableQCrouzeixRaviartElement<2>,
1189 RefineableQPVDElement<2,3> > >
1190 problem(nup, ncollapsible, ndown, ny,
1191 lup, lcollapsible, ldown, ly);
1206 problem.time_pt()->time()=t_min;
1207 problem.initialise_dt(dt);
1210 problem.set_initial_condition();
1214 doc_info.set_directory(
"RESLT");
1217 ofstream trace_file;
1219 sprintf(filename,
"%s/trace.dat",doc_info.directory().c_str());
1220 trace_file.open(filename);
1223 problem.doc_solution(doc_info, trace_file);
1226 doc_info.number()++;
1229 unsigned nstep = unsigned((t_max-t_min)/dt);
1230 if (CommandLineArgs::Argc>1)
1237 problem.bulk_mesh_pt()->max_permitted_error()=1.0e-3;
1238 problem.bulk_mesh_pt()->min_permitted_error()=1.0e-5;
1241 if (CommandLineArgs::Argc>1)
1244 problem.bulk_mesh_pt()->max_permitted_error()=0.5e-2;
1245 problem.bulk_mesh_pt()->min_permitted_error()=0.5e-4;
1250 unsigned max_adapt=3;
1254 for (
unsigned istep=0;istep<nstep;istep++)
1257 problem.unsteady_newton_solve(dt,max_adapt,first);
1260 problem.doc_solution(doc_info, trace_file);
1263 doc_info.number()++;
Upgrade mesh to solid mesh.
void create_lagrange_multiplier_elements()
Create elements that enforce prescribed boundary motion by Lagrange multiplilers.
void actions_after_newton_solve()
Update the problem after solve (empty)
void actions_before_adapt()
Actions before adapt: Wipe the mesh of prescribed traction elements.
void delete_lagrange_multiplier_elements()
Delete elements that enforce prescribed boundary motion by Lagrange multiplilers.
~FSICollapsibleChannelProblem()
Destructor (empty)
ElasticRefineableCollapsibleChannelMesh< ELEMENT > * bulk_mesh_pt()
Access function for the specific bulk (fluid) mesh.
ElasticRefineableCollapsibleChannelMesh< ELEMENT > * Bulk_mesh_pt
Pointer to the "bulk" mesh.
void actions_before_newton_convergence_check()
Update no slip before Newton convergence check.
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.
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 Nu
Pseudo-solid Poisson ratio.
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 Lambda_sq
Pseudo-solid mass density.
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.