rectangle_groove.md

Flow in a rectangle microfluidics groove

Poisson equation

We compute the velocity of a laminar flow in a rectangular, free-surface channel. The (dimensionless) velocity field u satisfies the two-dimensional Poisson equation:

with u=0 on the bottom and on the sides, and

at the surface.

Build the mesh

We first need to build a rectangular mesh. We could do this using a FreeFem++ script, like here. For a change, here, we design it in Python, and refine it with FreeFem++. Here is the full code.

Initial mesh

We first choose the rectangle's dimensions, and defines its four corners:

from pylab import *

D = 0.3
W = 1

x = array([ .5,-.5,-.5,.5 ])*W
y = array([0,0,-1,-1])*D

We then define a triangular mesh that fits into our rectangle:

import pyFreeFem as pyff

Th = pyff.TriMesh( x, y )

Th.plot_triangles()
Th.plot_nodes( labels = 'index' )

The result is:

Boundaries

At this point, the mesh Th has no boundaries. To define them, we only need to provide the list of the nodes that lie on a specific boundary. Here, a solid boundary is labelled with '1', whereas the free surface is labelled with '2'.

Th.add_boundary_edges( [1,2,3,0] )
Th.add_boundary_edges( [0,1] )

Th.plot_boundaries()
legend()

The mesh now has labelled boundaries:

Mesh refinement

Finally, we refine our mesh with the adaptmesh function of FreeFem++:

Th = pyff.adaptmesh(Th, hmax = .03, iso = 1)

The mesh now looks like so:

Solve the Poisson equation

To approximate the velocity field, we use finite elements and FreeFem++. One way to to this would be to import the build the stiffness matrix with FreeFem++, and to invert it with python, like here. This, however, is unnecessary if we solve the problem only once, and are only interested in the result.

Here, instead, we solve the entire problem with FreeFem++, and import the results with pyFreeFem.

Export the mesh to FreeFem++

We first write an edpScript object that uses our mesh:

script = pyff.InputScript( Th = Th )

print(script.get_edp())

This translates into FreeFem++ language as

mesh Th;
Th = readmesh( "/tmp/tmpf0z7l38w.msh" ) ;

So far, we haven't run FreeFem++.

Solve the finite-elements problem

We simply use a standard FreeFem++ script to do so. This one, for instance---or a slightly adapted version of it. We just need to add it to our script:

script += '''
// Fespace
fespace Vh(Th, P1); //P2 conforming triangular FEM
Vh phi, w;

// Solve
solve Poisson(phi, w)
    = int2d(Th)(
          dx(phi)*dx(w)
        + dy(phi)*dy(w)
    )
    - int2d(Th)(
          w
    )
    + on(1, phi=0)
    ;
'''

Once run in FreeFem++, this script will yield a P1 approximation of the velocity field, inconveniently named phi here.

Export the solution to Python

A last line will allow us to use the result of the FreeFem++ computation. PyFreefem can't guess the type of data we want to export, so we need to specify it.

script += pyff.OutputScript( phi = 'vector' )

We now only need to run FreeFem++ and plot the result:

u = script.get_output()['phi']

tricontourf(Th,u)

This is what we finally get:

Influence of aspect ratio on discharge

When the groove aspect ratio is large enough, the shallow-water approximation applies to most of it. Neglecting the neighborhood of its sides, the total (dimensionless) discharge reads

We now check numerically how good this approximation is.

Compute the finite-element discharge

We simply add the coresponding integral at the end of our script, and make sure we declare the result as an output:

script += '''
real discharge;
discharge = int2d(Th)(phi);
'''
script += pyff.OutputScript( discharge = 'real' )

Turn the FreeFem++ script into a function

We can call the FreeFem script multiple times with different input parameters. Here, the initial mesh will change at each call. This allow us to nest the call to the script into a function (see full code):

def get_discharge( W, D = 1 ) :

    x = array([ .5,-.5,-.5,.5 ])*W
    y = array([0,0,-1,-1])*D

    Th = pyff.TriMesh( x, y )

    Th.add_boundary_edges( [1,2,3,0] )
    Th.add_boundary_edges( [0,1] )

    FF_out = script.get_output( Th = Th )

    return FF_out['discharge']

We may then call this function at will:

W = logspace( log10(.3), 1.5, 15 )
Q = [get_discharge( W ) for W in W ]

The result looks like this:

We find that the goove needs to be about 20 times wider than it is deep for the shallow-water approximation to yield a predicion with an error of less than 10%.