unstructured_two_d_fluid.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// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7//LIC//
8//LIC// This library is free software; you can redistribute it and/or
9//LIC// modify it under the terms of the GNU Lesser General Public
10//LIC// License as published by the Free Software Foundation; either
11//LIC// version 2.1 of the License, or (at your option) any later version.
12//LIC//
13//LIC// This library is distributed in the hope that it will be useful,
14//LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15//LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16//LIC// Lesser General Public License for more details.
17//LIC//
18//LIC// You should have received a copy of the GNU Lesser General Public
19//LIC// License along with this library; if not, write to the Free Software
20//LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21//LIC// 02110-1301 USA.
22//LIC//
23//LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24//LIC//
25//LIC//====================================================================
26// Driver for flow past odd-shaped obstacle -- domain meshed with triangle.
27// This is a warm-up problem for an unstructured fsi problem.
28
29//Generic includes
30#include "generic.h"
31#include "navier_stokes.h"
32#include "solid.h"
33#include "constitutive.h"
34
35// The mesh
36#include "meshes/triangle_mesh.h"
37
38
39using namespace std;
40
41using namespace oomph;
42
43
44//===========start_mesh====================================================
45/// Triangle-based mesh upgraded to become a (pseudo-) solid mesh
46//=========================================================================
47template<class ELEMENT>
48class ElasticTriangleMesh : public virtual TriangleMesh<ELEMENT>,
49 public virtual SolidMesh
50{
51
52public:
53
54 /// Constructor:
55 ElasticTriangleMesh(const std::string& node_file_name,
56 const std::string& element_file_name,
57 const std::string& poly_file_name,
58 TimeStepper* time_stepper_pt=
59 &Mesh::Default_TimeStepper) :
60 TriangleMesh<ELEMENT>(node_file_name,element_file_name,
61 poly_file_name, time_stepper_pt)
62 {
63 //Assign the Lagrangian coordinates
64 set_lagrangian_nodal_coordinates();
65
66 //Identify the domain boundaries
67 this->identify_boundaries();
68 }
69
70
71 /// Function used to identify the domain boundaries
73 {
74 // We will have three boundaries
75 this->set_nboundary(3);
76
77 unsigned n_node=this->nnode();
78 for (unsigned j=0;j<n_node;j++)
79 {
80 Node* nod_pt=this->node_pt(j);
81
82 // Boundary 1 is left (inflow) boundary
83 if (nod_pt->x(0)<0.226)
84 {
85 // If it's not on the upper or lower channel walls remove it
86 // from boundary 0
87 if ((nod_pt->x(1)<4.08)&&(nod_pt->x(1)>0.113))
88 {
89 this->remove_boundary_node(0,nod_pt);
90 }
91 this->add_boundary_node(1,nod_pt);
92 }
93
94 // Boundary 2 is right (outflow) boundary
95 if (nod_pt->x(0)>8.28)
96 {
97 // If it's not on the upper or lower channel walls remove it
98 // from boundary 0
99 if ((nod_pt->x(1)<4.08)&&(nod_pt->x(1)>0.113))
100 {
101 this->remove_boundary_node(0,nod_pt);
102 }
103 this->add_boundary_node(2,nod_pt);
104 }
105 }
106 this->setup_boundary_element_info();
107 }
108
109 /// Empty Destructor
111
112
113};
114
115
116
117//==start_namespace==============================
118/// Namespace for physical parameters
119//==================================================
121{
122
123 /// Reynolds number
124 double Re=0.0;
125
126 /// Pseudo-solid Poisson ratio
127 double Nu=0.3;
128
129 /// Constitutive law used to determine the mesh deformation
130 ConstitutiveLaw *Constitutive_law_pt=0;
131
132} // end_of_namespace
133
134
135
136/// /////////////////////////////////////////////////////////////////////
137/// /////////////////////////////////////////////////////////////////////
138/// /////////////////////////////////////////////////////////////////////
139
140
141
142//==start_of_problem_class============================================
143/// Unstructured fluid problem
144//====================================================================
145template<class ELEMENT>
146class UnstructuredFluidProblem : public Problem
147{
148
149public:
150
151 /// Constructor
153
154 /// Destructor (empty)
156
157 /// Access function for the fluid mesh
159 {
160 return Fluid_mesh_pt;
161 }
162
163 /// Set the boundary conditions
165
166 /// Doc the solution
167 void doc_solution(DocInfo& doc_info);
168
169private:
170
171 /// Fluid mesh
173
174
175}; // end_of_problem_class
176
177
178//==start_constructor=====================================================
179/// Constructor
180//========================================================================
181template<class ELEMENT>
183{
184
185 //Create fluid mesh
186 string fluid_node_file_name="fluid.fig.1.node";
187 string fluid_element_file_name="fluid.fig.1.ele";
188 string fluid_poly_file_name="fluid.fig.1.poly";
189 Fluid_mesh_pt = new ElasticTriangleMesh<ELEMENT>(fluid_node_file_name,
190 fluid_element_file_name,
191 fluid_poly_file_name);
192
193 //Set the boundary conditions
194 this->set_boundary_conditions();
195
196 // Add fluid mesh
197 add_sub_mesh(fluid_mesh_pt());
198
199 // Build global mesh
200 build_global_mesh();
201
202 // Setup equation numbering scheme
203 cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;
204
205} // end_of_constructor
206
207
208
209//============start_of_set_boundary_conditions=============================
210/// Set the boundary conditions for the problem and document the boundary
211/// nodes
212//==========================================================================
213template<class ELEMENT>
215{
216 // Doc pinned nodes
217 std::ofstream solid_bc_file("pinned_solid_nodes.dat");
218 std::ofstream u_bc_file("pinned_u_nodes.dat");
219 std::ofstream v_bc_file("pinned_v_nodes.dat");
220
221 // Set the boundary conditions for fluid problem: All nodes are
222 // free by default -- just pin the ones that have Dirichlet conditions
223 // here.
224 unsigned nbound=Fluid_mesh_pt->nboundary();
225 for(unsigned ibound=0;ibound<nbound;ibound++)
226 {
227 unsigned num_nod=Fluid_mesh_pt->nboundary_node(ibound);
228 for (unsigned inod=0;inod<num_nod;inod++)
229 {
230 // Get node
231 Node* nod_pt=Fluid_mesh_pt->boundary_node_pt(ibound,inod);
232
233 // Pin velocity everywhere apart from outlet where we
234 // have parallel outflow
235 if (ibound!=2)
236 {
237 // Pin it
238 nod_pt->pin(0);
239 // Doc it as pinned
240 u_bc_file << nod_pt->x(0) << " " << nod_pt->x(1) << std::endl;
241 }
242 // Pin it
243 nod_pt->pin(1);
244 // Doc it as pinned
245 v_bc_file << nod_pt->x(0) << " " << nod_pt->x(1) << std::endl;
246
247 // Pin pseudo-solid positions everywhere
248 for(unsigned i=0;i<2;i++)
249 {
250 // Pin the node
251 SolidNode* nod_pt=Fluid_mesh_pt->boundary_node_pt(ibound,inod);
252 nod_pt->pin_position(i);
253
254 // Doc it as pinned
255 solid_bc_file << nod_pt->x(i) << " ";
256 }
257 solid_bc_file << std::endl;
258 }
259 } // end loop over boundaries
260
261 // Close
262 solid_bc_file.close();
263
264 // Complete the build of all elements so they are fully functional
265 unsigned n_element = fluid_mesh_pt()->nelement();
266 for(unsigned e=0;e<n_element;e++)
267 {
268 // Upcast from GeneralisedElement to the present element
269 ELEMENT* el_pt = dynamic_cast<ELEMENT*>(fluid_mesh_pt()->element_pt(e));
270
271 //Set the Reynolds number
272 el_pt->re_pt() = &Global_Physical_Variables::Re;
273
274 // Set the constitutive law
275 el_pt->constitutive_law_pt() =
277 }
278
279
280 // Apply fluid boundary conditions: Poiseuille at inflow
281 //------------------------------------------------------
282
283 // Find max. and min y-coordinate at inflow
284 unsigned ibound=1;
285 //Initialise y_min and y_max to y-coordinate of first node on boundary
286 double y_min=fluid_mesh_pt()->boundary_node_pt(ibound,0)->x(1);
287 double y_max=y_min;
288 //Loop over the rest of the nodes
289 unsigned num_nod= fluid_mesh_pt()->nboundary_node(ibound);
290 for (unsigned inod=1;inod<num_nod;inod++)
291 {
292 double y=fluid_mesh_pt()->boundary_node_pt(ibound,inod)->x(1);
293 if (y>y_max)
294 {
295 y_max=y;
296 }
297 if (y<y_min)
298 {
299 y_min=y;
300 }
301 }
302 double y_mid=0.5*(y_min+y_max);
303
304 // Loop over all boundaries
305 for (unsigned ibound=0;ibound<fluid_mesh_pt()->nboundary();ibound++)
306 {
307 unsigned num_nod= fluid_mesh_pt()->nboundary_node(ibound);
308 for (unsigned inod=0;inod<num_nod;inod++)
309 {
310 // Parabolic inflow at the left boundary (boundary 1)
311 if (ibound==1)
312 {
313 double y=fluid_mesh_pt()->boundary_node_pt(ibound,inod)->x(1);
314 double veloc=1.5/(y_max-y_min)*
315 (y-y_min)*(y_max-y)/((y_mid-y_min)*(y_max-y_mid));
316 fluid_mesh_pt()->boundary_node_pt(ibound,inod)->set_value(0,veloc);
317 fluid_mesh_pt()->boundary_node_pt(ibound,inod)->set_value(1,0.0);
318 }
319 // Zero flow elsewhere apart from x-velocity on outflow boundary
320 else
321 {
322 if(ibound!=2)
323 {
324 fluid_mesh_pt()->boundary_node_pt(ibound,inod)->set_value(0,0.0);
325 }
326 fluid_mesh_pt()->boundary_node_pt(ibound,inod)->set_value(1,0.0);
327 }
328 }
329 }
330} // end_of_set_boundary_conditions
331
332
333
334
335//==start_of_doc_solution=================================================
336/// Doc the solution
337//========================================================================
338template<class ELEMENT>
340{
341 ofstream some_file;
342 char filename[100];
343
344 // Number of plot points
345 unsigned npts;
346 npts=5;
347
348 // Output solution
349 sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),
350 doc_info.number());
351 some_file.open(filename);
352 fluid_mesh_pt()->output(some_file,npts);
353 some_file.close();
354} //end_of_doc_solution
355
356
357/// /////////////////////////////////////////////////////////////////////
358/// /////////////////////////////////////////////////////////////////////
359/// /////////////////////////////////////////////////////////////////////
360
361
362
363
364
365
366
367//==start_of_main======================================================
368/// Driver for unstructured fluid test problem
369//=====================================================================
370int main()
371{
372 // Label for output
373 DocInfo doc_info;
374
375 //Set the constitutive law for the pseudo-elasticity
377 new GeneralisedHookean(&Global_Physical_Variables::Nu);
378
379 //Taylor--Hood formulation
380 {
381 // Set output directory
382 doc_info.set_directory("RESLT_TH");
383
384 // Step number
385 doc_info.number()=0;
386
387 // Build the problem with TTaylorHoodElements
388 UnstructuredFluidProblem<PseudoSolidNodeUpdateElement<
389 TTaylorHoodElement<2>,
390 TPVDElement<2,3> > > problem;
391
392 // Output boundaries
393 problem.fluid_mesh_pt()->output_boundaries("RESLT_TH/boundaries.dat");
394
395 // Output the initial guess for the solution
396 problem.doc_solution(doc_info);
397 doc_info.number()++;
398
399 // Parameter study
400 double re_increment=5.0;
401 unsigned nstep=2; // 10;
402 for (unsigned i=0;i<nstep;i++)
403 {
404 // Solve the problem
405 problem.newton_solve();
406
407 // Output the solution
408 problem.doc_solution(doc_info);
409 doc_info.number()++;
410
411 // Bump up Re
412 Global_Physical_Variables::Re+=re_increment;
413 } //end_of_parameter_study
414 }
415
416
417
418 //Crouzeix--Raviart formulation
419 {
420 // Set output directory
421 doc_info.set_directory("RESLT_CR");
422
423 // Step number
424 doc_info.number()=0;
425
426 // Reset Reynolds number
428
429 // Build the problem with TCrouzeixRaviartElements
430 UnstructuredFluidProblem<PseudoSolidNodeUpdateElement<
431 TCrouzeixRaviartElement<2>,
432 TPVDBubbleEnrichedElement<2,3> > > problem;
433
434 // Output boundaries
435 problem.fluid_mesh_pt()->output_boundaries("RESLT_CR/boundaries.dat");
436
437 // Output the initial guess for the solution
438 problem.doc_solution(doc_info);
439 doc_info.number()++;
440
441 // Parameter study
442 double re_increment=5.0;
443 unsigned nstep=2; // 10;
444 for (unsigned i=0;i<nstep;i++)
445 {
446 // Solve the problem
447 problem.newton_solve();
448
449 // Output the solution
450 problem.doc_solution(doc_info);
451 doc_info.number()++;
452
453 // Bump up Re
454 Global_Physical_Variables::Re+=re_increment;
455 }
456 }
457
458
459} // end_of_main
460
461
462
463
464
465
466
467
468
469
Triangle-based mesh upgraded to become a (pseudo-) solid mesh.
ElasticTriangleMesh(const std::string &node_file_name, const std::string &element_file_name, const std::string &poly_file_name, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor:
void identify_boundaries()
Function used to identify the domain boundaries.
virtual ~ElasticTriangleMesh()
Empty Destructor.
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
~UnstructuredFluidProblem()
Destructor (empty)
void set_boundary_conditions()
Set the boundary conditions.
ElasticTriangleMesh< ELEMENT > * Fluid_mesh_pt
Fluid mesh.
void doc_solution(DocInfo &doc_info)
Doc the solution.
ElasticTriangleMesh< ELEMENT > *& fluid_mesh_pt()
Access function for the fluid mesh.
Namespace for physical parameters.
ConstitutiveLaw * Constitutive_law_pt
Constitutive law used to determine the mesh deformation.
double Nu
Pseudo-solid Poisson ratio.
int main()
///////////////////////////////////////////////////////////////////// ///////////////////////////////...