40 #include "navier_stokes.h"
43 #include "meshes/cylinder_with_flag_mesh.h"
47 using namespace oomph;
99 if (t<=0.0) tmp_ampl=0.0;
100 Vector<double> uppertip(2);
103 tmp_ampl*sin(2.0*MathematicalConstants::Pi*t/
Period);
111 if (t<=0.0) tmp_ampl=0.0;
112 Vector<double> lowertip(2);
115 tmp_ampl*sin(2.0*MathematicalConstants::Pi*t/
Period);
139 void position(
const unsigned& t,
const Vector<double> &xi, Vector<double> &r)
143 Vector<double> point_on_circle(2);
152 1.0/3.0*sin((r[0]-point_on_circle[0])/
154 point_on_circle[0])*MathematicalConstants::Pi)
155 *sin(2.0* MathematicalConstants::Pi*
Time_pt->time(t)/
Period);
161 void position(
const Vector<double> &xi, Vector<double> &r)
const
163 return position(0,xi,r);
190 void position(
const unsigned& t,
const Vector<double> &xi, Vector<double> &r)
194 Vector<double> point_on_circle(2);
203 1.0/3.0*sin((r[0]-point_on_circle[0])/
205 point_on_circle[0])*MathematicalConstants::Pi)
211 void position(
const Vector<double> &xi, Vector<double> &r)
const
213 return position(0,xi,r);
241 void position(
const unsigned& t,
const Vector<double> &xi, Vector<double> &r)
247 r[1]= point_bottom[1]+(xi[0]+
H/2.0)/
H*(point_top[1]-point_bottom[1]);
248 r[0]= point_bottom[0]+(xi[0]+
H/2.0)/
H*(point_top[0]-point_bottom[0]);
252 void position(
const Vector<double> &xi, Vector<double> &r)
const
254 return position(0,xi,r);
310 template<
class ELEMENT>
318 const double &height);
327 void actions_after_adapt();
330 void actions_before_implicit_timestep();
333 #ifdef USE_MACRO_ELEMENTS
336 RefineableCylinderWithFlagMesh<ELEMENT>*
mesh_pt()
338 return dynamic_cast<RefineableCylinderWithFlagMesh<ELEMENT>*
>
339 (Problem::mesh_pt());
345 RefineableAlgebraicCylinderWithFlagMesh<ELEMENT>*
mesh_pt()
347 return dynamic_cast<RefineableAlgebraicCylinderWithFlagMesh<ELEMENT>*
>
348 (Problem::mesh_pt());
354 void doc_solution(DocInfo& doc_info);
370 template<
class ELEMENT>
372 const double &length,
const double &height) : Height(height)
376 Max_residuals=100.0;;
377 Max_newton_iterations=50;
380 add_time_stepper_pt(
new BDF<2>);
385 #ifdef USE_MACRO_ELEMENTS
390 Problem::mesh_pt()=
new RefineableCylinderWithFlagMesh<ELEMENT>
403 std::cout <<
"Using Domain/MacroElement-based node update" << std::endl;
408 Problem::mesh_pt()=
new RefineableAlgebraicCylinderWithFlagMesh<ELEMENT>
421 std::cout <<
"Using AlgebraicMesh-based node update" << std::endl;
430 Z2ErrorEstimator* error_estimator_pt=
new Z2ErrorEstimator;
431 mesh_pt()->spatial_error_estimator_pt()=error_estimator_pt;
438 unsigned num_bound =
mesh_pt()->nboundary();
439 for(
unsigned ibound=0;ibound<num_bound;ibound++)
441 unsigned num_nod=
mesh_pt()->nboundary_node(ibound);
442 for (
unsigned inod=0;inod<num_nod;inod++)
447 mesh_pt()->boundary_node_pt(ibound,inod)->pin(0);
448 mesh_pt()->boundary_node_pt(ibound,inod)->pin(1);
452 mesh_pt()->boundary_node_pt(ibound,inod)->pin(1);
459 unsigned nelem=
mesh_pt()->nelement();
460 for (
unsigned e=0;e<nelem;e++)
463 ELEMENT* el_pt =
dynamic_cast<ELEMENT*
>(
mesh_pt()->element_pt(e));
475 unsigned num_nod=
mesh_pt()->nboundary_node(ibound);
476 for (
unsigned inod=0;inod<num_nod;inod++)
478 double ycoord =
mesh_pt()->boundary_node_pt(ibound,inod)->x(1);
483 mesh_pt()->boundary_node_pt(ibound,inod)->set_value(0,uy);
484 mesh_pt()->boundary_node_pt(ibound,inod)->set_value(1,0.0);
489 RefineableNavierStokesEquations<2>::
490 pin_redundant_nodal_pressures(
mesh_pt()->element_pt());
493 cout <<
"Number of equations: " << assign_eqn_numbers() << endl;
502 template <
class ELEMENT>
506 RefineableNavierStokesEquations<2>::
507 unpin_all_pressure_dofs(mesh_pt()->element_pt());
510 RefineableNavierStokesEquations<2>::
511 pin_redundant_nodal_pressures(mesh_pt()->element_pt());
519 template <
class ELEMENT>
524 mesh_pt()->node_update();
528 for(
unsigned ibound=5;ibound<8;ibound++)
530 unsigned num_nod=mesh_pt()->nboundary_node(ibound);
531 for (
unsigned inod=0;inod<num_nod;inod++)
534 Node* node_pt=mesh_pt()->boundary_node_pt(ibound,inod);
537 FSI_functions::apply_no_slip_on_moving_wall(node_pt);
546 template<
class ELEMENT>
558 sprintf(filename,
"%s/soln%i.dat",doc_info.directory().c_str(),
560 some_file.open(filename);
561 mesh_pt()->output(some_file,npts);
576 int main(
int argc,
char* argv[])
584 doc_info.set_directory(
"RESLT");
591 #ifdef USE_MACRO_ELEMENTS
595 problem(length,height);
601 <AlgebraicElement<RefineableQCrouzeixRaviartElement<2> > >
602 problem(length, height);
608 unsigned nsteps_per_period=40;
614 unsigned nstep=nsteps_per_period*nperiod;
615 if (CommandLineArgs::Argc>1)
626 problem.initialise_dt(dt);
630 unsigned max_adapt=3;
631 if (CommandLineArgs::Argc>1)
643 problem.steady_newton_solve(max_adapt);
657 for (
unsigned istep=0;istep<nstep;istep++)
660 problem.unsteady_newton_solve(dt,max_adapt,first);
GeomObject that defines the lower boundary of the flag.
void position(const unsigned &t, const Vector< double > &xi, Vector< double > &r) const
Return the position along the bottom of the flag (xi[0] varies between 0 and Lx)
BottomOfFlag()
Constructor:
unsigned ngeom_data() const
Number of geometric Data in GeomObject: None.
void position(const Vector< double > &xi, Vector< double > &r) const
Current position.
~BottomOfFlag()
Destructor (empty)
GeomObject that defines the tip of the flag.
unsigned ngeom_data() const
Number of geometric Data in GeomObject: None.
~TipOfFlag()
Destructor (empty)
void position(const Vector< double > &xi, Vector< double > &r) const
Current position.
void position(const unsigned &t, const Vector< double > &xi, Vector< double > &r) const
Return the position This object links the tips of the top and bottom by a straight line whilst xi[0] ...
GeomObject that defines the upper boundary of the flag.
void position(const Vector< double > &xi, Vector< double > &r) const
Current position.
TopOfFlag()
Constructor: It's a 1D object in 2D space.
unsigned ngeom_data() const
Number of geometric Data in GeomObject: None.
void position(const unsigned &t, const Vector< double > &xi, Vector< double > &r) const
Return the position along the top of the flag (xi[0] varies between 0 and Lx)
~TopOfFlag()
Destructor (empty)
///////////////////////////////////////////////////////////////// ///////////////////////////////////...
void actions_after_newton_solve()
Update the problem specs after solve (empty)
void actions_before_implicit_timestep()
Update the velocity boundary condition on the flag.
void doc_solution(DocInfo &doc_info)
Doc the solution.
RefineableAlgebraicCylinderWithFlagMesh< ELEMENT > * mesh_pt()
Access function for the specific mesh.
RefineableCylinderWithFlagMesh< ELEMENT > * mesh_pt()
Access function for the specific mesh.
double Height
Height of channel.
void actions_after_adapt()
After adaptation: Unpin pressures and pin redudant pressure dofs.
TurekNonFSIProblem(const double &length, const double &height)
Constructor: Pass length and height of domain.
void actions_before_newton_solve()
Update the problem specs before solve (empty)
Namespace for definition of flag boundaries.
double Centre_y
y position of centre of cylinder
TipOfFlag * Tip_flag_pt
Pointer to GeomObject that bounds the tip edge of the flag.
double Period
Period of prescribed flag oscillation.
double Centre_x
x position of centre of cylinder
void setup(Time *time_pt)
Create all GeomObjects needed to define the cylinder and the flag.
Vector< double > upper_tip(const double &t)
Time-dependent vector to upper tip of the "flag".
Circle * Cylinder_pt
Pointer to GeomObject of type Circle that defines the central cylinder.
Vector< double > lower_tip(const double &t)
Time-dependent vector to bottom tip of the "flag".
double Radius
Radius of cylinder.
double Amplitude
Amplitude of tip deflection.
Time * Time_pt
Pointer to the global time object.
BottomOfFlag * Bottom_flag_pt
Pointer to GeomObject that bounds the bottom edge of the flag.
TopOfFlag * Top_flag_pt
Pointer to GeomObject that bounds the upper edge of the flag.
double Re
Reynolds number.
int main(int argc, char *argv[])
//////////////////////////////////////////////////////////////////// ////////////////////////////////...