Creeping Flow in a 2D Lid-Driven Cavity Tutorial

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. In the FEAScript implementation below, we set the viscosity coefficient to \(\mu = 1.0\) with the viscosity coefficient of the model (a number, or a function of \(x\) and \(y\) for a varying viscosity).

The above schematic illustrates the problem domain and outlines the associated boundary conditions. To satisfy the Ladyzhenskaya–Babuška–Brezzi (LBB) inf-sup stability condition, FEAScript uses the Taylor-Hood mixed finite element formulation. In this approach, velocity is discretized with 9-node biquadratic (Q2) elements and pressure with 4-node bilinear (Q1) elements.

Below is a demonstration of how to use the FEAScript library to solve the Stokes lid-driven cavity problem in your web browser. You only need a simple HTML page to run this example where the following code snippets should be included. First, load the required external libraries:

<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>

We should then define the problem parameters, such as the model type, the mesh configuration and the boundary conditions. This is performed using the FEAScript API directly in the HTML file:

<body>
  <!-- ...body region... -->
  <script type="module">
    import {
      FEAScriptModel,
      plotSolution,
      printVersion,
    } from "https://core.feascript.com/dist/feascript.esm.js";

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

      // 1. Create the model and select the physics
      const model = new FEAScriptModel();
      model.setModelConfig("creepingFlowScript", {
        coefficientFunctions: {
          viscosity: 1, // μ
        },
      });

      // 2. Define the mesh (creeping flow needs quadratic elements)
      model.setMeshConfig({
        meshDimension: "2D",
        elementOrder: "quadratic",
        numElementsX: 12,
        numElementsY: 6,
        maxX: 4,
        maxY: 2,
      });

      // 3. Apply the boundary conditions
      model.addBoundaryCondition("0", ["constantVelocity", 0, 0]); // Bottom wall
      model.addBoundaryCondition("1", ["constantVelocity", 0, 0]); // Left wall
      model.addBoundaryCondition("2", ["constantVelocity", 1, 0]); // Moving lid
      model.addBoundaryCondition("3", ["constantVelocity", 0, 0]); // Right wall

      // 4. Solve
      model.setSolverMethod("lusolve");
      const result = model.solve();

      // 5. Plot the results
      plotSolution(model, result, "contour", "velocityMagnitudeCanvas");
      plotSolution(model, result, "contour", "uVelocityCanvas", { field: "u" });
      plotSolution(model, result, "contour", "vVelocityCanvas", { field: "v" });
    });
  </script>
  <!-- ...rest of body region... -->
</body>

In the boundary condition definition, the numbers at the left side (from "0" to "3") indicate the boundaries of the geometry (the mesh generator of FEAScript has assigned numbers to the boundaries, starting from the bottom boundary and proceeding clockwise). The "constantVelocity" condition sets the velocity components \(u\) and \(v\) at the boundary nodes (the second and third arguments respectively).

The mesh defined in step 2 has 12 × 6 quadratic elements. The figure also shows the boundary numbering used in the boundary conditions:

Since the Stokes solver uses a mixed formulation, the result of solve() also contains the named fields velocityMagnitude, u, v and p, each with one value per mesh node (the pressure, solved at the corner nodes, is interpolated linearly to the remaining nodes). plotSolution plots the velocity magnitude by default, and any other field with the field option.

After solving the case, the results are shown in 2D contour plots. To visualize them, include HTML containers where the plots will render:

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

Below are the 2D contour plots of the computed velocity magnitude \(|\mathbf{u}| = \sqrt{u^2+v^2}\), horizontal velocity component \(u\), and vertical velocity component \(v\). These plots are generated in real time using FEAScript.

• FEAScript Website
• All Tutorials