Skip to content
Permalink
Branch: master
Find file Copy path
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
128 lines (88 sloc) 4.48 KB
// This code simulates the Stokes flow of a viscous fluid (water) in a microvalve.
// Reynold's number rho*v*L/mu is very low so that the inertia terms can be neglected.
// The fluid is considered incompressible.
//
// This code can be validated with analytical solutions of a flow between two
// infinite parallel plates with a thin gap. Only the geometry has to be adapted.
#include "sparselizardbase.h"
using namespace mathop;
mesh createmesh(void);
void sparselizard(void)
{
// Define the physical regions that will be used:
int support = 1, fluid = 2, pillar = 3, membrane = 4, inlet = 7, outlet = 8;
// Create the geometry and the mesh:
mesh mymesh = createmesh();
// Write the mesh for display:
mymesh.write("microvalve.msh");
// Define the fluid boundary as the intersection between the solid and the fluid regions:
int solid = regionunion({pillar,support,membrane});
int fluidboundary = regionintersection({fluid,solid});
// Dynamic viscosity of water [Pa.s] and density [kg/m3]:
double mu = 8.9e-4, rho = 1000;
// Field v is the flow velocity. It uses nodal shape functions "h1" with two components.
// Field p is the pressure.
field v("h1xy"), p("h1");
// Force the flow velocity to 0 at the solid interface:
v.setconstraint(fluidboundary);
// Set a pressure of 100 Pa at the valve inlet and 0 at the outlet:
p.setconstraint(inlet, 100);
p.setconstraint(outlet);
// Use an order 1 interpolation for p and 2 for v on the fluid region:
p.setorder(fluid, 1);
v.setorder(fluid, 2);
// Define the weak formulation of the flow problem.
// The strong form can be found at https://en.wikipedia.org/wiki/Stokes_flow
formulation viscousflow;
viscousflow += integral(fluid, predefinedstokesflow(dof(v), tf(v), dof(p), tf(p), mu, rho) );
// Generate, solve and save:
solve(viscousflow);
// Write the p and v fields to file. Write v with an order 2 interpolation.
p.write(fluid, "p.vtk");
v.write(fluid, "v.vtk", 2);
// Output the flowrate for a unit width:
double flowrate = (normal(inlet)*v).integrate(inlet, 4);
std::cout << std::endl << "Flowrate for a unit width: " << flowrate << " m^3/s" << std::endl;
// Code validation line. Can be removed.
std::cout << (flowrate < 1.4845e-7 && flowrate > 1.4844e-7);
}
mesh createmesh(void)
{
// Give names to the physical region numbers:
int support = 1, fluid = 2, pillar = 3, membrane = 4, inlet = 7, outlet = 8;
// Define the x and y geometrical dimensions:
double y1 = 15e-6, y2 = 25e-6, y3 = 30e-6;
double x1 = -100e-6, x2 = -50e-6, x3 = -25e-6, x4 = 25e-6, x5 = 50e-6, x6 = 100e-6;
// Define the mesh finesse:
int n = 10;
shape q1("quadrangle", support, {x1,0,0, x2,0,0, x2,y1,0, x1,y1,0}, {n,n,n,n});
shape q2("quadrangle", support, {x5,0,0, x6,0,0, x6,y1,0, x5,y1,0}, {n,n,n,n});
shape q3("quadrangle", fluid, {x2,0,0, x3,0,0, x3,y1,0, x2,y1,0}, {n,n,n,n});
shape q4("quadrangle", fluid, {x4,0,0, x5,0,0, x5,y1,0, x4,y1,0}, {n,n,n,n});
shape q5("quadrangle", pillar, {x3,0,0, x4,0,0, x4,y1,0, x3,y1,0}, {n,n,n,n});
shape q6("quadrangle", support, {x1,y1,0, x2,y1,0, x2,y2,0, x1,y2,0}, {n,n,n,n});
shape q7("quadrangle", support, {x5,y1,0, x6,y1,0, x6,y2,0, x5,y2,0}, {n,n,n,n});
shape q8("quadrangle", fluid, {x2,y1,0, x3,y1,0, x3,y2,0, x2,y2,0}, {n,n,n,n});
shape q9("quadrangle", fluid, {x4,y1,0, x5,y1,0, x5,y2,0, x4,y2,0}, {n,n,n,n});
shape q10("quadrangle", fluid, {x3,y1,0, x4,y1,0, x4,y2,0, x3,y2,0}, {n,n,n,n});
shape q11("quadrangle", membrane, {x1,y2,0, x2,y2,0, x2,y3,0, x1,y3,0}, {n,n,n,n});
shape q12("quadrangle", membrane, {x5,y2,0, x6,y2,0, x6,y3,0, x5,y3,0}, {n,n,n,n});
shape q13("quadrangle", membrane, {x2,y2,0, x3,y2,0, x3,y3,0, x2,y3,0}, {n,n,n,n});
shape q14("quadrangle", membrane, {x4,y2,0, x5,y2,0, x5,y3,0, x4,y3,0}, {n,n,n,n});
shape q15("quadrangle", membrane, {x3,y2,0, x4,y2,0, x4,y3,0, x3,y3,0}, {n,n,n,n});
// Get the inlet and outlet lines:
shape l1 = q3.getsons()[0];
shape l2 = q4.getsons()[0];
l1.setphysicalregion(inlet);
l2.setphysicalregion(outlet);
// Provide to the mesh all shapes of interest:
mesh mymesh({q1,q2,q3,q4,q5,q6,q7,q8,q9,q10,q11,q12,q13,q14,q15,l1,l2});
return mymesh;
}
int main(void)
{
SlepcInitialize(0,{},0,0);
sparselizard();
SlepcFinalize();
return 0;
}
You can’t perform that action at this time.