Creeping flow through a 2D pillar array.

The steady Stokes equations for an incompressible Newtonian fluid with no body forces are: \(\nabla \cdot \mathbf{u} = 0\) and \(-\mu \nabla^{2} \mathbf{u} + \nabla p = \mathbf{0}\), where \(\mathbf{u} = (u,v)\) is the velocity vector, \(p\) is the pressure and \(\mu\) is the dynamic viscosity. They hold when the Reynolds number is small, so that viscous forces dominate inertia, as in most microchannels. We set \(\mu = 1\) and work in nondimensional units.

The fluid enters at the left with a uniform velocity \(u = 1, v = 0\) and leaves at the right, where the traction is zero (stress-free outlet). The velocity is zero on the channel walls and on the pillars (no-slip). The pillars are placed in five columns at \(x = 2\) to \(6\), alternating two pillars at \(y = 0.5\) and \(1.5\) and one pillar at \(y = 1\), so the fluid has to weave between them.

As in the lid-driven cavity tutorial, FEAScript satisfies the Ladyzhenskaya–Babuška–Brezzi (LBB) inf-sup condition with Taylor-Hood elements. On triangles, the velocity uses 6-node quadratic (P2) elements and the pressure 3-node linear (P1) elements on the corner nodes.

The mesh is generated with Gmsh from the geometry file pillar.geo. It subtracts the pillars from the channel, finds each boundary by its position, and puts the boundaries in four physical groups:

// 4. Physical groups: FEAScript keys the boundary conditions by these names
Physical Curve("inlet", 1) = {inlet()};          // FEAScript boundary "inlet"
Physical Curve("walls", 2) = {bottom(), top()};  // FEAScript boundary "walls"
Physical Curve("outlet", 3) = {outlet()};        // FEAScript boundary "outlet"
Physical Curve("pillars", 4) = {pillars()};      // FEAScript boundary "pillars"
Physical Surface("domain", 5) = {10};            // Domain

// 5. Mesh: quadratic triangles (P2), refined on the pillars
MeshSize{ PointsOf{ Surface{10}; } } = lc;
MeshSize{ PointsOf{ Curve{pillars()}; } } = lp;
Mesh.MeshSizeFromCurvature = 12; // Elements per 2π of curvature on the pillars
Mesh.ElementOrder = 2;

Run gmsh pillar.geo -2 -format msh41 -o pillar_triang_quadratic.msh to write the mesh in the .msh version 4.1 format that FEAScript reads, or download pillar_triang_quadratic.msh. FEAScript keys each boundary by the name of its physical group, and these names are used in the boundary conditions below.

The model runs in a simple HTML page. First, load Plotly in the <head>:

<head>
  <!-- ...head region... -->
  <script src="https://cdnjs.cloudflare.com/ajax/libs/plotly.js/2.35.3/plotly.min.js"></script>
  <!-- ...rest of head region... -->
</head>

Then fetch the mesh, define the model and solve. The mesh file must be served from the same site as the page (or from a server that allows cross-origin requests):

<body>
  <!-- ...body region... -->
  <script type="module">
    import {
      FEAScriptWorker,
      importGmshMesh,
      plotInterpolatedSolution,
      version,
    } from "https://core.feascript.com/dist/feascript.esm.js";

    window.addEventListener("DOMContentLoaded", async () => {
      console.log("FEAScript version:", version);

      // 1. Create the model in a web worker, so the solve does not block the page, and select the physics
      const model = new FEAScriptWorker();
      await model.setModelConfig("creepingFlow", {
        solver: { method: "lusolve" },
        coefficientFunctions: {
          viscosity: 1, // μ
        },
      });

      // 2. Load the Gmsh mesh (P2 triangles)
      const response = await fetch("./pillar_triang_quadratic.msh");
      const meshFile = new File([await response.text()], "pillar_triang_quadratic.msh");
      await model.setMesh(await importGmshMesh(meshFile));

      // 3. Apply the boundary conditions, keyed by the names of the Gmsh physical groups
      await model.addBoundaryCondition("inlet", { type: "constantVelocity", u: 1, v: 0 });   // Inlet
      await model.addBoundaryCondition("walls", { type: "constantVelocity", u: 0, v: 0 });   // Channel walls
      await model.addBoundaryCondition("outlet", { type: "stressFree" });                    // Outlet
      await model.addBoundaryCondition("pillars", { type: "constantVelocity", u: 0, v: 0 }); // Pillars

      // 4. Solve in the worker
      const result = await model.solve();

      // 5. Plot the velocity magnitude, with the model information of the worker
      const modelInfo = await model.getModelInfo();
      plotInterpolatedSolution(modelInfo, result, "contour", "velocityMagnitudeCanvas");

      // 6. Terminate the worker
      model.terminate();
    });
  </script>
  <!-- ...rest of body region... -->
</body>

The "stressFree" condition at the outlet is the natural boundary condition of the formulation, so it adds nothing to the system; it is listed so that every boundary has a condition. Because the outlet is open, the pressure level is fixed by the outlet and no pressure node needs to be pinned. plotInterpolatedSolution evaluates the quadratic shape functions inside each triangle, so the plot shows the curved contours of the P2 solution. The solve runs in a web worker (FEAScriptWorker), so the page stays responsive while it runs. Add the container where the plot will render:

<body>
  <!-- ...body region... -->
  <div id="velocityMagnitudeCanvas"></div>
  <!-- ...rest of body region... -->
</body>

Below is the 2D contour plot of the velocity magnitude \(|\mathbf{u}| = \sqrt{u^2+v^2}\), computed in your browser when the page loads. The velocity components and the pressure are in result.fields (u, v and p), and plot the same way, e.g. with { field: "p" }. A Node.js version of this tutorial, which also prints the flow rates, is in the example directory.

The flow accelerates through the gaps between the pillars. The fastest jets, with a velocity of about 5.2 (five times the inflow), pass through the 0.4-wide gaps between the paired pillars at \(x = 2\), \(4\) and \(6\). The vertical velocity changes sign from one column to the next as the flow weaves around the pillars. The pressure falls from about 545 near the inlet to zero at the outlet, and almost all of this drop is across the pillars: this is the resistance of the array, the quantity a microfluidic designer would size a pump for. With \(\mu = 1\) and an inflow velocity of 1, the pressure is in units of \(\mu U / \ell\), with \(\ell\) the length unit of the geometry.

There is no analytical solution for this geometry, but two checks show whether the result is sound, and one point needs care:

• Mass conservation. The flow rate through the outlet should equal the flow rate through the inlet. Summed from the boundary nodes, they are 1.875 and 1.909, within 2%. Taylor-Hood elements conserve mass only in the weak (integral) sense, so a small difference is expected, and it shrinks as the mesh is refined. The inlet rate is below 2 because the corner nodes belong to the no-slip walls.
• Mesh convergence. A finer mesh (element size 0.15 in the channel and 0.05 on the pillars, 11,265 nodes) raises the largest velocity from 5.18 to 5.28, a change of 2%. Refine the mesh in pillar.geo to see how the result changes.
• Corner pressure. At the two inlet corners the velocity jumps from 1 to 0, and the pressure is singular there: its peak (564 on this mesh) grows as the mesh is refined instead of converging. Read the pressure drop away from the corners. Docs: Pressure level.