---
title: "Laminar Flow"
language: "en"
type: "Monograph"
summary: "The analysis and behavior of fluids are of fundamental importance in science and engineering. This monograph gives an introduction to modeling fluids with partial differential equations. Equations and boundary conditions that are relevant for performing fluid mechanics analysis are derived and explained. Fluids are substances in gaseous or liquid phase, like air or water, respectively. The term fluid comprises both liquids and gases. The opposite of fluids is solids, and the distinguishing attribute is that fluids cannot withstand shear stress when at rest. A fluid will take the shape of a container. The borderline between what is a fluid and a solid is not always clear. For example, asphalt is usually considered to be solid but on a longer timescale, it starts to behave as a fluid. A distinction is typically made between the following fluid regimes:"
keywords: 
- flow
- Laminar flow
- laminar
- Navier Stokes
- Stokes
- newtonian flow
- non newtonian flow
- compressible
- incompressible
- mach number
- Reynolds number
- viscosity
- apparent viscosity
- Boussinesq
- Boussinesq Approximation
- Navier Stokes Fourier
- Rayleigh–Benard
- Rayleigh–Bénard
- Rayleigh–Bénard convection
- convection
- axisymmetric flow
- power law
- carreau
- Bingham
- Papanasasiou
- Herschel
- Bulkley
- swirl
- swirling flow
canonical_url: "https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.html"
source: "Wolfram Language Documentation"
related_tutorials: 
  - 
    title: "PDEModels Overview"
    link: "https://reference.wolfram.com/language/PDEModels/tutorial/PDEModelsOverview.en.md"
  - 
    title: "Fluid Dynamics Model Verification Tests"
    link: "https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/FluidDynamicsVerificationTests.en.md"
---
# Laminar Flow	[image]

Contents

* [Introduction](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1723663506)

* [Overview Example and Analysis Types](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1313601044)

	* [Stationary Analysis](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1172177949)

	* [Post-Processing](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1045542984)

		* [Common visualization techniques](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1348458401)

	* [Parametric Analysis](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1014914524)

	* [Time-Dependent Analysis](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#181369525)

* [Equations](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1212896928)

	* [Cylindrical Coordinates](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1593586754)

		* [Axisymmetric models](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1212358490)

		* [Swirling flow](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1483426158)

* [Newtonian versus Non-Newtonian Flow](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1789304502)

	* [Power Law](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1686395575)

		* [Power law example](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#433589366)

	* [Carreau](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1432993996)

	* [Cross](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#137409490)

	* [Bingham–Papanastasiou](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#787419883)

	* [Herschel–Bulkley–Papanastasiou](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#270042841)

	* [Custom Apparent Viscosity Model](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1025599629)

	* [Custom Viscous Stress Tensor](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1275114574)

* [Computing viscosity and stress](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#166104088)

* [Viscoelastic Flow](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#84842597)

* [Energy Transport—Nonisothermal Flow](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#2054167427)

	* [Boussinesq Approximation](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#449698672)

	* [Rayleigh–Bénard Convection](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1601703969)

* [Flow Boundary Conditions](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1856021341)

	* [Inflow Conditions](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#308138130)

	* [Outflow Conditions](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#54168846)

	* [Wall Conditions](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#2018602509)

		* [No-slip](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1696064518)

		* [No-slip example](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1792824947)

		* [Slip](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1768541953)

	* [Traction](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#683422966)

		* [Traction example](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#425609897)

	* [Initial Conditions](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#2043187725)

* [Convergence of Fluid Flow Models](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1014610487)

* [Initial Seeding](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1443519582)

	* [Iterative Stepping](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1571490221)

* [Equation Modification](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1769874127)

	* [Stokes Equation](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1148220986)

* [References](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#480587635)

Introduction

The analysis and behavior of fluids are of fundamental importance in science and engineering. This monograph gives an introduction to modeling fluids with partial differential equations. Equations and boundary conditions that are relevant for performing fluid mechanics analysis are derived and explained.

Fluids are substances in gaseous or liquid phase, like air or water, respectively. The term "fluid" comprises both liquids and gases. The opposite of fluids is solids, and the distinguishing attribute is that fluids cannot withstand shear stress when at rest. A fluid will take the shape of a container. The borderline between what is a fluid and a solid is not always clear. For example, asphalt is usually considered to be solid but on a longer timescale, it starts to behave as a fluid.

A distinction is typically made between the following fluid regimes:

* Creep flow

* Laminar flow

* Turbulent flow

The difference between these three flow regimes is given by the Reynolds number $ℛℯ$, which is the ratio of inertia forces versus viscous forces. It is used to predict whether or not turbulence will appear in a fluid flow. Turbulent flow is characterized by a high Reynolds number, while laminar flow occurs when the Reynolds number is small. The term creep flow is reserved for flows with a very small Reynolds number. The concept of the Reynolds number will be introduced more formally later, and what high and low Reynolds numbers mean will also be discussed.

This monograph focuses on laminar flow.

Modeling fluids with partial differential equations (PDEs) is not the only way to model fluids. Other techniques include setting up ordinary differential equations (ODEs). This approach is followed by [Wolfram System Modeler](https://reference.wolfram.com/system-modeler/). Roughly speaking, the system modeler approach is more suitable for large systems of flows interacting, like flow though pipes, while the partial differential equation approach is more suitable for a fine-grained analysis of a specific fluid flow behavior. In some cases, it is beneficial to use a combination of the two approaches.

The approach taken here is that in an introductory section, a single fluid flow example, a driven cavity example, is used to introduce various fluid flow analysis types and the functionality available. This will be followed by a more theoretical explanation of the underlying ideas and concepts. The theoretical background is much easier to understand once an intuition for the various analysis types exists. After that, the available boundary conditions are discussed.

The goal of a fluid flow analysis is to find the vector-valued fluid flow velocities of a fluid under constraints. Another goal might be to find the scalar pressure field or the scalar density field, depending on whether an incompressible or a compressible fluid is considered. The analysis and interpretation of the fluid behavior is useful to create a better-quality engineering design of the flow behavior in question. For example, through some parts of an object there might not be sufficient flow. These areas can then be identified and improved upon.

A fluid flow analysis is typically done in stages. First, for the body in which the fluid flows, a geometric model needs to be created. The geometric model is typically created within a computer-aided design (CAD) process. CAD models can either be imported or created in product. To import geometries, common file formats like [`DXF`](https://reference.wolfram.com/language/ref/format/DXF.en.md), [`STL`](https://reference.wolfram.com/language/ref/format/STL.en.md) or [STEP](https://reference.wolfram.com/language/OpenCascadeLink/tutorial/UsingOpenCascadeLink.en.md#977262051) are supported. These geometries can be imported with ``Import``. The alternative is to create the geometrical models in-product, for example, by using [OpenCascadeLink](https://reference.wolfram.com/language/OpenCascadeLink/tutorial/UsingOpenCascadeLink.en.md). Once the geometric model is made available, some thought needs to be put into what type of analysis is to be performed. Currently supported analysis types are static and time-dependent fluid flow analysis. The next step is the setup of boundary conditions and constraints. Material properties of the fluid in question further specify the PDE model. Once the PDE model is fully specified, the subsequent finite element method will then compute the desired quantities under investigation. These quantities are then post-processed, either by visualizing them or some derived quantities are computed. This notebook shows the necessary steps for everything except the CAD model generation.

The modeling process as such results in a system of partial differential equations (PDEs) that can be solved with ``NDSolve`` and ``ParametricNDSolve``.

Many of the animations of the simulation results shown in this tutorial are generated with a call to ``Rasterize``. This is to reduce the disk space required. The downside is that the visual quality of the animations will not be as crisp as without it. To obtain high-quality graphics, remove or comment out the call to ``Rasterize``.

To get high-fidelity visualizations, comment out the rasterization process.

```wl
(*frames=Rasterize[#1,"Image",ImageResolution->35]&/@frames;*)
```

The symbols and corresponding units used throughout this tutorial are summarized in the [Nomenclature](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#658755935) section.

Load the Finite Element package and set the ``\$HistoryLength`` :

```wl
In[2]:=
Needs["NDSolve`FEM`"]
$HistoryLength = 0;
```

Overview Example and Analysis Types

This monograph starts with an overview example, which focuses on the workflow more than the actual equations and theory. This example considers the flow of glycerol in a narrowing tube.

Creating a fluid flow model always comprises the same steps:

* Geometry creation

* Material selection

* Boundary condition specification

* Mesh generation

* Solving the partial differential equations

* Post-processing the solution

Set up and visualize the geometry:

```wl
In[18]:=
Ω = RegionUnion[Rectangle[{0, 0}, {1, 1 / 2}], Rectangle[{1, 1 / 10}, {2, 2 / 5}]];
RegionPlot[Ω, AspectRatio -> Automatic]

Out[19]= [image]
```

Set up model variables:

```wl
In[20]:= vars = {{u[x, y], v[x, y], p[x, y]}, {x, y}};
```

The dependent velocity variables $u$, $v$ represent the displacement in the $x$ and $y$ directions, respectively. $p$ is the pressure.

The fluid is glycerol:

```wl
In[6]:= pars = <|"Material" -> Entity["Chemical", "Glycerol"]|>;
```

Here a chemical entity was used for the material specification. Alternatively, it is also possible to specify material properties.

Specify fluid parameters through material properties:

```wl
In[6]:= <|"DynamicViscosity" -> Quantity[0.934, "Pascals"*"Seconds"], "MassDensity" -> Quantity[1.25, "Grams"/"Centimeters"^3]|>;
```

Set up the fluid flow PDE component:

```wl
In[7]:= MatrixForm[op = FluidFlowPDEComponent[vars, pars]]

Out[7]//MatrixForm=
(⁠|     |
| -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
| {1250. u[x, y], 1250. v[x, y]}. ... . u[x, y], 1250. v[x, y]}.Inactive[Grad][v[x, y], {x, y}] + Inactive[Div][{0, p[x, y]}, {x, y ... vative[0, 1][u][x, y] + Derivative[1, 0][v][x, y]),    -1.868*Derivative[0, 1][v][x, y]}, {x, y}] |
| 1250. v^(0, 1)[x, y] + 1250. u^(1, 0)[x, y] |⁠)
```

The inflow velocity in the $x$ direction is prescribed through boundary conditions. This is done by specifying an inflow profile for the $u$ dependent variable on the left-hand side of the region. On the bottom and upper parts of the region, a wall boundary condition is specified. This is done by setting both the $u$ and $v$ components on those parts to $0$. An outflow condition is set on the right-hand side of the domain. Here the pressure $p$ is arbitrarily set to $0$. In general, each dependent variable should either have a [Dirichlet condition](https://reference.wolfram.com/language/ref/DirichletCondition.en.md#1759349358) or a [Robin-type Neumann](https://reference.wolfram.com/language/ref/NeumannValue.en.md#1397642515) specified. A pure Neumann value like Neumann $0$ is not sufficient. In that case, the solution to the equation can only be correct up to a constant.

Set up boundary conditions:

```wl
In[8]:=
bcs = {DirichletCondition[
	{u[x, y] == 4 * 0.3 * y * (0.5 - y) / (0.41) ^ 2, v[x, y] == 0.}, x == 0.], DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, 0 < x < 2], 
	DirichletCondition[p[x, y] == 0., x == 2]};
```

A stable solution can be found if the velocities are interpolated with a higher order than the pressure. ``NDSolve`` allows an interpolation order for each dependent variable to be specified. Most of the time, specifying the ``"InterpolationOrder"`` is not necessary, but this is a standard procedure when using the finite element method for fluid flow, and this is how it can be done in the Wolfram Language. This setup corresponds to what is called P2P1- or Q2Q1-type elements, which are also known as Taylor–Hood elements. Leaving out the ``"InterpolationOrder"`` will result in using a second-order interpolation for all variables, which can lead to instabilities.

Solve the PDE:

```wl
In[12]:= {uVel, vVel, pressure} = NDSolveValue[{op == {0, 0, 0}, bcs}, {u, v, p}, {x, y}∈Ω, Method -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}]

Out[12]=
{InterpolatingFunction[{{0., 2.0000000000000284}, {0., 0.5000000000000071}}, 
 {5, 4225, 0, {1052, 0}, {3, 0}, 0, 0, 0, 0, Indeterminate & , {}, {}, False}, 
 {NDSolve`FEM`ElementMesh[CompressedData["«16627»"], {NDSolve`FEM`TriangleElement[Compress ...  3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 11, 16, 8, 8, 4, 4, 13, 15, 6, 6, 6, 6, 
      6, 6, 6, 14, 2, 2, 2, 2, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 5, 5, 5, 5, 5, 5, 5, 5, 5, 
      5, 5, 5, 5, 5, 7, 5}]}]}, CompressedData["«3021»"], {Automatic}]}
```

In many cases, there is no need to generate a mesh; passing the geometry to the solver is sufficient.

Visualize the fluid velocity:

```wl
In[26]:= VectorPlot[{uVel[x, y], vVel[x, y]}, {x, y}∈Ω, AspectRatio -> Automatic]

Out[26]= [image]
```

Visualize the contours of the pressure distribution:

```wl
In[27]:= ContourPlot[pressure[x, y], {x, y}∈Ω, AspectRatio -> Automatic]

Out[27]= [image]
```

This gives an overview of the workflow. The next sections describe the steps in more detail.

Stationary Analysis

In this example, the driven cavity example, a classic benchmark in fluid dynamics [6] is investigated. Here a rectangular region is filled with a fluid. On the top there is a mechanism that drives the fluid with a fixed flow velocity in the positive $x$ direction, much like a conveyor belt would. The remaining sides are walls and have a no-slip condition. No-slip means that the velocity at the wall boundary is $0$ in both the $x$ and $y$ directions. In the lower-left corner, there is a reference pressure condition.

[image]

In a rectangular region filled with a fluid, the top is driven by a velocity $Subscript[u, 0]$ in the $x$ direction. All other edges are walls, where a no-slip, $u = v = 0$ boundary condition is set.

This driven cavity will be modeled as a time-independent, stationary model. The dependent velocity variables $u$ and $v$ and the dependent pressure variable $p$ are set up. The independent variables are the spatial variables $x$ and $y$.

Set up the stationary model variables:

```wl
In[5]:= vars = {{u[x, y], v[x, y], p[x, y]}, {x, y}};
```

Create a rectangular region:

```wl
In[6]:= Ω = Rectangle[{0, 0}, {1, 1}];
```

Standard texts present this example with various Reynolds numbers $ℛℯ$. With everything else remaining the same, the flow pattern changes, depending only on the values of the Reynolds number $ℛℯ$.

Set up the fluid flow parameters ``pars`` with $ℛℯ = 1000$ :

```wl
In[26]:= pars = <|"ReynoldsNumber" -> 1000|>;
```

Set up the fluid flow PDE:

```wl
In[27]:= pde = FluidFlowPDEComponent[vars, pars] == {0, 0, 0}

Out[27]=
{{u[x, y], v[x, y]}.Inactive[Grad][u[x, y], {x, y}] + Inactive[Div][{p[x, y], 0}, {x, y}] + Inactive[Div][{(-(1/500))*Derivative[1, 0][u][x, y], 
  (-Derivative[0, 1][u][x, y] - Derivative[1, 0][v][x, y])/1000}, {x, y}], {u[x, y], v[x, y]}.Inactive[Grad][v[x, y], {x, y}] + Inactive[Div][{0, p[x, y]}, {x, y}] + Inactive[Div][{(-Derivative[0, 1][u][x, y] - Derivative[1, 0][v][x, y])/1000, 
  (-(1/500))*Derivative[0, 1][v][x, y]}, {x, y}], v^(0, 1)[x, y] + u^(1, 0)[x, y]} == {0, 0, 0}
```

Set up the top boundary condition:

```wl
In[28]:= top = DirichletCondition[{u[x, y] == 1, v[x, y] == 0}, y == 1]

Out[28]= DirichletCondition[{u[x, y] == 1, v[x, y] == 0}, y == 1]
```

Set up wall boundary conditions:

```wl
In[29]:= wall = DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, y < 1]

Out[29]= DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, y < 1]
```

A reference pressure point should be applied in a "calm" area, or in other words, far away from the interesting behavior.

Pressure point boundary condition:

```wl
In[30]:= referencePressure = DirichletCondition[p[x, y] == 0, x == 0 && y == 0]

Out[30]= DirichletCondition[p[x, y] == 0, x == 0 && y == 0]
```

Combine the boundary conditions:

```wl
In[31]:= bcs = {top, wall, referencePressure};
```

Solve the fluid flow PDE:

```wl
In[32]:= {uVel, vVel, pressure} = NDSolveValue[{pde, bcs}, vars[[1]], {x, y}∈Ω, Method -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}]

Out[32]=
{InterpolatingFunction[{{0., 1.}, {0., 1.}}, {5, 4225, 0, {1281, 0}, {3, 0}, 0, 0, 0, 0, 
  Indeterminate & , {}, {}, False}, 
 {NDSolve`FEM`ElementMesh[CompressedData["«4656»"], 
   {NDSolve`FEM`QuadElement[CompressedData["«6192»"]]}, 
   {NDSolve ...  7, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 
      3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 6, 4, 4, 4, 
      4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 8}]}]}, 
 CompressedData["«4796»"], {Automatic}][x, y]}
```

The solution returned constrains the $u$ velocity profile, which is the velocity profile in the $x$ direction; the $v$ velocity profile, which is the velocity profile in the $y$ direction; and a pressure distribution $p$. There are various ways to visualize the solution.

The ``"InterpolationOrder"`` method option given to ``NDSolve`` specifies that both velocities are interpolated with second order and that the pressure is interpolated with first order. The ``"InterpolationOrder"`` option specification is equivalent to what is known in the finite element world as a Taylor–Hood element, and it helps stabilize the solution. Unfortunately, there is no known way to automatically set this option for specific PDEs, and it is recommended that the option be set for fluid flow PDEs that are based on the Navier–Stokes equations, like the ``FluidFlowPDEComponent``.

Post-Processing

Common visualization techniques

The computed pressure is a scalar field, and as such, can be visualized like any other scalar field. Contour or density plots are a good way to visualize these.

Visualize the pressure field with a contour plot:

```wl
In[33]:= ContourPlot[pressure, {x, y}∈Ω, PlotRange -> All]

Out[33]= [image]
```

The velocity field is a vector-valued field. Here, vector, stream or line integral convolution plots are a good choice for visualization.

Visualize the velocity field as a vector plot:

```wl
In[23]:= VectorPlot[{uVel, vVel}, {x, y}∈Ω]

Out[23]= [image]
```

Visualize the velocity field as a stream plot:

```wl
In[18]:= StreamPlot[{uVel, vVel}, {x, y}∈Ω]

Out[18]= [image]
```

Visualize the velocity field as a line integral convolution plot:

```wl
In[19]:= LineIntegralConvolutionPlot[{uVel, vVel}, {x, y}∈Ω, ColorFunction -> "RoseColors", LineIntegralConvolutionScale -> 1]

Out[19]= [image]
```

In two dimensions, a vorticity plot can also provide useful information.

Create a function to compute the vorticity of the flow field:

```wl
In[20]:= vorticity = Function[{x, y}, Evaluate[Curl[{uVel, vVel}, {x, y}]]];
```

Visualize the contours in the range of $[-10, 10]$ of the vorticity:

```wl
In[21]:=
cc = Subdivide[-10, 10, 50];
colors = ColorData["TemperatureMap"][Rescale[#, MinMax[cc]]]& /@ cc;
Show[
	RegionPlot[Ω], 
	ContourPlot[vorticity[x, y], {x, y}∈Ω, PlotRange -> All, ContourShading -> None, Contours -> cc, ContourStyle -> colors, Frame -> None]]

Out[23]= [image]
```

Parametric Analysis

As a next step, it would be interesting to compare the flow field for various Reynolds numbers. ``ParametricNDSolve`` works just like ``NDSolve``, with the addition that parameters can be specified. ``ParametricNDSolve`` then returns a function that, when fed numerical values for the parameters, computes the numerical solution. In this case, the Reynolds number is the parameter.

Set a symbolic Reynolds number ``re`` :

```wl
In[24]:= pars = <|"ReynoldsNumber" -> re|>;
```

Set up the fluid flow PDE:

```wl
In[25]:= pde = FluidFlowPDEComponent[vars, pars] == {0, 0, 0}

Out[25]=
{{u[x, y], v[x, y]}.Inactive[Grad][u[x, y], {x, y}] + Inactive[Div][{p[x, y], 0}, {x, y}] + Inactive[Div][{-((2*Derivative[1, 0][u][x, y])/re), 
  -((Derivative[0, 1][u][x, y] + Derivative[1, 0][v][x, y])/re)}, {x, y}], {u[x, y], v[x, y]}.Inactive[Grad][v[x, y], {x, y}] + Inactive[Div][{0, p[x, y]}, {x, y}] + Inactive[Div][{-((Derivative[0, 1][u][x, y] + Derivative[1, 0][v][x, y])/re), 
  -((2*Derivative[0, 1][v][x, y])/re)}, {x, y}], v^(0, 1)[x, y] + u^(1, 0)[x, y]} == {0, 0, 0}
```

Set up a ``ParametricFunction`` for various Reynolds numbers:

```wl
In[26]:= pfun = ParametricNDSolveValue[{pde, bcs}, vars[[1]], {x, y}∈Ω, re, Method -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}]

Out[26]= ParametricFunction[<>]
```

Compute the solution for several Reynolds numbers:

```wl
In[27]:= solutions = pfun /@ {1, 10, 100, 1000, 2000, 4000};
```

In the previous example, where ``NDSolve`` was used, the solution was stored in a list of the form ``{uVel, vVel, pressure}``. Now, the solution of all simulation runs is stored in the variable ``solutions``. Each entry in the solutions vector now contains an interpolation function for the $u$ field, the $v$ field and the pressure $p$.

Visualize the velocity fields of each of the solutions:

```wl
In[28]:= StreamPlot[#, {x, y}∈Ω]& /@ solutions[[All, {1, 2}]]

Out[28]= [image]
```

Too large a Reynolds number will eventually prohibit convergence of the nonlinear PDE, as the PDE then becomes too strongly nonlinear. A failure in convergence is typically announced by a message from ``FindRoot``, as ``FindRoot`` is the function that eventually solves the nonlinear equations within ``NDSolve``. Not all hope is lost then, however. There are various measures that can be taken to extend the range in which a model converges.

Not all Reynolds numbers will converge right away:

```wl
In[29]:= pfun[7000]
```

FindRoot::dfmin: The minimal damping factor of 1/10000 has been reached.

ParametricNDSolveValue::fempsf: PDESolve could not find a solution.

```wl
Out[29]= ParametricFunction[<>][7000]
```

There are various ways to circumvent the loss of convergence, which are discussed in the section [Convergence of Fluid Flow Models](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1014610487). Here a closer look is taken at the mesh used for the analysis. In the above example, a mesh was generated automatically. One option is to make use of a mesh that is more appropriate for the problem at hand. The specifics of what "more appropriate" means depend on the problem at hand. For this example, note that the gradients of the $u$- and $v$-velocity fields are largest at the top and at the right wall.

Find the positions of the largest gradients of the $u$- and $v$-velocity fields for the last computed Reynolds number:

```wl
In[30]:= GraphicsRow[ContourPlot[Evaluate[Sqrt[Total[Grad[#, {x, y}] ^ 2]]], {x, y}∈Ω, PlotRange -> All]& /@ solutions[[-1, {1, 2}]]]

Out[30]= [image]
```

The generated mesh should reflect the steep gradients by having a finer resolution at steep gradients. Since the flow gradient is higher at the walls, a possible way to improve on that is to make use of a graded mesh. This can be done with the function ``ToGradedMesh``. This allows for more refinement at the boundary, while not increasing the overall number of elements too much. Other refinement options like using a ``MeshRefinementFunction`` are also possible. More information about the mesh generation process can be found in the [ElementMesh generation tutorial](https://reference.wolfram.com/language/FEMDocumentation/tutorial/ElementMeshCreation.en.md).

Generate a mesh with refinement at the walls and top:

```wl
In[3]:=
m1 = ToGradedMesh[Line[{{0}, {1}}], <|"Alignment" -> "BothEnds", "ElementCount" -> 75|>];
mesh = ElementMeshRegionProduct[m1, m1]

Out[4]= ElementMesh[{{0., 1.}, {0., 1.}}, {QuadElement[<5625>]}]
```

Visualize the mesh:

```wl
In[5]:= mesh["Wireframe"]

Out[5]= [image]
```

With the graded mesh, the overall number of elements is reduced, and at the same time, the element density is increased near the walls.

Create a new parametric function, now using the constructed mesh:

```wl
In[34]:= pfun = ParametricNDSolveValue[{pde, bcs}, vars[[1]], {x, y}∈mesh, re, Method -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}]

Out[34]= ParametricFunction[<>]
```

Compute solutions for various Reynolds numbers:

```wl
In[35]:=
reynoldsNumbers = {100, 400, 1000, 3200, 5000, 7500, 10000};
solutions = pfun /@ reynoldsNumbers;
```

Extract a few mesh points to be used as ``StreamPoints`` :

```wl
In[37]:= pts = Take[Flatten[m1["Coordinates"]], {1, -1, 10}];
```

Visualize the solutions at the various Reynolds numbers:

```wl
In[40]:=
frames = MapThread[(StreamPlot[#1, {x, y}∈mesh, StreamPoints -> {Tuples[pts, 2], Automatic, 1}, StreamMarkers -> None, RegionFillingStyle -> None, PlotLabel -> "Reynolds Number: " <> ToString[#2], Frame -> None] /. Arrow -> Line)&, {solutions[[All, {1, 2}]], reynoldsNumbers}];frames = Rasterize[#, "Image", ImageResolution -> 120]& /@ frames;
ListAnimate[frames]

Out[41]= DynamicModule[«8»]
```

See [this note](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1211049090) about improving the visual quality of the animation.

Since the driven cavity is a common fluid dynamics benchmark problem [6], literature data exists that can be used to compare numerical simulations.

Import the benchmark results:

```wl
In[42]:=
uyData = Import[FileNameJoin[{$InstallationDirectory, "SystemFiles", "Components", "PDEModels", "SupportFiles", "LidDrivenCavityBenchmarkUY.wdx"}]];
vxData = Import[FileNameJoin[{$InstallationDirectory, "SystemFiles", "Components", "PDEModels", "SupportFiles", "LidDrivenCavityBenchmarkVX.wdx"}]];
```

Compare the computed solution at $x = 1 / 2$ with the benchmark results for various Reynolds numbers:

```wl
In[44]:=
Show[
	ListPlot[uyData], 
	Plot[Evaluate[Through[solutions[[All, 1, 0]][0.5, y]]], {y, 0, 1}, PlotLegends -> LineLegend[reynoldsNumbers, LegendLabel -> "TraditionalForm\`ℛℯ"]], ImageSize -> Medium]

Out[44]= [image]
```

Compare the computed solution at $y = 1 / 2$ with the benchmark results for various Reynolds numbers:

```wl
In[45]:=
Show[
	ListPlot[vxData], 
	Plot[Evaluate[Through[solutions[[All, 2, 0]][x, 0.5]]], {x, 0, 1}, PlotLegends -> LineLegend[reynoldsNumbers, LegendLabel -> "TraditionalForm\`ℛℯ"]], ImageSize -> Medium]

Out[45]= [image]
```

The solution matches the literature values.

Inspect the velocity field at the lower-right corner for Reynolds number $10000$ :

```wl
In[46]:= StreamPlot[solutions[[-1, {1, 2}]], {x, 0.6, 1}, {y, 0, 0.3}, AspectRatio -> Automatic, PlotRange -> All, StreamPoints -> Fine]

Out[46]= [image]
```

Inspect the velocity field at the lower-left corner for Reynolds number $10000$ :

```wl
In[47]:= StreamPlot[solutions[[-1, {1, 2}]], {x, 0, 0.4}, {y, 0, 0.3}, AspectRatio -> Automatic, StreamPoints -> Fine]

Out[47]= [image]
```

Inspect the velocity field at the top-left corner for Reynolds number $10000$ :

```wl
In[48]:= StreamPlot[solutions[[-1, {1, 2}]], {x, 0, 0.4}, {y, 0.7, 1}, AspectRatio -> Automatic, StreamPoints -> Fine]

Out[48]= [image]
```

For computing even higher Reynolds numbers of the driven cavity example, see the section [Convergence of Fluid Flow Models](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1014610487).

Time-Dependent Analysis

As a next workflow example, the time-dependent driven cavity example is solved. Now the fluid is initially at rest and the driving mechanism will ramp up with time. Everything else remains the same.

The time-dependent variables are set up as:

```wl
In[34]:= tvars = {{u[t, x, y], v[t, x, y], p[t, x, y]}, t, {x, y}};
```

The transient Navier–Stokes operator is set for a Reynolds number of $1000$ :

```wl
In[35]:=
tpars = <|"ReynoldsNumber" -> 1000|>;
op = FluidFlowPDEComponent[tvars, tpars]

Out[36]= {{u[t, x, y], v[t, x, y]}.Inactive[Grad][u[t, x, y], {x, y}] + Inactive[Div][(-(1/500))*{{0, 1/4}, {1/4, 0}} . Inactive[Grad][v[t, x, y], {x, y}], {x, y}] + Inactive[Div][(-(1/500))*{{1, 0}, {0, 1/2}} . Inactive[Grad][u[t, x, y], {x, y}], {x, y}] + ... . Inactive[Grad][u[t, x, y], {x, y}], {x, y}] + Inactive[Div][(-(1/500))*{{1/2, 0}, {0, 1}} . Inactive[Grad][v[t, x, y], {x, y}], {x, y}] + Inactive[Div][{0, p[t, x, y]}, {x, y}] + v^(1, 0, 0)[t, x, y], v^(0, 0, 1)[t, x, y] + u^(0, 1, 0)[t, x, y]}
```

For the inflow boundary condition to be consistent with initial conditions, a helper function that smoothly ramps up the inflow over time is created.

A function to ramp up the inflow boundary condition over time:

```wl
In[37]:=
rampFunction[min_, max_, c_, r_] := Function[t, (min * Exp[c * r] + max * Exp[r * t]) / (Exp[c * r] + Exp[r * t])]
sf = rampFunction[0, 1, 4, 5];
Plot[sf[t], {t, -1, 10}, PlotRange -> All]

Out[39]= [image]
```

Set up the top boundary condition:

```wl
In[40]:= ttop = DirichletCondition[{u[t, x, y] == 1 * sf[t], v[t, x, y] == 0}, y == 1]

Out[40]= DirichletCondition[{u[t, x, y] == (E^5 t/E^20 + E^5 t), v[t, x, y] == 0}, y == 1]
```

Set up wall boundary conditions:

```wl
In[41]:= twall = DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, y < 1]

Out[41]= DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, y < 1]
```

Pressure point boundary condition:

```wl
In[42]:= treferencePressure = DirichletCondition[p[t, x, y] == 0, x == 0 && y == 0]

Out[42]= DirichletCondition[p[t, x, y] == 0, x == 0 && y == 0]
```

Combine the boundary conditions:

```wl
In[43]:= bcs = {ttop, twall, treferencePressure};
```

Initially the fluid is at rest. The initial conditions are set to 0:

```wl
In[44]:= ic = {u[0, x, y] == 0, v[0, x, y] == 0, p[0, x, y] == 0};
```

It is possible to monitor the time integration.

Time-integrate the Navier–Stokes equations:

```wl
In[50]:=
tEnd = 20;
AbsoluteTiming[Monitor[{uVel, vVel, pressure} = NDSolveValue[{op == {0, 0, 0}, bcs, ic}, {u, v, p}, {x, y}∈mesh, {t, 0, tEnd}, 
	Method -> {
	"PDEDiscretization" -> {"MethodOfLines", 
	"SpatialDiscretization" -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}
	}}}, EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])];, monitor]]

Out[51]= {67.4989, Null}
```

Visualize the velocity field for various times:

```wl
In[52]:=
frames = StreamPlot[{uVel[#, x, y], vVel[#, x, y]}, {x, y}∈Ω, PlotLabel -> "Time: " <> ToString[#] <> " s", Frame -> False]& /@ Range[0, tEnd, 1];
ListAnimate[frames]

Out[53]= DynamicModule[«8»]
```

See [this note](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1211049090) about improving the visual quality of the animation.

Equations

When planning a fluid flow simulation, a first step is to classify what type of flow is to be expected. For this, characteristic numbers that capture various aspects of flow behavior are used. One characterization is done through the Reynolds number $ℛℯ$. The Reynolds number is the ratio of inertial force in the fluid to the fluid's viscosity.

```wl
ℛℯ = (inertial force/viscosity) = (ρ v L/μ)
```

Here $μ$ [$Quantity[1, "Pascals" "Seconds"]$] is the dynamic viscosity, $ρ$ [$Quantity[1, ("Kilograms"/"Meters"^3)]$] is the density, $v$ [$Quantity[1, ("Meters"/"Seconds")]$] is a typical flow velocity and $L$ [$Quantity[1, "Meters"]$] is a characteristic length. The Reynolds number gives an approximate indication of what type of flow can be expected. The Reynolds number will be introduced more formally in the section [Dimensionless Form of the Navier-Stokes Equations](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1651516342). They are a variant of the fluid flow equations based on the Reynolds number. For now the Reynolds number servers as guide to classify flow.

Flows with $ℛℯ≲1$ are in the creeping flow regime. For flows where the Reynolds number is larger than, say, $ℛℯ≳5000$, turbulence can appear. At what particular Reynolds number a flow transitions from laminar to turbulent depends on the specifics of that particular flow. For example, in tubes, a value of $ℛℯ≳2000$ can induce turbulence flow. Contrary to that, flow along a flat plate will become turbulent only for $ℛℯ > 10 ^ 5$.

A distinction is typically made between the following fluid flow regimes:

* Creep flow: $ℛℯ≲1$

* Laminar flow: $1≲ℛℯ≲2000$

* Turbulent flow: $ℛℯ≳2000$

Honey, compared to water, is a viscous fluid. The higher the viscosity of a fluid, the more likely it will show laminar flow behavior. Gases, on the other hand have very low viscosity, and the Reynolds number approaches infinity. In most applications, gases are therefore considered to be inviscid, which means that for simplicity, their viscosity is assumed to be $0$. Gases are also called inviscid fluids. In some applications, such an assumption is not possible, however.

Viscous flow is described by a set of vector-valued equations called the Navier–Stokes equations. A wide range of fluid flow phenomena can be described with the Navier–Stokes equations, even beyond laminar flow. The Navier–Stokes equations are the overarching equations in fluid dynamics. The Navier–Stokes equations are a set of partial differential equations and exists in various forms, depending on the flow regime one is interested in.

Inviscid flow can be described by the Euler equations, which are not treated here.

The Navier–Stokes equations are a combination of a momentum equation and the continuity equation.

The momentum equation for a compressible fluid flow is a vector-valued equation expressed with the viscous stress tensor $**τ**$ in [$Pa$]:

```wl
ρ((∂**v**/∂t) + **v**·∇**v**) + ∇·(-**τ** + p**I**) - **F**	 = 	0
```

where

* $ρ$ is the density in [$kg / m^3$]

* $**v**$ is the fluid flow velocity vector in [$m / s$]

* $t$ is time in [$s$]

* $**τ**$ is the viscous stress tensor in [$Pa$]

* $p$ is the fluid pressure in [$Pa$]

* $**F**$ is a body force density in [$N / m^3$]

* $**I**$ is the identity matrix

The viscous stress tensor $**τ**$ is given as:

```wl
**τ**	 = 	2μ**Overscript[ϵ, .]** - (2/3)μ(∇·**v**)**I**
```

where

* $μ$ is the dynamic viscosity in [$Pa s$]

* $Overscript[**ϵ**, .]$ is the strain-rate tensor in [$1 / s$]

$Overscript[**ϵ**, .]$ is given by:

```wl
**Overscript[ϵ, .]**	 = (1/2)(∇**v** + (∇**v**)^T)
```

Inserting the definition of the viscous stress tensor $**τ**$ into the momentum equation gives:

```wl
ρ((∂**v**/∂t) + **v**·∇**v**) + ∇·(Underscript[-μ(∇**v** + (∇**v**)^T) + (2/3)μ(∇·**v**)**I**, Underscript[︸,                      viscous stress tensor **τ**                     ]] + p**I**) - **F**	 = 	0
```

Some formulations of the Navier–Stokes momentum equation contain an additional term $ζ$ called second viscosity, or volume viscosity, such that

```wl
**τ**	 = 	2μ**Overscript[ϵ, .]** - (2/3)μ(∇·v)I + ζ(∇·v)I
```

but for most practical purposes, $ζ$ is negligible and $ζ = 0$.

The unit of the momentum equation is a force density in [$N / m^3$]. The conservation of momentum is essentially Newton's second law of motion for fluids. This can be seen when the equation is slightly rearranged:

```wl
Underscript[ρ, Underscript[︸, mass / vol]]Underscript[((∂**v**/∂t) + **v**·∇**v**), Underscript[︸, acceleration / vol]]	 = 	Underscript[∇·(**τ** - p**I**) + **F**, Underscript[︸,          force / vol           ]]
```

The momentum equation is almost always solved together with the mass continuity equation, which describes the conservation of mass:

```wl
(∂ρ/∂t) + ∇·(**v** ρ)	 = 	0
```

with density $ρ$ the dependent variable and the fluid flow velocity vector $**v**$.

The Navier–Stokes equations are only valid when the smallest geometrical feature of the domain is still much larger then the mean free path of the fluid's molecules.

The Navier–Stokes equations are suitable to describe compressible, weakly compressible and incompressible flow. Weakly compressible flow means that a pressure dependence of density can be neglected and a is density evaluated at a reference pressure is made use of.

For constant values of the mass density $ρ$, the mass continuity equation simplifies to a volumetric continuity equation $∇·**v**  = 0$, and with that, the viscous stress tensor simplifies to $τ = 2μ**Overscript[ϵ, .]**$. When $ρ$ is constant, it is possible to write $∇·**v** ρ = 0$ as $ρ∇·v  = 0$, since in that case the divergence $∇·ρ  = 0$.

The incompressible form of the Navier–Stokes equations is:

```wl
|                                                                              |     |   |
| ---------------------------------------------------------------------------- | --- | - |
| ρ((∂**v**/∂t) + **v**·∇**v**) + ∇·(-μ(∇**v** + (∇**v**)^T) + p**I**) - **F** |  =  | 0 |
| ∇·**v**                                                                      |  =  | 0 |
```

The incompressible Navier–Stokes equations written out in component form are given as:

```wl
|                                                                                                 |     |   |
| ----------------------------------------------------------------------------------------------- | --- | - |
| ρ((∂u/∂t) + u(∂u/∂x) + v(∂u/∂y) + …) + (∂/∂x)(-μ(∂u/∂x)) + (∂/∂y)(-μ(∂u/∂y)) + … + (∂p/∂x) - fx |  =  | 0 |
| ρ((∂v/∂t) + u(∂v/∂x) + v(∂v/∂y) + …) + (∂/∂x)(-μ(∂v/∂x)) + (∂/∂y)(-μ(∂v/∂y)) + … + (∂p/∂y) - fy |  =  | 0 |
| ρ((∂w/∂t) + u(∂w/∂x) + v(∂w/∂y) + …) + (∂/∂x)(-μ(∂w/∂x)) + (∂/∂y)(-μ(∂w/∂y)) + … + (∂p/∂z) - fz |  =  | 0 |
| (∂u/∂x) + (∂v/∂y) + (∂w/∂z)                                                                     |  =  | 0 |
```

For a constant dynamics viscosity $μ$, this can be further simplified to:

```wl
|                                                                                                  |     |   |
| ------------------------------------------------------------------------------------------------ | --- | - |
| ρ((∂u/∂t) + u(∂u/∂x) + v(∂u/∂y) + …) - μ(∂^2u/∂x^2) - μ(∂^2u/∂y^2) - μ(∂^2u/∂z^2) + (∂p/∂x) - fx |  =  | 0 |
| ρ((∂v/∂t) + u(∂v/∂x) + v(∂v/∂y) + …) - μ(∂^2v/∂x^2) - μ(∂^2v/∂y^2) - μ(∂^2v/∂z^2) + (∂p/∂y) - fy |  =  | 0 |
| ρ((∂w/∂t) + u(∂w/∂x) + v(∂w/∂y) + …) - μ(∂^2w/∂x^2) - μ(∂^2w/∂y^2) - μ(∂^2w/∂z^2) + (∂p/∂z) - fz |  =  | 0 |
| (∂u/∂x) + (∂v/∂y) + (∂w/∂z)                                                                      |  =  | 0 |
```

As an example, a 3D time-dependent flow through connected pipes is modeled.

Create a geometry:

```wl
In[6]:=
radius = 1 / 5;
Ω = RegionUnion[Ball[], Cylinder[{{0, 0, 0}, {0, 0, 2}}, radius], Cylinder[{{0, 0, 0}, {0, 2, 0}}, radius]];
```

Create a somewhat refined mesh from the geometry:

```wl
In[8]:= mesh = ToElementMesh[Ω, MaxCellMeasure -> 0.001]

Out[8]= ElementMesh[{{-0.999996, 0.994936}, {-0.999996, 2.}, {-1., 2.}}, {TetrahedronElement[<10965>]}]
```

Visualize the mesh:

```wl
In[9]:= mesh["Wireframe"]

Out[9]= [image]
```

Specify the variables and material parameters:

```wl
In[11]:=
vars = {{u[t, x, y, z], v[t, x, y, z], w[t, x, y, z], p[t, x, y, z]}, t, {x, y, z}};
pars = <|"DynamicViscosity" -> 10 ^ -3, "MassDensity" -> 1|>;
```

Set up the equation:

```wl
In[13]:= pde = FluidFlowPDEComponent[vars, pars] == {0, 0, 0, 0};
```

Create a function to slowly ramp up the inflow velocity:

```wl
In[8]:= ramp = Function[t, Exp[5 * t] / (Exp[20] + Exp[5 * t])];
```

Create and visualize a parabolic inflow profile:

```wl
In[9]:=
flowProfile[x_, y_] := -(x ^ 2 / radius ^ 2 + y ^ 2 / radius ^ 2) + 1
Plot3D[flowProfile[x, y], {x, y}∈Disk[{0, 0}, radius]]

Out[10]= [image]
```

Create the inflow boundary condition:

```wl
In[14]:=
velocity = 3 / 4;
inflowBC = DirichletCondition[{u[t, x, y, z] == 0, v[t, x, y, z] == 0, w[t, x, y, z] == -ramp[t] * velocity * flowProfile[x, y]}, z == 2];
```

Set up the pressure at the outflow:

```wl
In[16]:= outflowBC = DirichletCondition[p[t, x, y, z] == 0, y == 2];
```

Specify a no-slip wall boundary condition:

```wl
In[17]:= wallBC = DirichletCondition[{u[t, x, y, z] == 0, v[t, x, y, z] == 0, w[t, x, y, z] == 0}, z != 2 && y != 2];
```

Collect the boundary conditions:

```wl
In[18]:= bcs = {inflowBC, outflowBC, wallBC};
```

Set up initial conditions at rest:

```wl
In[19]:= ics = {u[0, x, y, z] == 0, v[0, x, y, z] == 0, w[0, x, y, z] == 0, p[0, x, y, z] == 0};
```

Time-integrate the equation and monitor its progress:

```wl
In[20]:= Monitor[AbsoluteTiming[{uVel, vVel, zVel, pressure} = NDSolveValue[{pde, bcs, ics}, {u, v, w, p}, {x, y, z}∈mesh, {t, 0, 10}, Method -> {"PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, w -> 2, p -> 1}}}}, EvaluationMonitor :> (currentTime = Row[{"t = ", CForm[t]}])];], currentTime]

Out[20]= {566.615, Null}
```

Find the minimal and maximal value of the flow velocities:

```wl
In[28]:= minmaxVel = MinMax[#["ValuesOnGrid"]& /@ {uVel, vVel, zVel}]

Out[28]= {-0.791355, 0.683604}
```

Visualize the flow vectors and slices of a contour plot through the geometry:

```wl
In[27]:= Show[mesh["Edgeframe"], SliceContourPlot3D[Evaluate[Sqrt[uVel[10, x, y, z] ^ 2 + vVel[10, x, y, z] ^ 2 + zVel[10, x, y, z] ^ 2]], {x, y, z}∈mesh, ...], VectorPlot3D[{uVel[10, x, y, z], vVel[10, x, y, z], zVel[10, x, y, z]}, {x, y, z}∈mesh, ...], Rule[...]]

Out[27]= [image]
```

Dimensionless Form of the Navier-Stokes Equations

Sometimes it is of interest to replicate fluid flow experiments that have a large spacial extent, for example the flow around a container ship, with a spatially scaled-down prototype. The question then arises as how to use the results from the small scale experiment to make statements about the actual, large scale fluid flow setup. One approach to this issue is to transform the Navier–Stokes equations into a dimensionless form. In this form, the equations are independent of units and make meaningful predictions in the large scale set up feasible.

In order to derive the dimensionless form, it is necessary to define a suitable set of dimensionless variables. These variables are based on an appropriate but arbitrary scale.

```wl
x = Lx^ * 
t = (L/V)t^ * 
Subscript[v, i] = VSubscript[v, i]^ * 
p = ρV^2p^ *
```

where $L$ [$Quantity[1, "Meters"]$] is the characteristic length and $V$ [$Quantity[1, ("Meters"/"Seconds")]$] is the characteristic flow speed. The asterisk $^ * $ represents a dimensionless quantity. The characteristic length and speed are related to the problem one is interested in. For example, the characteristic length may be the size of the domain and the characteristic flow speed may be the speed at which fluid flows into the domain. Note that even though the length, time and velocity have a natural scale, this is not the case for pressure.

The specific choice of scaling for pressure presented here is rather arbitrary and it possible to consider a different pressure scale that would then result in slightly different system of equations.

Plugging these identities into the Navier–Stokes equations leads to the dimensionless form of the Navier–Stokes equations:

```wl
|     |     |     |
| --- | --- | --- |
| ρ((∂**v**^ * /∂t^ * ) + **v^ * **·∇^ * **v^ * **) + ∇^ * ·(-(μ/ρVL)(∇^ * **v^ * ** + (∇^ * **v^ * **)^T) + p^ * **I**) |  =  | **0** |
| ∇^ * ·**v^ * ** |  =  | 0   |
```

The symbol ∇^ *  represents gradient with respect to the dimensionless coordinates $x^ * $. Also note that the external force $**F**$ has been set to $0$. In case it is non-zero a $**F**^ * $ needs to be introduced and a new dimensionless number, the Froude number will be needed.

The dimensionless incompressible Navier–Stokes equations written out in component form are given as:

```wl
|     |     |   |
| --- | --- | - |
| ρ((∂u^ * /∂t^ * ) + u^ * (∂u^ * /∂x^ * ) + v^ * (∂u^ * /∂y^ * ) + …) - (μ/ρVL)((∂^2u^ * /∂x^ * ^2) + (∂^2u^ * /∂y^ * ^2) + (∂^2u^ * /∂z^ * ^2)) + (∂p^ * /∂x^ * ) |  =  | 0 |
| ρ((∂v^ * /∂t^ * ) + u^ * (∂v^ * /∂x^ * ) + v^ * (∂v^ * /∂y^ * ) + …) - (μ/ρVL)((∂^2v^ * /∂x^ * ^2) + (∂^2v^ * /∂y^ * ^2) + (∂^2v^ * /∂z^ * ^2)) + (∂p^ * /∂y^ * ) |  =  | 0 |
| ρ((∂w^ * /∂t^ * ) + u^ * (∂w^ * /∂x^ * ) + v^ * (∂w^ * /∂y^ * ) + …) - (μ/ρVL)((∂^2w^ * /∂x^ * ^2) + (∂^2w^ * /∂y^ * ^2) + (∂^2w^ * /∂z^ * ^2)) + (∂p^ * /∂z^ * ) |  =  | 0 |
| (∂u^ * /∂x^ * ) + (∂v^ * /∂y^ * ) + (∂w^ * /∂z^ * ) |  =  | 0 |
```

The inverse of the fraction μ/(ρ V L) present in the dimensionless Navier–Stokes equations is called the Reynolds number:

```wl
ℛℯ = (ρVL/μ) = (inertial force/viscous force)
```

Since the Reynolds number is dimensionless, two different flows can have the same behavior only if their Reynolds numbers are the same. However, the Reynolds number is not the only dimensionless number that can occur in fluid dynamics; therefore, the same Reynolds numbers might not be sufficient for flow comparison.

Cylindrical Coordinates

When modeling a fluid flow problem, sometimes it is not convenient to describe the model in Cartesian coordinates $(x, y, z)$. In such cases, the equations may also be expressed using a cylindrical coordinate system.

[image]

A graphic showing cylindrical coordinates in relation to Cartesian coordinates.

The Navier–Stokes equations can be expressed in cylindrical coordinates. There are a few simplifications, such as axisymmetric fluid flow, that are derived from the Navier–Stokes equations in cylindrical coordinates, and for this reason, the equations are presented next. For simplicity consider, but without restriction to generality, the incompressible Navier–Stokes equations [(9)](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#528541537). They have the following form in the cylindrical coordinates:

ρ((∂Subscript[v, r]/∂t) + Subscript[v, r](∂Subscript[v, r]/∂r) + (Subscript[v, θ]/r)(∂Subscript[v, r]/∂θ) + Subscript[v, z](∂Subscript[v, r]/∂z) - (Subscript[v, θ]^2/r)) + (∂p/∂r) - μ((1/r)(∂/∂r)(r(∂Subscript[v, r]/∂r)) + (1/r^2)(∂^2Subscript[v, r]/∂^2θ) + (∂^2Subscript[v, r]/∂^2z) - (Subscript[v, r]/r^2) - (2/r^2)(∂Subscript[v, θ]/∂θ)) - ρSubscript[g, r] = 0
ρ((∂Subscript[v, θ]/∂t) + Subscript[v, r](∂Subscript[v, θ]/∂r) + (Subscript[v, θ]/r)(∂Subscript[v, θ]/∂θ) + Subscript[v, z](∂Subscript[v, θ]/∂z) + (Subscript[v, r]Subscript[v, θ]/r)) + (1/r)(∂p/∂θ) - μ((1/r)(∂/∂r)(r(∂Subscript[v, θ]/∂r)) + (1/r^2)(∂^2Subscript[v, θ]/∂^2θ) + (∂^2Subscript[v, θ]/∂^2z) - (Subscript[v, θ]/r^2) - (2/r^2)(∂Subscript[v, r]/∂θ)) - ρSubscript[g, θ] = 0
ρ((∂Subscript[v, z]/∂t) + Subscript[v, r](∂Subscript[v, z]/∂r) + (Subscript[v, θ]/r)(∂Subscript[v, z]/∂θ) + Subscript[v, z](∂Subscript[v, z]/∂z)) + (∂p/∂z) - 
		       μ((1/r)(∂/∂r)(r(∂Subscript[v, z]/∂r)) + (1/r^2)(∂^2Subscript[v, z]/∂^2θ) + (∂^2Subscript[v, z]/∂^2z)) - ρSubscript[g, z] = 0
(1/r)(∂v/∂r)(rSubscript[v, r]) + (1/r)(∂Subscript[v, θ]/∂θ) + (∂Subscript[v, z]/∂z) = 0

where $r$ denotes the radial component, $θ$ is the angular component, and $z$ is the same component as in the Cartesian coordinates.

In terms of the Cartesian coordinates $(x, y, z)$, the cylindrical coordinates $(r, θ, z)$ are defined by:

```wl
{|   |     |                 |
| - | --- | --------------- |
| r |  =  | Sqrt[x^2 + y^2] |
| θ |  =  | arctan(x, y)    |
| z |  =  | z               |
```

Convert from Cartesian to cylindrical coordinates:

```wl
In[1]:= CoordinateTransform[ "Cartesian" -> "Cylindrical", {x, y, z}]

Out[1]= {Sqrt[x^2 + y^2], ArcTan[x, y], z}
```

or

```wl
{|   |     |          |
| - | --- | -------- |
| x |  =  | r·cos(θ) |
| y |  =  | r·sin(θ) |
| z |  =  | z        |
```

Convert from cylindrical to Cartesian coordinates:

```wl
In[2]:= CoordinateTransform[ "Cylindrical" -> "Cartesian", {r, θ, z}]

Out[2]= {r Cos[θ], r Sin[θ], z}
```

These are the axisymmetric incompressible Navier–Stokes equations. They can be generated by ``FluidFlowPDEComponent`` by setting ``"CoordinateChart"`` to ``"Cylindrical"``.

Create a symbolic stationary axisymmetric fluid dynamics PDE with a dynamic viscosity of $μ$ and a mass density of $ρ$ :

```wl
In[9]:= FluidFlowPDEComponent[{{u[r, θ, z], v[r, θ, z], w[r, θ, z], p[r, θ, z]}, {r, θ, z}}, <|"DynamicViscosity" -> μ, "MassDensity" -> ρ, "CoordinateChart" -> "Cylindrical"|>]//Short

Out[9]//Short= {«1», «1», «1», ρ w^(0, 0, 1)[r, θ, z] + (ρ u[r, θ, z] + ρ «1»[«1»]/r) + ρ u^(1, 0, 0)[r, θ, z]}
```

The time-dependent equations can be created in a similar manner.

Create a symbolic time-dependent axisymmetric fluid dynamics PDE with a dynamic viscosity of $μ$ and a mass density of $ρ$ :

```wl
In[8]:= FluidFlowPDEComponent[{{u[t, r, θ, z], v[t, r, θ, z], w[t, r, θ, z], p[t, r, θ, z]}, t, {r, θ, z}}, <|"DynamicViscosity" -> μ, "MassDensity" -> ρ, "CoordinateChart" -> "Cylindrical"|>]//Short

Out[8]//Short= {«10» + ρ u^(1, «3»)[t, r, θ, z], «1», «1», ρ w^(0, 0, 0, 1)[t, r, θ, z] + (ρ «1» + «1»/r) + ρ u^(0, 1, 0, 0)[t, r, θ, z]}
```

Axisymmetric models

Equations arising from fluid dynamics are often highly nonlinear, which makes them difficult and expensive to solve. Therefore, any possible simplification of the equations is highly attractive, as it allows the equations to be solved with fewer resources. A common approach is to assume that the domain has some kind of symmetry. Pipes and similar geometries usually have an axis of revolution around which they are symmetric. This leads to a notion of axisymmetric flows that are assumed to be independent of the angular component. This allows the number of dimensions to be reduced from $3$ to $2$, which makes the computation less expensive both in time and memory.

[image]

The illustration shows a body that has an axis of revolution around the dashed $z$ axis. The body can be represented in an axisymmetric setting by making use of a 2D area that is rotated around the $z$ axis, following the angular direction $θ$.

The derivation of the axisymmetric equations starts from the Navier–Stokes equations in cylindrical coordinates. The axisymmetric equations can be derived from this system of equations by assuming that $∂ / ∂θ = 0$, which means no gradients in the $θ$ direction. Additionally $Subscript[v, ϕ] = 0$ is used. Also, $Subscript[v, r]$, $Subscript[v, z]$ and $p$ do not depend on $ϕ$, which leads to the following truncated set of equations:

ρ((∂Subscript[v, r]/∂t) + Subscript[v, r](∂Subscript[v, r]/∂r) + Subscript[v, z](∂Subscript[v, r]/∂z)) + (∂p/∂r) - μ((1/r)(∂/∂r)(r(∂Subscript[v, r]/∂r)) + (∂^2Subscript[v, r]/∂^2z) - (Subscript[v, r]/r^2)) - ρSubscript[g, r] = 0
ρ((∂Subscript[v, z]/∂t) + Subscript[v, r](∂Subscript[v, z]/∂r) + Subscript[v, z](∂Subscript[v, z]/∂z)) + (∂p/∂z) - μ((1/r)(∂/∂r)(r(∂Subscript[v, z]/∂r)) + (∂^2Subscript[v, z]/∂^2z)) - ρSubscript[g, z] = 0
(1/r)(∂v/∂r)(rSubscript[v, r]) + (∂Subscript[v, z]/∂z) = 0

These are the axisymmetric incompressible Navier–Stokes equations. They can be generated by ``FluidFlowPDEComponent`` by setting ``"RegionSymmetry"`` to ``"Axisymmetric"``.

Create a symbolic stationary axisymmetric fluid dynamics PDE with a dynamic viscosity of $μ$ and a mass density of $ρ$ :

```wl
In[4]:= FluidFlowPDEComponent[{{u[r, z], 0, w[r, z], p[r, z]}, {r, θ, z}}, <|"DynamicViscosity" -> μ, "MassDensity" -> ρ, "RegionSymmetry" -> "Axisymmetric"|>]

Out[4]= {{ρ u[r, z], ρ w[r, z]}.Inactive[Grad][u[r, z], {r, z}] + p^(1, 0)[r, z] + ((2 μ u[r, z]/r) - 2 μ u^(1, 0)[r, z]/r) - μ (u^(0, 2)[r, z] + w^(1, 1)[r, z]) - 2 μ u^(2, 0)[r, z], {ρ u[r, z], ρ w[r, z]}.Inactive[Grad][w[r, z], {r, z}] + p^(0, 1)[r, z] - 2 μ w^(0, 2)[r, z] - (μ (u^(0, 1)[r, z] + w^(1, 0)[r, z])/r) - μ (u^(1, 1)[r, z] + w^(2, 0)[r, z]), (ρ u[r, z]/r) + ρ w^(0, 1)[r, z] + ρ u^(1, 0)[r, z]}
```

The time-dependent equations can be created in a similar manner.

Create a symbolic time-dependent axisymmetric fluid dynamics PDE with a dynamic viscosity of $μ$ and a mass density of $ρ$ :

```wl
In[5]:= FluidFlowPDEComponent[{{u[t, r, z], 0, w[t, r, z], p[t, r, z]}, t, {r, θ, z}}, <|"DynamicViscosity" -> μ, "MassDensity" -> ρ, "RegionSymmetry" -> "Axisymmetric"|>]

Out[5]= {{ρ u[t, r, z], ρ w[t, r, z]}.Inactive[Grad][u[t, r, z], {r, z}] + p^(0, 1, 0)[t, r, z] + ((2 μ u[t, r, z]/r) - 2 μ u^(0, 1, 0)[t, r, z]/r) - μ (u^(0, 0, 2)[t, r, z] + w^(0, 1, 1)[t, r, z]) - 2 μ u^(0, 2, 0)[t, r, z] + ρ u^(1, 0, 0)[t, r, z], {ρ u[ ... + p^(0, 0, 1)[t, r, z] - 2 μ w^(0, 0, 2)[t, r, z] - (μ (u^(0, 0, 1)[t, r, z] + w^(0, 1, 0)[t, r, z])/r) - μ (u^(0, 1, 1)[t, r, z] + w^(0, 2, 0)[t, r, z]) + ρ w^(1, 0, 0)[t, r, z], (ρ u[t, r, z]/r) + ρ w^(0, 0, 1)[t, r, z] + ρ u^(0, 1, 0)[t, r, z]}
```

The following example shows the setup of an axisymmetric fluid flow model. It considers a pipe composed of three different sections with different lengths and radii, a parabolic inflow at the bottom and the outlet at the top.

Specify the geometry:

```wl
In[10]:=
geometryData = <|l1 -> 3, l2 -> 2, l3 -> 4, r1 -> 1.5, r2 -> 1, r3 -> 2|>;
Ω = RegionUnion[Rectangle[{0, 0}, {r1, l1}], Rectangle[{0, l1}, {r2, l1 + l2}], Rectangle[{0, l1 + l2}, {r3, l1 + l2 + l3}]] /. geometryData;
```

Create and visualize the mesh:

```wl
In[12]:=
mesh = ToElementMesh[Ω];
mesh["Wireframe"]

Out[13]= [image]
```

Specify the variables:

```wl
In[10]:= vars = {{u[r, z], 0, w[r, z], p[r, z]}, {r, θ, z}};
```

Note how the dependent variable at the position of the angular direction is set to $0$.

Specify model parameters:

```wl
In[9]:= pars = <|"DynamicViscosity" -> 0.01, "MassDensity" -> 1, "RegionSymmetry" -> "Axisymmetric"|>;
```

Set up the inflow boundary condition:

```wl
In[12]:= inflow = DirichletCondition[{u[r, z] == 0, w[r, z] == 0.1(r1 ^ 2 - r ^ 2)}, z == 0] /. geometryData;
```

Set up the wall boundary condition:

```wl
In[13]:= wall = DirichletCondition[{u[r, z] == 0, w[r, z] == 0}, r > 0 && z != 0 && z != l1 + l2 + l3] /. geometryData;
```

Set up an outflow boundary condition:

```wl
In[14]:= outflow = DirichletCondition[u[r, z] == 0, z == l1 + l2 + l3] /. geometryData;
```

Set up the axisymmetric symmetry boundary condition:

```wl
In[15]:= symmetry = DirichletCondition[u[r, z] == 0, r == 0];
```

Set up the pressure boundary condition to make the solution unique:

```wl
In[16]:= pressure = DirichletCondition[p[r, z] == 0, r == 0 && z == l1 + l2 + l3] /. geometryData;
```

Combine the boundary conditions

```wl
In[17]:= bcs = {inflow, wall, outflow, symmetry, pressure};
```

Create the PDE:

```wl
In[18]:= pde = FluidFlowPDEComponent[vars, pars] == {0, 0, 0};
```

Solve the equations and measure the time it takes:

```wl
In[19]:=
AbsoluteTiming[{rVel, zVel, pressure} = NDSolveValue[{pde, 
	bcs}, {u[r, z], w[r, z], p[r, z]}, {r, z}∈mesh, 
	Method -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, w -> 2, p -> 1}}];]

Out[19]= {1.03219, Null}
```

Use ``VectorPlot`` to plot the solution:

```wl
In[20]:= VectorPlot[{rVel, zVel}, {r, z}∈mesh]

Out[20]= [image]
```

Use ``StreamPlot`` to plot the solution:

```wl
In[21]:= StreamPlot[{rVel, zVel}, {r, z}∈mesh]

Out[21]= [image]
```

Use`` ContourPlot`` to plot the pressure distribution:

```wl
In[23]:= ContourPlot[pressure, {r, z}∈mesh]

Out[23]= [image]
```

Sometimes it might be useful to visualize the flow in 3D.

Create the domain in 3D:

```wl
In[33]:= domain3D = RegionUnion[Cylinder[{{0, 0, 0}, {0, 0, l1}}, r1], Cylinder[{{0, 0, l1}, {0, 0, l1 + l2}}, r2], Cylinder[{{0, 0, l1 + l2}, {0, 0, l1 + l2 + l3}}, r3]] /. geometryData

Out[33]= [image]
```

Visualize the flow in 3D:

```wl
In[38]:=
Show[ToElementMesh[domain3D]["Edgeframe"], 
	StreamPlot3D[Evaluate[{Cos[ArcTan[x, y]] * rVel /. {r -> Sqrt[x ^ 2 + y ^ 2]}, Sin[ArcTan[x, y]] * rVel /. {r -> Sqrt[x ^ 2 + y ^ 2]}, zVel /. {r -> Sqrt[x ^ 2 + y ^ 2]}}], {x, y, z}∈domain3D, Axes -> None, Boxed -> False]]

Out[38]= [image]
```

The 3D visualization is done in two steps. First the $x$, $y$ and $z$ coordinates, which the plot uses, are transformed to the cylindrical coordinates, which are then plugged into the solution. The results are then transformed back into the Cartesian coordinates.

The time-dependent problem can be solved similarly.

Specify the variables:

```wl
In[10]:= tvars = {{u[t, r, z], 0, w[t, r, z], p[t, r, z]}, t, {r, θ, z}};
```

Define and plot the ramp function:

```wl
In[11]:=
rampFunction[min_, max_, c_, r_] := Function[t, (min * Exp[c * r] + max * Exp[r * t]) / (Exp[c * r] + Exp[r * t])]
ramp = rampFunction[0, 1, 4, 5];
Plot[ramp[t], {t, -1, 10}, PlotRange -> All]

Out[13]= [image]
```

Set up the inflow boundary condition:

```wl
In[14]:= tinflow = DirichletCondition[{u[t, r, z] == 0, w[t, r, z] == ramp[t] * 0.1 * (r1 ^ 2 - r ^ 2)}, z == 0] /. geometryData;
```

Set up the wall boundary condition:

```wl
In[15]:= twall = DirichletCondition[{u[t, r, z] == 0, w[t, r, z] == 0}, r > 0 && z != 0 && z != l1 + l2 + l3] /. geometryData;
```

Set up an outflow boundary condition:

```wl
In[16]:= toutflow = DirichletCondition[u[t, r, z] == 0, z == l1 + l2 + l3] /. geometryData;
```

Set up the axisymmetric symmetry boundary condition:

```wl
In[17]:= tsymmetry = DirichletCondition[u[t, r, z] == 0, r == 0];
```

Set up the pressure boundary condition to make the solution unique:

```wl
In[18]:= tpressure = DirichletCondition[p[t, r, z] == 0, r == 0 && z == l1 + l2 + l3] /. geometryData;
```

Combine the boundary conditions:

```wl
In[19]:= tbcs = {tinflow, twall, toutflow, tsymmetry, tpressure};
```

Set up the initial condition:

```wl
In[20]:= ic = {u[0, r, z] == 0, w[0, r, z] == 0, p[0, r, z] == 0};
```

Create the PDE:

```wl
In[21]:= pde = FluidFlowPDEComponent[tvars, pars] == {0, 0, 0};
```

Solve the equations and measure the time it takes:

```wl
In[22]:=
tEnd = 10;
AbsoluteTiming[{rVel, zVel, pressure} = NDSolveValue[{pde, 
	tbcs, ic}, {u, w, p}, {r, z}∈mesh, {t, 0, tEnd}, 
	Method -> {"PDEDiscretization" -> {"MethodOfLines", 
	"SpatialDiscretization" -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, w -> 2, p -> 1}}}}];]

Out[23]= {4.62708, Null}
```

Use ``StreamPlot`` to plot the solution:

```wl
In[24]:=
frames = StreamPlot[{rVel[#, r, z], zVel[#, r, z]}, {r, z}∈mesh, PlotLabel -> "Time: " <> ToString[#] <> " s", Frame -> False]& /@ Range[1, tEnd, 1];
frames = Rasterize[#, "Image", ImageResolution -> 120]& /@ frames;
ListAnimate[frames]

Out[26]= DynamicModule[«8»]
```

Use ``ContourPlot`` to plot the pressure distribution:

```wl
In[30]:=
frames = ContourPlot[pressure[#, r, z], {r, z}∈mesh, PlotLabel -> "Time: " <> ToString[#] <> " s", Frame -> False]& /@ Range[1, tEnd, 1];
frames = Rasterize[#, "Image", ImageResolution -> 120]& /@ frames;
ListAnimate[frames]

Out[32]= DynamicModule[«8»]
```

Visualize the flow in 3D:

```wl
In[34]:=
frames = Show[StreamPlot3D[{Cos[ArcTan[x, y]] * rVel[#, Sqrt[x ^ 2 + y ^ 2], z], Sin[ArcTan[x, y]] * rVel[#, Sqrt[x ^ 2 + y ^ 2], z], zVel[#, Sqrt[x ^ 2 + y ^ 2], z]}, {x, y, z}∈domain3D, ...], ToElementMesh[domain3D]["Edgeframe"]]& /@ Range[1, tEnd, 1];
frames = Rasterize[#, "Image", ImageResolution -> 120]& /@ frames;
ListAnimate[frames]

Out[36]= DynamicModule[«8»]
```

Swirling flow

The axisymmetric case was derived based on the assumption that $∂ / ∂θ = 0$ and by setting $Subscript[v, ϕ] = 0$. This allowed the $θ$ equation to be removed. When only the $∂ / ∂θ = 0$ assumption is used, it is possible to model swirl flow. This means that the $θ$ equation is maintained and allows an angular velocity to be specified. A rotating mixer in a beaker can be modeled with such a setup.

Specify the geometry of a beaker and a mixer:

```wl
In[14]:=
bw = 0.02;bh = 0.05;
radius = 0.001;height = 0.015;
beaker = Rectangle[{0, 0}, {bw, bh}];
shaft = Rectangle[{0, height}, {radius, bh}];
rotorThickness = 2radius;
rotor1 = Rectangle[{0, height}, {1 / 2 * bw, height + rotorThickness}];
rotor2 = Rectangle[{0, height + 2 * rotorThickness}, {1 / 3 * bw, height + 3rotorThickness}];
rotor3 = Rectangle[{0, height + 4 * rotorThickness}, {1 / 2 * bw, height + 5rotorThickness}];
mixer = RegionUnion[shaft, rotor1, rotor2, rotor3];
Ω = RegionDifference[beaker, mixer]

Out[23]=
Polygon[{{0.02, 0.05}, {0.02, 0.}, {0., 0.}, {0., 0.015}, {0.001, 0.015}, {0.01, 0.015}, 
  {0.01, 0.017}, {0.001, 0.017}, {0.001, 0.019}, {0.006666666666666666, 0.019}, 
  {0.006666666666666666, 0.020999999999999998}, {0.001, 0.020999999999999998}, {0.001, 0.023}, 
  {0.01, 0.023}, {0.01, 0.025}, {0.001, 0.025}, {0.001, 0.05}}]
```

Create a mesh and visualize the refined mesh:

```wl
In[24]:=
mesh = ToElementMesh[Ω, "MaxCellMeasure" -> 0.0000005];
mesh["Wireframe"]

Out[25]= [image]
```

In this case, the velocity in the $θ$ direction is not zero, and thus the dependent variable is specified in the set of dependent variables.

Specify the variables:

```wl
In[16]:= vars = {{u[r, z], v[r, z], w[r, z], p[r, z]}, {r, θ, z}};
```

Set the material properties and the coordinate chart.

```wl
In[17]:= pars = <|"DynamicViscosity" -> 10 ^ -3, "MassDensity" -> 10 ^ 3, "CoordinateChart" -> "Cylindrical"|>;
```

Create the operator:

```wl
In[18]:= op = FluidFlowPDEComponent[vars, pars];
```

In this version of the Wolfram Language, the finite element parser can not process terms like ``Inactive[Grad][u[r, z], {r, θ, z}]``, but the remedy is easy: The equation can be activated, and then all terms can be parsed.

Activate the equation:

```wl
In[19]:= op = Activate[op];
```

Set up the parametric boundary conditions for the mixer with the angular velocity $ω$ as a parameter:

```wl
In[29]:= mixerBC = DirichletCondition[{u[r, z] == 0, v[r, z] == r * ω, w[r, z] == 0}, height <= z < bh && radius <= r < bw];
```

Set up no-slip boundary conditions at the beaker wall:

```wl
In[30]:= wall = DirichletCondition[{u[r, z] == 0, v[r, z] == 0, w[r, z] == 0}, r == bw || z == 0];
```

On the axis of symmetry, the $u$ and $v$ velocity are set to 0:

```wl
In[31]:= symmetryAxis = DirichletCondition[{u[r, z] == 0, v[r, z] == 0}, r == 0];
```

The top of the beaker is open, and there is no fluid velocity in the $z$ direction:

```wl
In[32]:= top = DirichletCondition[{w[r, z] == 0}, z == bh];
```

A reference pressure condition is set at the top-right corner:

```wl
In[33]:= refPressure = DirichletCondition[{p[r, z] == 0}, r == bw && z == bh];
```

The parametric function is created:

```wl
In[34]:=
fun = ParametricNDSolveValue[{op == {0, 0, 0, 0}, mixerBC, wall, symmetryAxis, top, refPressure}, {u, v, w, p}, {r, z}∈mesh, ω, Method -> {"PDEDiscretization" -> {"FiniteElement"
	, "InterpolationOrder" -> {u -> 2, v -> 2, w -> 2, p -> 1}
	}}]

Out[34]= ParametricFunction[<>]
```

Compute the solution for an angular velocity of $π$ [$Quantity[1, ("Meters"/"Seconds")]$]:

```wl
In[35]:= {ufun, vfun, wfun, pfun} = fun[π]

Out[35]=
{InterpolatingFunction[{{0., 0.02}, {0., 0.05}}, {5, 4225, 0, {6013, 0}, {3, 0}, 0, 0, 0, 0, 
  Indeterminate & , {}, {}, False}, 
 {NDSolve`FEM`ElementMesh[CompressedData["«119831»"], {NDSolve`FEM`TriangleElement[CompressedData["«44957»"]]}, 
   { ... , 13, 13, 13, 14, 11, 11, 11, 13, 13, 13, 13, 11, 13, 13, 14, 14, 14, 14, 14, 
      14, 14, 14, 14, 14, 14, 14, 14, 14, 11, 11, 11, 14, 14, 11, 11, 11, 11, 11, 11, 14, 11, 11, 
      5, 11, 11, 9, 11}]}]}, CompressedData["«16763»"], {Automatic}]}
```

Compute the magnitude of the velocity:

```wl
In[36]:= velMag = EvaluateOnElementMesh[{r, z}, Sqrt[ufun[r, z] ^ 2 + vfun[r, z] ^ 2 + wfun[r, z] ^ 2], mesh]

Out[36]=
InterpolatingFunction[{{0., 0.02}, {0., 0.05}}, {5, 4225, 0, {6013, 0}, {3, 0}, 0, 0, 0, 0, 
  Automatic, {}, {}, False}, {NDSolve`FEM`ElementMesh[CompressedData["«119829»"], {NDSolve`FEM`TriangleElement[CompressedData["«44957»"]]}, {NDSolve`FEM`Li ... 4, 14, 
      14, 14, 14, 14, 14, 14, 11, 11, 11, 14, 11, 11, 14, 11, 11, 11, 11, 11, 14, 11, 11, 5, 11, 
      11, 9, 11}]}, {NDSolve`FEM`PointElement[CompressedData["«975»"], CompressedData["«311»"]]}]}, 
 CompressedData["«63821»"], {Automatic}]
```

Visualize the magnitude of the velocity and contours of the $v$ component of the velocity:

```wl
In[37]:=
Show[
	HighlightMesh[mixer, {Style[1, Red], Style[2, LightGray]}], 
	DensityPlot[velMag[r, z], {r, z}∈mesh, ...], ...]

Out[37]= [image]
```

Extract the minimal and maximal values of the velocity magnitude:

```wl
In[38]:= {minMag, maxMag} = MinMax[velMag["ValuesOnGrid"]]

Out[38]= {-4.058614618256235`*^-6, 0.0314267}
```

Create animation frames for increased angular velocity of the mixer:

```wl
In[39]:=
legendBar = BarLegend[{"TemperatureMap", {minMag, maxMag}}, ...];
frames = Table[{ufun, vfun, wfun, pfun} = fun[i * π];velMag = EvaluateOnElementMesh[{r, z}, Sqrt[ufun[r, z] ^ 2 + vfun[r, z] ^ 2 + wfun[r, z] ^ 2], mesh];Show[HighlightMesh[mixer, {Style[1, Red], Style[2, LightGray]}], Legended[DensityPlot[velMag[r, z], {r, z}∈mesh, ...], legendBar], ...], {i, 0.1, 1, 0.1}];
frames = Rasterize[#, "Image", ImageResolution -> 120]& /@ frames;
ListAnimate[frames]

Out[42]= DynamicModule[«8»]
```

Newtonian versus Non-Newtonian Flow

The viscosity of a non-Newtonian fluid changes with the applied forces. For a Newtonian fluid, such as water or honey, the viscosity does not change depending on how much the fluid is stirred. For non-Newtonian fluids, the amount of stirring changes the viscosity. Some materials become harder when they are stirred, while others become easier to stir. In practice, all gases and most liquids can be considered Newtonian. Non-Newtonian liquids are things like paint, blood, water-starch mixtures or liquid polymers.

In the following section, various non-Newtonian fluid flow models are introduced. For the derivation, consider the compressible momentum equation:

```wl
ρ((∂**v**/∂t) + **v**·∇**v**) + ∇·(-μ(∇**v** + (∇**v**)^T) + (2/3)μ(∇·**v**)**I** + p**I**)	 = 	**F**
```

and introduce a factor:

```wl
ρ((∂**v**/∂t) + **v**·∇**v**) + ∇·(-2μ(1/2)(∇**v** + (∇**v**)^T) + (2/3)μ(∇·**v**)**I** + p**I**)	 = 	**F**
```

Now it is possible to define the strain rate tensor $**Overscript[ϵ, .]**$ in units of [1/s]:

```wl
Overscript[ϵ, .]	 = 	(1/2)(∇**v** + (∇**v**)^T)
```

Note that this strain rate looks very much like the [infinitesimal strain measure](https://reference.wolfram.com/language/PDEModels/tutorial/StructuralMechanics/SolidMechanics.en.md#1798204897) defined in solid mechanics. In solid mechanics, the dependent variable $**u**$ is the displacement in [m] and $**ϵ**$ defines the infinitesimal strain. In fluid dynamics, the dependent variable $**v**$ is the velocity in [m/s] and $**Overscript[ϵ, .]**$ is a strain rate tensor. The strain rate tensor is also referred to as the rate or strain tensor.

The viscous stress tensor $**τ**$ in [Pa] for a Newtonian fluid is:

```wl
**τ**	 = 	2μ**Overscript[ϵ, .]** - (2/3)μ(∇·**v**)**I**
```

This gives a different formulation for the momentum equation:

```wl
ρ((∂**v**/∂t) + **v**·∇**v**) + ∇·(-**τ** + p**I**)	 = 	**F**
```

For non-Newtonian fluids, the stress strain rate relation is nonlinear. This nonlinearity can be expressed in terms of an apparent viscosity $Subscript[μ, app]$, where the apparent viscosity is a function of the shear rate $Overscript[γ, .]$ in [1/s]. Sometimes the apparent viscosity $Subscript[μ, app]$ is also called the effective viscosity.

For simplicity, but without restriction to generality, assume an incompressible non-Newtonian fluid, such that the $2 / 3μ(∇·**v**)**I**$ term is not present. This assumption leads to:

```wl
**τ**	 = 	2Subscript[μ, app](Overscript[γ, .])**Overscript[ϵ, .]**
```

The shear rate is computed as:

```wl
Overscript[γ, .]	 = 	Sqrt[2**Overscript[ϵ, .]** : **Overscript[ϵ, .]**]
```

where $ : $ is the contraction and computed by:

```wl
a : b	 = 	Subscript[∑, n]Subscript[∑, m]Subscript[a, nm]Subscript[b, nm]
```

A function to compute the shear rate is implanted as follows and is in the ``PDEModels``` context:

```
ShearRate[strainRate_] := 
	Sqrt[2 * TensorContract[strainRate * strainRate, {{1}, {2}}]];
```

The dynamic viscosity may depend on temperature and pressure, and for non-Newtonian fluids the dynamic viscosity additionally depends on the shear rate.

Different non-Newtonian constitutive models are used to express the apparent viscosity $Subscript[μ, app]$. The following non-Newtonian fluid flow models are implemented:

* Power law

* Carreau

* Bingham–Papanastasiou

* Herschel–Bulkley–Papanastasiou

* Custom apparent viscosity function

* Custom viscous stress tensor function

All models presented below the parameters ``pars`` need an entry for the ``"FluidDynamicsMaterialModel"`` :

```wl
pars = <|"FluidDynamicsMaterialModel" -> <|...|>, ...|>
```

Non-Newtonian fluids are also supported in the axisymmetric case.

Power Law

The power law model is a commonly used non-Newtonian fluid model.

The power law implements:

```wl
Subscript[μ, app]	 = 	m((max(Overscript[γ, .], Subscript[Overscript[γ, .], min])/Subscript[Overscript[γ, .], ref]))^n - 1
```

where $m$ is a factor that defaults to the dynamic viscosity $μ$. The reference shear rate $Subscript[Overscript[γ, .], ref]$ has a default of 1 [1/s]. The default minimal shear rate $Subscript[Overscript[γ, .], min]$ is set to $10^-2$ [1/s]. The scalar $n$ is an exponent. For $n = 1$, a Newtonian fluid is obtained, and as such, the power law implements a generalized Newtonian fluid model. For $n < 1$, the power law model results in a shear thinning fluid, also called a pseudo-plastic fluid. Paint, blood or ketchup falls in this category. For $n > 1$, the result is a shear thickening fluid, also called a dilatant fluid. An example for a shear thickening fluid is a mixture of water and cornstarch or quicksand.

The model given here deviates a bit from what is normally found in the literature. The original model predicts an infinite viscosity as the shear rate goes to zero. In real fluids, the shear rate tends to some constant value $Subscript[Overscript[γ, .], min]$. This shortcoming of the original formulation is avoided by the formulation given here. A disadvantage of the power law model is that it cannot describe a Newtonian plateau. The Carreau model does that.

Parameter names for the power law model:

```wl
In[4]:= "FluidDynamicsMaterialModel" -> <|"PowerLaw" -> <|"Exponent" -> n, "MinimalShearRate" -> Subscript[Overscript[γ, .], min], "ReferenceShearRate" -> Subscript[Overscript[γ, .], ref], "PowerLawViscosity" -> m|>|>;
```

The power law is often given as an expression of the kinematic viscosity:

```wl
ν	 = 	kOverscript[γ, .]^n - 1
```

where $k$ is called the consistency index.

An illustration of shear flow versus Newtonian flow is given below.

[image]

Compare the strain rate versus shear stress for shear thickening and thinning flows with Newtonian flows.

It can be seen that for shear thinning fluids, the apparent viscosity is reduced for an increased strain rate, and for a shear thickening fluid, the apparent viscosity is increased for an increased strain rate.

The actual implementation of the power law model for non-Newtonian flow is given below, and the way this model is constructed can be used as an example for constructing other non-Newtonian flow models.

Implementation of the power law non-Newtonian flow model:

```
ApparentViscosity["PowerLaw", vars_, pars_, data__]:=                           
Module[                                                                         
    {model, strainRate, shearRate, m, apparentViscosity, minShearRate, n, refShearRate},
	model = pars["FluidDynamicsMaterialModel"]["PowerLaw"];
	minShearRate = model["MinimalShearRate"];
	refShearRate = model["ReferenceShearRate"];
	n = model["Exponent"];
	m = model["PowerLawViscosity"];
	strainRate = "StrainRate" /. data;
	shearRate = PDEModels`ShearRate[strainRate];
	apparentViscosity = m * (Max[shearRate, minShearRate]/refShearRate)^(n-1);
	apparentViscosity
]
```

Power law example

The following section gives an example of a power law non-Newtonian flow.

Set up a widening region:

```wl
In[34]:= r = RegionUnion[Rectangle[{0, -1 / 2}, {5, 1 / 2}], Rectangle[{5, -3 / 2}, {30, 3 / 2}]];
```

Mesh and visualize the region:

```wl
In[35]:=
mesh = ToElementMesh[r, MaxCellMeasure -> 0.1];
mesh["Wireframe"]

Out[36]= [image]
```

In the given example, the kinematic viscosity is specified. Since the density is given as $1$, the kinematic viscosity equals the dynamic viscosity $ν = μ / ρ$.

To experiment, the power law exponent $n$ and the power law factor $m$ are set as symbolic parameters.

Set up variables, parameters and the parametric power law PDE:

```wl
In[8]:=
vars = {{u[x, y], v[x, y], p[x, y]}, {x, y}};
pars = <|"MassDensity" -> 1, 
	"FluidDynamicsMaterialModel" -> <|"PowerLaw" -> <|"Exponent" -> n, "PowerLawViscosity" -> m|>|>
	|>;
pde = FluidFlowPDEComponent[vars, pars] == {0, 0, 0};
```

The inflow profile is given as ${(1/2), 0}$, and at the outflow the pressure is set to $0$. The walls are no-slip walls.

Set up boundary conditions:

```wl
In[11]:=
bcs = {DirichletCondition[
	{u[x, y] == 1 / 2, v[x, y] == 0}, x == 0.], DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, 0 < x < 30], 
	DirichletCondition[p[x, y] == 0, x == 30]};
```

Create the parametric function:

```wl
In[12]:= pfun = ParametricNDSolveValue[{pde, bcs}, {u, v, p}, {x, y}∈mesh, {n, m}, Method -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}]

Out[12]= ParametricFunction[<>]
```

Evaluate the parametric model for various exponents $n$ and factors $m$ and measure the time it takes:

```wl
In[13]:= AbsoluteTiming[solutions = pfun@@@{{0.5, 0.008838835}, {1, 0.0125}, {1.5, 0.01767767}, {2, 0.025}};]

Out[13]= {9.99741, Null}
```

Extract the $x$-direction velocities:

```wl
In[14]:= uVelocities = solutions[[All, 1]];
```

Plot the flow profile from the middle of the channel to the top scaled to $1$ :

```wl
In[15]:= Plot[Evaluate[Through[uVelocities[28, y * 1.5]]], {y, 0, 1}]

Out[15]= [image]
```

Carreau

The Carreau model also implements a generalized Newtonian flow model and is useful to model polymer or blood flows. The Carreau model implements:

```wl
Subscript[μ, app]	 = 	Subscript[μ, ∞] + (Subscript[μ, 0] - Subscript[μ, ∞])(1 + (λOverscript[γ, .])^a)^(n - 1/a)
```

where $Subscript[μ, ∞]$ is the infinite shear rate viscosity, $Subscript[μ, 0]$ is the zero shear rate viscosity, $λ$ is a relaxation time with units [s], $n$ is a power law exponent, and $a$ is a transition exponent, which controls the transition between Newtonian and power law models. For $λ = 0$ or $n = 1$, this model reduces to the standard model of Newtonian fluid. For $λOverscript[γ, .]≪1$, the Carreau model reduces to a Newtonian flow. For $λOverscript[γ, .]≫1$, the Carreau model approaches the power law model.

The model presented here is a generalization of the Carreau model usually called the Carreau–Yasuda model. The original Carreau model can be obtained through the default setting of $a = 2$.

The model enables the description of Newtonian plateaus for very small and very large shear rates. A Newtonian plateau is a phenomenon that is associated with some non-Newtonian fluids. Such fluids show constant viscosities for very small shear rates. This is called the primary plateau. The secondary plateau can be observed for very high shear rates.

Parameter names for the Carreau model:

```wl
"FluidDynamicsMaterialModel" -> <| "Carreau" -> <|"PowerLawExponent" -> n, "TransitionExponent" -> a, "ZeroShearRateViscosity" -> Subscript[μ, 0], "InfiniteShearRateViscosity" -> Subscript[μ, ∞], "Lambda" -> λ|>|>;
```

The ``"TransitionExponent"`` parameter is optional. If not given a default value, $a = 2$ is used.

The actual implementation of the Carreau model for non-Newtonian flow is given below, and the way this model is constructed can be used as an example for constructing other non-Newtonian flow models.

Implementation of the Carreau non-Newtonian flow model:

```
ApparentViscosity["Carreau", vars_, pars_, data__] :=
 Module[
  	{model, mu0, muInf, n, a, lambda, strainRate, shearRate, apparentViscosity},
  	model = pars["FluidDynamicsMaterialModel"]["Carreau"];
  	mu0 = model["ZeroShearRateViscosity"];
  	muInf = model["InfiniteShearRateViscosity"];
  	n = model["PowerLawExponent"];
  	a = model["TransitionExponent"];
  	lambda = model["Lambda"];
  	strainRate = "StrainRate" /. data;
  	shearRate = ShearRate[strainRate + $MachineEpsilon];
  	apparentViscosity = muInf + (mu0 - muInf) *
     		(1 + (lambda * shearRate)^a )^((n - 1)/a);
  	apparentViscosity
  ]
```

A more detailed look at the parameters of the model could be useful. In what follows, several plots that show how the choice of parameters affects the apparent viscosity are presented. The model is considered with some fixed default parameters, the parameters are varied one by one, and the resulting apparent viscosities are compared. This will give a good idea how the choice of parameters affects the model. As the default choice of parameters, the following values are chosen:

* $Subscript[μ, 0] = 2$

* $Subscript[μ, ∞] = 1$

* $n = 1 / 2$

* $λ = 1$

* $a = 2$

Define the apparent viscosity function:

```wl
In[1]:= apparentViscosity = Subscript[μ, ∞] + (Subscript[μ, 0] - Subscript[μ, ∞]) * (1 + (λ * Overscript[γ, .]) ^ a) ^ ((n - 1) / a);
```

The first parameter to vary is $Subscript[μ, 0]$. This parameter is used to prescribe the viscosity of the fluid when the shear rate is equal to zero. The zero shear rate means that the fluid velocity is constant.

Set base and test parameter values $Subscript[μ, 0]$ :

```wl
In[2]:=
baseParameters = {Subscript[μ, ∞] -> 1, n -> 1 / 2, λ -> 1, a -> 2};
μ0Test = {3, 2, 3 / 2};
```

Plot a comparison of the apparent viscosities for different values of $Subscript[μ, 0]$ in the interval of [0,3]:

In[7]:= Plot[Evaluate[Table[(apparentViscosity/.baseParameters)/.Subscript[\[Mu], 0]->v,{v,\[Mu]0Test}]],{Overscript[\[Gamma], .],0,3},PlotLegends->(("Subscript[\[Mu], 0] = "<>ToString[N[\#]])&/@\[Mu]0Test),]

Out[7]=
Subscript[\[Mu], 0] = 3.
	Subscript[\[Mu], 0] = 2.
	Subscript[\[Mu], 0] = 1.5

It can be seen from the plot that $Subscript[μ, 0]$ has an impact on the apparent viscosity $Subscript[μ, app]$ for small values of the shear rate $Overscript[γ, .]$. Note that up until a shear rate of about 1/4, the plots exhibit a plateau, though not very pronounced. After that, $Subscript[μ, app]$ drops. For higher values of the shear rate $Overscript[γ, .]$, the impact of $Subscript[μ, 0]$ on the apparent viscosity $Subscript[μ, app]$ diminishes and becomes increasingly negligible, as can be seen in the plot below once the shear rate is expanded out further.

Plot a comparison of the apparent viscosities for different values of $Subscript[μ, 0]$ in the interval of [0, 100]:

In[5]:= Plot[Evaluate[Table[(apparentViscosity/.baseParameters)/.Subscript[\[Mu], 0]->v,{v,\[Mu]0Test}]],{Overscript[\[Gamma], .],0,100},PlotLegends->(("Subscript[\[Mu], 0] = "<>ToString[N[\#]])&/@\[Mu]0Test),]

Out[5]=
Subscript[\[Mu], 0] = 3.
	Subscript[\[Mu], 0] = 2.
	Subscript[\[Mu], 0] = 1.5

The secondary plateau is seen in the plot for the higher values of the shear rate $Overscript[γ, .]$. The zero shear rate viscosity $Subscript[μ, 0]$ also influences how quickly the second plateau is approached.

The parameter $Subscript[μ, ∞]$, on the other hand, describes how the viscosity behaves for higher values of the shear rate. It is the value to which the apparent viscosity converges if the shear rate goes to infinity.

Set base and test parameter values $Subscript[μ, ∞]$ :

```wl
In[6]:=
baseParameters = {Subscript[μ, 0] -> 2, n -> 1 / 2, λ -> 1, a -> 2};
μInfTest = {1.5, 1, 0.5};
```

Plot a comparison of the apparent viscosities for different values of $Subscript[μ, ∞]$ in the interval [0, 3]:

In[8]:= Plot[Evaluate[Table[(apparentViscosity/.baseParameters)/.Subscript[\[Mu], \[Infinity]]->v,{v,\[Mu]InfTest}]],{Overscript[\[Gamma], .],0,3},PlotLegends->(("Subscript[\[Mu], \[Infinity]] = "<>ToString[N[\#]])&/@\[Mu]InfTest),]

Out[8]=
Subscript[\[Mu], \[Infinity]] = 1.5
	Subscript[\[Mu], \[Infinity]] = 1.
	Subscript[\[Mu], \[Infinity]] = 0.5

The effect of $Subscript[μ, ∞]$ on the apparent viscosity $Subscript[μ, app]$ increases as the shear rate $Overscript[γ, .]$ increases. A better insight can be obtained by considering larger values of the shear rate.

Plot a comparison of the apparent viscosities for different values of $Subscript[μ, ∞]$ in the interval [0, 100]:

In[9]:= Plot[Evaluate[Table[(apparentViscosity/.baseParameters)/.Subscript[\[Mu], \[Infinity]]->v,{v,\[Mu]InfTest}]],{Overscript[\[Gamma], .],0,100},PlotLegends->(("Subscript[\[Mu], \[Infinity]] = "<>ToString[N[\#]])&/@\[Mu]InfTest),]

Out[9]=
Subscript[\[Mu], \[Infinity]] = 1.5
	Subscript[\[Mu], \[Infinity]] = 1.
	Subscript[\[Mu], \[Infinity]] = 0.5

The value of $Subscript[μ, ∞]$ plays a crucial role if the shear rate is large enough. It controls the value of the apparent viscosity at infinity.

The parameter $n$ represents the power law exponent. For values close to $1$, the model resembles Newtonian fluid behavior. For values of $n$ smaller than 1, the apparent viscosity goes to $Subscript[μ, ∞]$. For values of $n$ larger than 1, the model mimics the power law model and the apparent viscosity goes to \[ImplicitPlus]$±∞$, where the sign depends on the sign of $Subscript[μ, 0] - Subscript[μ, ∞]$.

Set base and test parameter values:

```wl
In[22]:=
baseParameters = {Subscript[μ, 0] -> 2, Subscript[μ, ∞] -> 1, λ -> 1, a -> 2};
nTest = {0.9, 0.5, 0.1};
```

Plot a comparison of the apparent viscosities for different values of $n$ in the interval of [0, 3]:

In[24]:= Plot[Evaluate[Table[(apparentViscosity/.baseParameters)/.n->v,{v,nTest}]],{Overscript[\[Gamma], .],0,3},PlotLegends->(("n = "<>ToString[N[\#]])&/@nTest),]

Out[24]=
n = 0.9
	n = 0.5
	n = 0.1

It can be seen that the parameter $n$ controls how fast the apparent viscosity increases or decreases. The effect of $n$ increases as the shear rate $Overscript[γ, .]$ increases.

Plot a comparison of the apparent viscosities for different values of $n$ in the interval of [0, 100]:

In[25]:= Plot[Evaluate[Table[(apparentViscosity/.baseParameters)/.n->v,{v,nTest}]],{Overscript[\[Gamma], .],0,100},PlotLegends->(("n = "<>ToString[N[\#]])&/@nTest),]

Out[25]=
n = 0.9
	n = 0.5
	n = 0.1

The plot shows that the power law exponent $n$ affects how fast the apparent viscosity $Subscript[μ, app]$ converges to $Subscript[μ, ∞]$. It should be noted that in the plots, all curves converge to $Subscript[μ, ∞] = 1$, even the curve associated with $n = 0.9$. The convergence is just slow because the value of $n$ is close to 1.

The parameter $λ$ is the relaxation time. It controls how fast the model approaches $Subscript[μ, ∞]$.

Set base and test parameter values $λ$ :

```wl
In[38]:=
baseParameters = {Subscript[μ, 0] -> 2, Subscript[μ, ∞] -> 1, n -> 1 / 2, a -> 2};
λTest = {10, 1, 0.2};
```

Plot a comparison of the apparent viscosities for different values of $λ$ in the interval [0, 3]:

In[40]:= Plot[Evaluate[Table[(apparentViscosity/.baseParameters)/.\[Lambda]->v,{v,\[Lambda]Test}]],{Overscript[\[Gamma], .],0,3},PlotLegends->(("\[Lambda] = "<>ToString[N[\#]])&/@\[Lambda]Test),]

Out[40]=
\[Lambda] = 10.
	\[Lambda] = 1.
	\[Lambda] = 0.2

The relaxation time $λ$ plays an important role even for very small values of the shear rate, unlike the power law index $n$.

Plot a comparison of the apparent viscosities for different values of $λ$ in the interval [0, 100]:

In[41]:= Plot[Evaluate[Table[(apparentViscosity/.baseParameters)/.\[Lambda]->v,{v,\[Lambda]Test}]],{Overscript[\[Gamma], .],0,100},PlotLegends->(("\[Lambda] = "<>ToString[N[\#]])&/@\[Lambda]Test),]

Out[41]=
\[Lambda] = 10.
	\[Lambda] = 1.
	\[Lambda] = 0.2

The relaxation time $λ$ controls how fast the apparent viscosity approaches $Subscript[μ, ∞]$. This is similar to the power law exponent $n$. The difference is that $λ$ is significant mostly for small values of the shear rate.

The last parameter $a$ controls how fast the model transitions from the Newtonian to the non-Newtonian behavior. Increasing the values of $a$ will increase the threshold for the shear rate at which the model transitions from Newtonian to non-Newtonian behavior.

Set base and test parameter values:

```wl
In[58]:=
baseParameters = {Subscript[μ, 0] -> 2, Subscript[μ, ∞] -> 1, n -> 1 / 2, λ -> 1};
aTest = {20, 2, 1};
```

Plot a comparison of the apparent viscosities for different values of $a$ in the interval of [0, 3]:

In[60]:= Plot[Evaluate[Table[(apparentViscosity/.baseParameters)/.a->v,{v,aTest}]],{Overscript[\[Gamma], .],0,3},PlotLegends->(("a = "<>ToString[N[\#]])&/@aTest),]

Out[60]=
a = 20.
	a = 2.
	a = 1.

It can be seen that the parameter $a$ controls the length of the primary plateau. But once the threshold is crossed and the fluid enters the non-Newtonian regime, then the apparent viscosities $Subscript[μ, app]$ quickly converge to the same specific value.

Plot a comparison of the apparent viscosities for different values of $a$ in the interval of [0, 100]:

In[61]:= Plot[Evaluate[Table[(apparentViscosity/.baseParameters)/.a->v,{v,aTest}]],{Overscript[\[Gamma], .],0,100},PlotLegends->(("a = "<>ToString[N[\#]])&/@aTest),]

```wl
Out[61]= [image]
```

The impact of $a$ for higher values of the shear rate is negligible.

Cross

The cross model is an empirical model and is useful for polymer melts or polymeric solutions. The Carreau model can be used to model a cross power law:

```wl
Subscript[μ, app]	 = 	Subscript[μ, ∞] + ((Subscript[μ, 0] - Subscript[μ, ∞])/(1 + λOverscript[γ, .])^p)
```

where $p = -(n - 1)$ that can be specified by setting $n = -p + 1$.

Bingham–Papanastasiou

The Bingham–Papanastasiou model is useful for viscoplastic material. The Bingham–Papanastasiou model implements:

```wl
Subscript[μ, app]	 = 	Subscript[μ, p] + (σ/Overscript[γ, .])(1 + E^-mOverscript[γ, .])
```

Parameter names for the Bingham–Papanastasiou model:

```wl
In[17]:= "FluidDynamicsMaterialModel" -> <|"Bingham-Papanastasiou" -> <|"PlasticViscosity" -> Subscript[μ, p], "YieldStress" -> σ, "ShearRateFactor" -> m|>|>;
```

The actual implementation of the Bingham–Papanastasiou model for non-Newtonian flow is given below, and the way this model is constructed can be used as an example for constructing other non-Newtonian flow models.

Implementation of the Bingham–Papanastasiou non-Newtonian flow model:

```
ApparentViscosity["Bingham-Papanastasiou", vars_, pars_, data__]:=              
Module[                                                                         
    {model, mup, sigma, m, strainRate, shearRate, apparentViscosity},
	model = pars["FluidDynamicsMaterialModel"]["Bingham-Papanastasiou"];
	mup = model["PlasticViscosity"];
	sigma = model["YieldStress"];
	m = model["ShearRateFactor"];
	strainRate = "StrainRate" /. data;
	shearRate = PDEModels`ShearRate[strainRate + $MachineEpsilon];
	apparentViscosity = mup + sigma/shearRate * (1 - Exp[-m * shearRate]);
	apparentViscosity
]
```

Herschel–Bulkley–Papanastasiou

The Herschel–Bulkley–Papanastasiou model implements:

```wl
Subscript[μ, app]	 = 	Subscript[μ, power law] + (σ/Overscript[γ, .])(1 + e^-mOverscript[γ, .])	 = 	m((Overscript[γ, .]/Subscript[Overscript[γ, .], ref]))^n - 1 + (σ/Overscript[γ, .])(1 + E^-mOverscript[γ, .])
```

The model uses a power law to compute the plastic viscosity of the Bingham–Papanastasiou model. All parameters for the power law and Bingham–Papanastasiou model can be specified, with the exception of the ``"PlasticViscosity"``, which is computed by the power law. The ``"MinimalShearRate"`` is set to $-∞$.

Parameter names for the Herschel–Bulkley–Papanastasiou model:

```wl
In[18]:= "FluidDynamicsMaterialModel" -> <|"Herschel-Bulkley-Papanastasiou" -> <|"Exponent" -> n, "ReferenceShearRate" -> Subscript[Overscript[γ, .], ref], "PowerLawViscosity" -> m, "YieldStress" -> σ, "ShearRateFactor" -> m|>|>;
```

Implementation of the Herschel–Bulkley–Papanastasiou non-Newtonian flow model:

```
ApparentViscosity["Herschel-Bulkley-Papanastasiou", vars_, pars_, data__]:=
Module[
{model, powerLawPars, mup, binghamPapanastasiouPars, apparentViscosity},
	model = pars["FluidDynamicsMaterialModel"]["Herschel-Bulkley-Papanastasiou"];
	powerLawPars = pars;
	powerLawPars["FluidDynamicsMaterialModel"] = <|"PowerLaw" -> model|>;
	mup = ApparentViscosity["PowerLaw", vars, powerLawPars, data];
	binghamPapanastasiouPars = pars;
	model["PlasticViscosity"] = mup;
	binghamPapanastasiouPars["FluidDynamicsMaterialModel"] =
		<|"Bingham-Papanastasiou"-> model|>;
	apparentViscosity = ApparentViscosity["Bingham-Papanastasiou", vars, 
		binghamPapanastasiouPars, data];
	apparentViscosity
]
```

Casson–Papanastasiou

The Casson–Papanastasiou model is useful for viscoplastic material or blood flow. The Casson–Papanastasiou model implements:

```wl
Subscript[μ, app] = (Sqrt[Subscript[μ, p]] + Sqrt[(τ/Overscript[γ, .])](1 + E^-Sqrt[mOverscript[γ, .]]))^2
```

Parameter names for the Casson–Papanastasiou model:

```wl
In[17]:= "FluidDynamicsMaterialModel" -> <|"Casson-Papanastasiou" -> <|"PlasticViscosity" -> Subscript[μ, p], "YieldStress" -> σ, "ShearRateFactor" -> m|>|>;
```

The actual implementation of the Casson–Papanastasiou model for non-Newtonian flow is given below, and the way this model is constructed can be used as an example for constructing other non-Newtonian flow models.

Implementation of the Casson–Papanastasiou non-Newtonian flow model:

```
ApparentViscosity["Casson-Papanastasiou", vars_, pars_, data__]:=
Module[
	{model, mup, tau, m, strainRate, shearRate, apparentViscosity},
	model = pars["FluidDynamicsMaterialModel"]["Casson-Papanastasiou"];
	mup = model["PlasticViscosity"];
	tau = model["YieldStress"];
	m = model["ShearRateFactor"];
	strainRate = "StrainRate" /. data;
	shearRate = ShearRate[strainRate + $MachineEpsilon];
	apparentViscosity = (Sqrt[mup] + Sqrt[tau/shearRate] * 
		(1 - Exp[-Sqrt[m * shearRate]]))^2;
	apparentViscosity
]
```

Custom Apparent Viscosity Model

If your favorite model is not in the list given so far, it can, nevertheless, be used by specifying a custom viscosity model.

Specifying a custom viscosity model:

```wl
In[19]:= "FluidDynamicsMaterialModel" -> <|"Custom" -> <|"ApparentViscosityFunction" -> apparentViscosity|>|>;
```

The function ``myApparentViscosity`` has the following signature:

```
apparentViscosity[_, vars_, pars_, data__]:=
Module[{},
....
]
```

The function ``myApparentViscosity`` should return a scalar value for the apparent viscosity. Other models given in this monograph can serve as a template for writing an apparent viscosity function.

Custom Viscous Stress Tensor

The entire viscous stress tensor $**τ**$ :

```wl
**τ**	 = 	2μ**Overscript[ϵ, .]** - (2/3)μ(∇·**v**)**I**
```

can be replaced by a function ``myViscousStressTensor``.

Specifying a custom viscous stress tensor:

```wl
In[1]:= "FluidDynamicsMaterialModelFunction" -> viscousStressTensor

Out[1]= "FluidDynamicsMaterialModelFunction" -> viscousStressTensor
```

The function ``myViscousStressTensor`` has the following signature:

```
viscousStressTensor[vars_, pars_, data__]:=
Module[{},
....
]
```

The function ``myViscousStressTensor`` should return a viscous stress tensor.

Axisymmetric Non-Newtonian Flow

Non-Newtonian fluids are also supported in the axisymmetric case.

Create a symbolic stationary axisymmetric fluid dynamics PDE of a power law fluid with a mass density of $ρ$, power law exponent $n$ and power law viscosity $m$ :

```wl
In[7]:=
FluidFlowPDEComponent[{{u[r, z], 0, w[r, z], p[r, z]}, {r, θ, z}}, <|"MassDensity" -> ρ, "RegionSymmetry" -> "Axisymmetric", 
	"FluidDynamicsMaterialModel" -> <|"PowerLaw" -> <|"Exponent" -> n, "PowerLawViscosity" -> m|>|>
	|>]
```

Computing viscosity and stress

The solution of the fluid dynamics equations can be used to obtain additional information about the fluid flow. This includes examples like computing the traction acting on the wall or computing the [apparent viscosity](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1356543500) of a non-Newtonian fluid. In order to compute the [viscosity](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1356543500) or the [viscous stress tensor](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1904112883), the functions ``FluidViscosity`` and ``FluidViscousStress`` may be used.

The following example shows how to make use of ``FluidViscosity`` and ``FluidViscousStress`` to analyze fluid flows.

Set up a domain:

```wl
In[5]:=
domain = RegionUnion[Rectangle[{0.5, -1}, {1.5, 5}], Disk[{1, 5}, 1], Rectangle[{1, 5.5}, {5, 4.5}], Disk[{5, 5}, 1], 
	Rectangle[{1, 0.75}, {5, 1.25}], Disk[{1, 1}, 1], 
	Rectangle[{4.5, 0.5}, {5.5, 5}], Disk[{5, 1}, 1], Rectangle[{5.5, 0.5}, {9, 1.5}], Disk[{9, 1}, 1], Rectangle[{5.5, 4.75}, {9, 5.25}], Disk[{9, 5}, 1], 
	Rectangle[{8.5, 1.5}, {9.5, 7}]]

Out[5]= [image]
```

Create and visualize the mesh:

```wl
In[6]:=
mesh = ToElementMesh[domain, MaxCellMeasure -> 0.1];
mesh["Wireframe"]

Out[7]= [image]
```

Set up variables, parameters and a power law material model:

```wl
In[8]:=
vars = {{u[x, y], v[x, y], p[x, y]}, {x, y}};
pars = <|"MassDensity" -> 1, "FluidDynamicsMaterialModel" -> <|"PowerLaw" -> <|"Exponent" -> 2, "PowerLawViscosity" -> 0.01|>|>|>;
pde = FluidFlowPDEComponent[vars, pars] == {0, 0, 0};
```

Set up boundary conditions:

```wl
In[11]:=
bcWall = DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, y > -1 && y < 7];
bcInlet = DirichletCondition[{u[x, y] == 0, v[x, y] == 2 * ((x - 0.5) - (x - 0.5) ^ 2)}, y == -1];
bcPressure = DirichletCondition[p[x, y] == 0, x == 9.5 && y == 7];
bcs = {bcWall, bcInlet, bcPressure};
```

Solve the equations:

```wl
In[15]:= {xVel, yVel, pressure} = NDSolveValue[{pde, bcs}, {u[x, y], v[x, y], p[x, y]}, {x, y}∈mesh, Method -> {"PDEDiscretization" -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}}]

Out[15]=
{InterpolatingFunction[{{0.011169173774871454, 10.}, {-1., 7.}}, 
 {5, 4225, 0, {1573, 0}, {3, 0}, 0, 0, 0, 0, Indeterminate & , {}, {}, False}, 
 {NDSolve`FEM`ElementMesh[CompressedData["«26170»"], 
   {NDSolve`FEM`TriangleElement[CompressedData[" ... 1, 31, 31, 33, 33, 33, 33, 35, 35, 35, 35, 29, 29, 29, 29, 23, 23, 23, 
      23, 25, 25, 25, 25, 27, 27, 27, 27, 21, 21, 21, 21, 22, 1, 36, 13, 22, 1, 36, 13, 15, 28, 3, 
      30, 23, 12, 2, 35}]}]}, CompressedData["«4750»"], {Automatic}][x, y]}
```

Visualize the fluid velocity using ``VectorPlot`` :

```wl
In[16]:= VectorPlot[Evaluate[{xVel, yVel}], {x, y}∈mesh, PlotLegends -> Automatic, VectorPoints -> Fine]

Out[16]= [image]
```

Visualize the pressure using ``DensityPlot`` :

```wl
In[17]:= DensityPlot[pressure, {x, y}∈mesh, PlotLegends -> Automatic]

Out[17]= [image]
```

Compute the viscosity using the ``FluidViscosity`` :

```wl
In[26]:= viscosity = FluidViscosity[vars, pars, {xVel, yVel}]

Out[26]=
InterpolatingFunction[{{0.011169173774871454, 10.}, {-1., 7.}}, 
 {5, 4225, 0, {1573, 0}, {3, 0}, 0, 0, 0, 0, Indeterminate & , {}, {}, False}, 
 {NDSolve`FEM`ElementMesh[CompressedData["«26168»"], 
   {NDSolve`FEM`TriangleElement[CompressedData["« ... , 35, 35, 35, 36, 36, 36, 36, 29, 
      29, 29, 29, 22, 1, 36, 13, 22, 1, 36, 13, 15, 28, 3, 30, 23, 12, 2, 35}]}, 
   {NDSolve`FEM`PointElement[CompressedData["«904»"], CompressedData["«418»"]]}]}, 
 CompressedData["«16762»"], {Automatic}][x, y]
```

``FluidViscosity`` computes the [viscosity](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#1356543500) of the fluid. The precise formula for the apparent viscosity of non-Newtonian fluids depends on the model. In the case of the power-law model, which this example uses, it is given by equation (\!\(\*CounterBox["NumberedEquation", "eqn: power law"]\)).

In order to compute viscosity or stress from the solution, the [strain rate](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#657917637) has to be computed first. This means that the strain rate is computed twice if both viscosity and stress are desired. In order to speed up the computations, the strain rate can be precomputed first and then passed to both ``FluidViscosity`` and ``FluidViscousStress``.

Compute the strain rate:

```wl
In[27]:= strainRate = PDEModels`StrainRateMeasure[{{xVel, yVel, pressure}, vars[[2]]}, pars]

Out[27]=
SymmetrizedArray[StructuredArray`StructuredData[{2, 2}, 
  {{{1, 1} -> InterpolatingFunction[{{0.011169173774871454, 10.}, {-1., 7.}}, 
       {5, 4225, 0, {1573, 0}, {3, 0}, 0, 0, 0, 0, Indeterminate & , {}, {}, False}, 
       {NDSolve`FEM`Elemen ... 9, 29, 29, 22, 1, 36, 13, 22, 1, 36, 13, 15, 
            28, 3, 30, 23, 12, 2, 35}]}, {NDSolve`FEM`PointElement[CompressedData["«908»"], 
           CompressedData["«422»"]]}]}, CompressedData["«17012»"], {Automatic}][x, y]}, Symmetric[{1, 2}]}]]
```

Compute the viscosity using the ``FluidViscosity`` using the precomputed strain rate:

```wl
In[28]:= viscosity = FluidViscosity[vars, pars, strainRate]

Out[28]=
InterpolatingFunction[{{0.011169173774871454, 10.}, {-1., 7.}}, 
 {5, 4225, 0, {1573, 0}, {3, 0}, 0, 0, 0, 0, Indeterminate & , {}, {}, False}, 
 {NDSolve`FEM`ElementMesh[CompressedData["«26168»"], 
   {NDSolve`FEM`TriangleElement[CompressedData["« ... , 35, 35, 35, 36, 36, 36, 36, 29, 
      29, 29, 29, 22, 1, 36, 13, 22, 1, 36, 13, 15, 28, 3, 30, 23, 12, 2, 35}]}, 
   {NDSolve`FEM`PointElement[CompressedData["«904»"], CompressedData["«418»"]]}]}, 
 CompressedData["«16762»"], {Automatic}][x, y]
```

Plot the viscosity using ``DensityPlot`` :

```wl
In[29]:= DensityPlot[viscosity, {x, y}∈mesh, PlotRange -> All, PlotLegends -> Automatic]

Out[29]= [image]
```

Compute the viscous stress with ``FluidViscousStress`` :

```wl
In[30]:= viscousStress = FluidViscousStress[vars, pars, {xVel, yVel}]

Out[30]=
SymmetrizedArray[StructuredArray`StructuredData[{2, 2}, 
  {{{1, 1} -> InterpolatingFunction[{{0.011169173774871454, 10.}, {-1., 7.}}, 
       {5, 4225, 0, {1573, 0}, {3, 0}, 0, 0, 0, 0, Indeterminate & , {}, {}, False}, 
       {NDSolve`FEM`Elemen ... 29, 22, 1, 36, 13, 22, 1, 36, 13, 15, 
            28, 3, 30, 23, 12, 2, 35}]}, {NDSolve`FEM`PointElement[CompressedData["«908»"], 
           CompressedData["«422»"]]}]}, CompressedData["«17074»"], {Automatic}][x, 
      y]}, Symmetric[{1, 2}]}]]
```

``FluidViscousStress`` computes the viscous stress tensor. The formula for the viscous stress tensor usually has has the form of the equation (\!\(\*CounterBox["NumberedEquation", "eqn: viscous stress"]\)).

``FluidViscousStress`` also supports the strain rate as a parameter.

Compute the viscous stress with ``FluidViscousStress`` using the precomputed strain rate:

```wl
In[32]:= viscousStress = FluidViscousStress[vars, pars, strainRate]

Out[32]=
SymmetrizedArray[StructuredArray`StructuredData[{2, 2}, 
  {{{1, 1} -> InterpolatingFunction[{{0.011169173774871454, 10.}, {-1., 7.}}, 
       {5, 4225, 0, {1573, 0}, {3, 0}, 0, 0, 0, 0, Indeterminate & , {}, {}, False}, 
       {NDSolve`FEM`Elemen ... 29, 22, 1, 36, 13, 22, 1, 36, 13, 15, 
            28, 3, 30, 23, 12, 2, 35}]}, {NDSolve`FEM`PointElement[CompressedData["«908»"], 
           CompressedData["«422»"]]}]}, CompressedData["«17074»"], {Automatic}][x, 
      y]}, Symmetric[{1, 2}]}]]
```

Plot the norm of the viscous stress using ``DensityPlot`` :

```wl
In[33]:= DensityPlot[Norm[viscousStress, "Frobenius"], {x, y}∈mesh, PlotLegends -> Automatic]

Out[33]= [image]
```

The total stress tensor of viscous fluid is equal to

Subscript[**τ**, tot]	 = 	 - p**I** + Subscript[**τ**, vis]

where $Subscript[τ, tot]$ is the total stress tensor, $Subscript[τ, vis]$ is the viscous stress tensor, $**I**$ is the identity matrix and $p$ is the pressure.

Compute the total stress:

```wl
In[34]:= totalStress = -pressure * IdentityMatrix[2] + viscousStress;
```

Plot the norm of the total stress using ``DensityPlot`` :

```wl
In[35]:= DensityPlot[Norm[totalStress, "Frobenius"], {x, y}∈mesh, PlotLegends -> Automatic]

Out[35]= [image]
```

Viscoelastic Flow

In viscoelastic flow, the viscous stress tensor $**τ**$ in [Pa] for incompressible flow is given by:

```wl
**τ**	 = 	2μ**S** + Subscript[**T**, ve]
```

where $**Subscript[T, ve]**$ is a viscoelastic stress tensor. The equations involving viscoelastic flow can quickly become numerically unstable, and the usage of stabilization techniques is necessary. Since this version of ``NDSolve`` does not offer stabilization techniques, viscoelastic flow is not implemented currently.

Energy Transport—Nonisothermal Flow

Often one wants to get an understanding of how temperature differences affect the flow of a fluid or how heat is transported with the fluid flow. To achieve this, the Navier–Stokes equations are coupled with a heat transfer equation:

```wl
ρ((∂**v**/∂t) + **v**·∇**v**) - ∇·(μ(∇**v** + (∇**v**)^T) - (2/3)μ(∇·**v**)**I**) + ∇p - **F**	 = 	0

(∂ρ/∂t) + ∇·(**v** ρ)	 = 	0

ρ Subscript[C, p](∂T/∂t) + ∇·(-k ∇T) + ρ Subscript[C, p]**v**·∇T - Q	 = 	0
```

where $T$ is absolute temperature in [$K$], $Subscript[C, p]$ the specific heat capacity at a specific pressure in [$J / (kg K)$], $k$ the thermal conductivity in [$W / (m K)$] and $Q$ the heat source densities, or if you will, energy densities, in [$W / m^3$]. A detailed explanation and derivation of the heat equation can be found in the [Heat Transfer Monograph](https://reference.wolfram.com/language/PDEModels/tutorial/HeatTransfer/HeatTransfer.en.md). Among the heat sources could be a viscous heating $**τ** : **Overscript[ϵ, .]**$ or terms that describe work done by changes in pressure. Both are often negligible but can be added if need be.

The key point is that a change in temperature induces a change in density. The volume increases and thus the weight decreases and the fluid rises: a temperature-dependent buoyancy force arises.

Boussinesq Approximation

A common simplification to model the coupling of the heat and Navier–Stokes equations is the Boussinesq approximation, which is based on the following assumptions [2, p.126]:

* The density $ρ$ is constant except in the buoyancy term.

* All other fluid properties are assumed constant.

* Viscous dissipation is negligible.

The first item allows the incompressible form of the continuity equation to be used:

```wl
∇·**v** 	 = 	0
```

The momentum equation is then given as:

```wl
Subscript[ρ, ref]((∂**v**/∂t) + **v**·∇**v**) + ∇·(-μ(∇**v** + (∇**v**)^T) + (2/3)μ(∇·**v**)**I** + p**I**) - ρ(T)**g**	 = 	0
```

Note the introduction of $Subscript[ρ, ref]$, and $**g**$ is the gravitational field. The product of $ρ(T) **g**$ is the buoyancy contribution. As a further simplification, the Boussinesq approximation assumes a linear relation between the density $ρ$ and temperature $T$ :

```wl
ρ(T) 	 = 	Subscript[ρ, ref](1 - α(T - Subscript[T, ref]))
```

where $α$ is the thermal expansion coefficient and $Subscript[ρ, ref]$ and $Subscript[T, ref]$ are reference values. In fact, the equation is a first-order Taylor expansion around $Subscript[T, ref]$. It should be noted that the linear relationship between $ρ$ and $T$ needs to be replaced with a close-to-quadratic relationship for special cases, like pure water at $4^o C$ [4].

The Boussinesq approximation provides a coupling between the flow equations and the heat transfer model. A temperature difference is now capable of driving the fluid flow: If a flow field is isothermal, that means $T = Subscript[T, ref]$; then the fluid will remain at rest. If, however, one wall is heated, a buoyancy force will drive the fluid flow.

Some caution needs to be exercised when making use of the Boussinesq approximation: The only material parameter that is varying is the density $ρ(T)$ in the buoyancy term. That means all other occurrences of the density and also the dynamic viscosity $μ$, the thermal expansion coefficient $α$, the specific heat capacity $Subscript[C, p]$ and the thermal conductivity $k$ are constant in the interval $T∈[Subscript[T, cool], Subscript[T, hot]]$. The Boussinesq approximation assumes that the pressure is constant. The density in the buoyancy term is a Taylor expansion where $Subscript[T, ref]∈[Subscript[T, cool], Subscript[T, hot]]$ is the expansion point.

The restrictions pointed out above are quite severe, and it is advisable to make sure that the conditions are satisfied. Some effort was made to derive the following criteria [2; 5]:

```wl
\[LeftBracketingBar](α\[LeftBracketingBar]**g**\[RightBracketingBar]LSubscript[T, ref]/Subscript[C, p]ΔT)\[RightBracketingBar]⩽(1/10)

\[LeftBracketingBar](α\[LeftBracketingBar]**g**\[RightBracketingBar]L/Subscript[C, p])\[RightBracketingBar]Pr = \[LeftBracketingBar](α\[LeftBracketingBar]**g**\[RightBracketingBar]L/Subscript[C, p])\[RightBracketingBar]((μ/Subscript[ρ, ref])) / α⩽(1/10)

\[LeftBracketingBar]αΔT\[RightBracketingBar] = \[LeftBracketingBar](1/ρ)(dρ/dT)ΔT\[RightBracketingBar]⩽(1/10)

\[LeftBracketingBar](1/μ)(dμ/dT)ΔT\[RightBracketingBar]⩽(1/10)

\[LeftBracketingBar](1/Subscript[C, p])(dSubscript[C, p]/dT)ΔT\[RightBracketingBar]⩽(1/10)

\[LeftBracketingBar](1/κ)(dκ/dT)ΔT\[RightBracketingBar]⩽(1/10)

\[LeftBracketingBar](1/α)(dα/dT)ΔT\[RightBracketingBar]⩽(1/10)
```

Here, $ΔT≈(Subscript[T, hot] - Subscript[T, cold]) / 2$ is the maximal temperature fluctuation around $Subscript[T, ref] = (Subscript[T, hot] + Subscript[T, cold]) / 2$, and $L$ is a characteristic length.

The following example illustrates this procedure for water [2, p. 133]. To make an estimate of when the approximation is valid, the inequalities from above are encoded.

Encode the Boussinesq approximation validity inequalities:

```wl
In[54]:=
ineq1 = α * g * L * Subscript[T, ref] / (Subscript[C, p] * ΔT) <= 1 / 10 && ΔT > 0;
ineq2 = α * g * L / Subscript[C, p] * Pr <= 1 / 10 && ΔT > 0;
ineq3 = α * ΔT <= 1 / 10 && ΔT > 0;
ineq4 = 1 / μ[T] D[μ[T], T] * ΔT <= 1 / 10;
ineq5 = 1 / Subscript[C, p][T] D[Subscript[C, p][T], T] * ΔT <= 1 / 10;
ineq6 = 1 / κ[T] D[κ[T], T] * ΔT <= 1 / 10;
ineq7 = 1 / α[T] D[α[T], T] * ΔT <= 1 / 10;
```

The water at $Subscript[T, ref] = 288 [K]$ and $Subscript[p, ref] = 1.01325 * 10^5 [Pa] = 1 [atm]$ has the following parameters [5]:

|                 |     |          |                     |                       |                     |                     |
| --------------- | --- | -------- | ------------------- | --------------------- | ------------------- | ------------------- |
| Cp [J / (Kg K)] | Pr  | α [K^-1] | (1/μ)(∂μ/∂T) [K^-1] | (1/Cp)(∂Cp/∂T) [K^-1] | (1/κ)(∂κ/∂T) [K^-1] | (1/α)(∂α/∂T) [K^-1] |
| 4200            | 8.1 | 15 10^-3 | -27 10^-3           | -24 10^-3             | 17 10^-4            | 8 10^-2             |

Set up the values as a list of rules:

```wl
In[61]:= valueRules = {Subscript[C, p] -> 4200, Pr -> 81 / 10, α -> 15 * 10 ^ -3, 1 / μ[T] D[μ[T], T] -> 27 * 10 ^ -3, 1 / Subscript[C, p][T] D[Subscript[C, p][T], T] -> 24 * 10 ^ -3, 1 / κ[T] D[κ[T], T] -> 14 * 10 ^ -4, 1 / α[T] D[α[T], T] -> 8 * 10 ^ -2, g -> 981 / 100, Subscript[T, ref] -> 288};
```

Write a helper function to simplify the insertion of values into the inequalities:

```wl
In[62]:= fun[ineq_] := Simplify[Reduce[ineq /. valueRules, {ΔT, L}]//N, Assumptions -> ΔT > 0 && L > 0]
```

Reduce the inequalities:

```wl
In[63]:= reducedInequalities = fun /@ {ineq1, ineq2, ineq3, ineq4, ineq5, ineq6, ineq7}

Out[63]= {L ≤ 9.91052 ΔT, L ≤ 352.374, ΔT ≤ 6.66667, ΔT ≤ 3.7037, ΔT ≤ 4.16667, ΔT ≤ 71.4286, ΔT ≤ 1.25}
```

Visualize the area spanned by $L$ and $ΔT$ values of water for which the Boussinesq approximation is valid:

```wl
In[64]:= RegionPlot[And@@reducedInequalities, {L, 10 ^ -2, 15}, {ΔT, 0.01, 1.4}, FrameLabel -> {"L [m]", "ΔT [K]"}, PlotPoints -> 31]

Out[64]= [image]
```

The same can also be done considering the pressure [5]. One final remark is that for integrating quantities that involve the material constants, the reference values like $Subscript[ρ, ref]$ should be used and not the temperature-dependent density $ρ(T)$ used in the buoyancy force term [4].

Rayleigh–Bénard Convection

A Rayleigh–Bénard convection is presented as an application of the Boussinesq approximation. In this example, a rectangular region is filled with a fluid. The side walls are insulated and the top and bottom walls are cooled and heated, respectively. First, a region and some parameters are specified.

The rectangle has an aspect ratio of 4:1.

Specify parameters and a rectangular region:

```wl
In[29]:=
sizes = {length -> 4, height -> 1};
Ω = Rectangle[{0, 0}, {length, height}] /. sizes;
```

Visualize the region:

```wl
In[31]:= RegionPlot[Ω, AspectRatio -> Automatic]

Out[31]= [image]
```

Create and visualize the mesh:

```wl
In[32]:=
mesh = ToElementMesh[Ω, MaxCellMeasure -> 0.002];
mesh["Wireframe"]

Out[33]= [image]
```

Set up Navier–Stokes equations that are coupled to a heat equation making use of a Boussinesq approximation. Use the material parameters specified.

Set up fluid and heat transfer variables:

```wl
In[9]:=
varsFluid = {{u[t, x, y], v[t, x, y], p[t, x, y]}, t, {x, y}};
varsHeat = {T[t, x, y], t, {x, y}};
```

Specify material parameters for air:

```wl
In[15]:= parameters = {Subscript[T, cold] -> 292.5, Subscript[T, hot] -> 293.5, Subscript[T, ref] -> 293, Subscript[ρ, ref] -> 1.205, μ -> 1.81 * 10 ^ -4, α -> 3.4 * 10 ^ -3, Subscript[g, z] -> -9.81, Subscript[c, p] -> 29.13, k -> 0.02586};
```

Set up the heat equation and laminar flow parameters:

```wl
In[16]:= pars = <|"MassDensity" -> Subscript[ρ, ref], "DynamicViscosity" -> μ, "HeatConvectionVelocity" -> {u[t, x, y], v[t, x, y]}, "SpecificHeatCapacity" -> Subscript[c, p], "ThermalConductivity" -> k|>;
```

Set up the Boussinesq approximation:

```wl
In[17]:= boussinesq = If[t < 1, t, 1] * Subscript[ρ, ref] * (1 - α * (T[t, x, y] - Subscript[T, ref])) * Subscript[g, z];
```

Set up Navier–Stokes equations coupled to a heat equation:

```wl
In[18]:= model = Flatten[{FluidFlowPDEComponent[varsFluid, pars] - {0, boussinesq, 0}, HeatTransferPDEComponent[varsHeat, pars]}]

Out[18]= {{Subscript[ρ, ref] u[t, x, y], Subscript[ρ, ref] v[t, x, y]}.Inactive[Grad][u[t, x, y], {x, y}] + Inactive[Div][-2*μ*{{0, 1/4}, {1/4, 0}} . Inactive[Grad][v[t, x, y], {x, y}], {x, y}] + Inactive[Div][-2*μ*{{1, 0}, {0, 1/2}} . Inactive[Grad][u[t, x ... bscript[ρ, ref] u[t, x, y], Subscript[c, p] Subscript[ρ, ref] v[t, x, y]}.Inactive[Grad][T[t, x, y], {x, y}] + Inactive[Div][{{-k, 0}, {0, -k}} . Inactive[Grad][T[t, x, y], {x, y}], {x, y}] + Subscript[c, p] Subscript[ρ, ref] T^(1, 0, 0)[t, x, y]}
```

Boundary conditions and initial conditions need to be set up.

Set up no-slip boundary conditions for the velocities on all boundary walls:

```wl
In[19]:= wall = DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, True];
```

Set up a reference pressure point:

```wl
In[20]:= reference = DirichletCondition[p[t, x, y] == 0, x == 0 && y == 0];
```

Specify a temperature difference between the top and bottom walls:

```wl
In[21]:= temperatures = {DirichletCondition[T[t, x, y] == Subscript[T, hot], y == 0], DirichletCondition[T[t, x, y] == Subscript[T, cold], y == height]};
```

Replace parameters in the boundary conditions:

```wl
In[22]:= bcs = {wall, reference, temperatures} /. sizes;
```

Set up initial conditions such that the system is at rest:

```wl
In[23]:= ic = {u[0, x, y] == 0, v[0, x, y] == 0, p[0, x, y] == 0, T[0, x, y] == Subscript[T, ref]};
```

Monitor the time integration progress and the total time it takes to solve the PDE while using a refined mesh and interpolating the velocities $u$ and $v$ and the temperature $T$ with second order and the pressure $p$ with first order.

Monitor the time integration of the equations:

```wl
In[24]:=
tEnd = 500;
pde = {model == {0, 0, 0, 0}, bcs, ic} /. parameters;
Monitor[AbsoluteTiming[
	result = NDSolveValue[pde, {u, v, p, T}, {x, y}∈mesh, {t, 0, tEnd}, 
	DependentVariables -> {u, v, p, T}, 
	Method -> {
	"PDEDiscretization" -> {"MethodOfLines", 
	"SpatialDiscretization" -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1, T -> 2}}}}
	, EvaluationMonitor :> (currentTime = Row[{"t = ", CForm[t]}])];
	], currentTime]

Out[26]= {324.242, Null}
```

The last step is post-processing.

Split the result according to the dependent variables:

```wl
In[27]:= {uVel, vVel, pressure, temperature} = result;
```

Construct a visualization of the region's boundary:

```wl
In[30]:= bmeshGr = uVel["ElementMesh"]["Edgeframe"]

Out[30]= [image]
```

Visualize the pressure distribution:

```wl
In[31]:= Show[ContourPlot[pressure[tEnd, x, y], {x, y}∈pressure["ElementMesh"], ...], bmeshGr]

Out[31]= [image]
```

Visualize the temperature distribution:

```wl
In[32]:= Show[ContourPlot[temperature[tEnd, x, y], {x, y}∈temperature["ElementMesh"], ..., PlotRange -> All], bmeshGr]

Out[32]= [image]
```

Visualize the velocity field:

```wl
In[33]:= Show[bmeshGr, VectorPlot[{uVel[tEnd, x, y], vVel[tEnd, x, y]}, {x, y}∈uVel["ElementMesh"], ...]]

Out[33]= [image]
```

Animate the change in temperature and the velocity streamlines:

```wl
In[34]:=
numberOfFrames = 30;
frames = Table[Show[bmeshGr, 
	ContourPlot[temperature[t, x, y], {x, y}∈temperature["ElementMesh"], ...], 
	StreamPlot[{uVel[t, x, y], vVel[t, x, y]}, {x, y}∈uVel["ElementMesh"], ...]
	], {t, 0, tEnd, tEnd / numberOfFrames}];ListAnimate[frames, SaveDefinitions -> True]

Out[35]= DynamicModule[«8»]
```

Further examples for heat transfer with fluid flow can be found in the model collection example [Buoyancy-Driven Flow](https://reference.wolfram.com/language/PDEModels/tutorial/Multiphysics/ModelCollection/BuoyancyDrivenFlow.en.md), where a Reynolds number–based model is used, and a second example, [Heat Exchange](https://reference.wolfram.com/language/PDEModels/tutorial/Multiphysics/ModelCollection/HeatExchanger.en.md), uses a standard model. The example [Thermal Decomposition](https://reference.wolfram.com/language/PDEModels/tutorial/Multiphysics/ModelCollection/ThermalDecomposition.en.md) shows a coupling of heat, mass and fluid transport.

Flow Boundary Conditions

The geometrical domain of a fluid flow model is finite. At the boundaries of the domain, boundary conditions are needed to model the surroundings in which the flow is embedded. Boundary conditions describe domain walls or inflows into or outflows out of a domain. To make it easier to talk about boundary conditions, it is useful to introduce some terminology. Let $Subscript[ϑ, n]$ denote a velocity component normal to a boundary and $Subscript[ϑ, t]$ denote a velocity component tangential to a boundary. $∂Subscript[ϑ, n] / ∂n$ and $∂Subscript[ϑ, t] / ∂n$ denote the derivatives of both components in the normal direction.

Consider for simplicity a 2D geometry that is parallel to the $x$ and $y$ axes, like the lid-driven cavity example from the introduction. There are two cases: upright boundaries and horizontal boundaries. For the upright boundary, there is:

```wl
|        |
| ------ |
| ϑn = u |
| ϑt = v |
```

and

```wl
|                    |
| ------------------ |
| ∂ϑn / ∂n = ∂u / ∂x |
| ∂ϑt / ∂n = ∂v / ∂x |
```

[image]

At a upright boundary, marked in light red, a normal fluid flow component is $Subscript[ϑ, n] = u$ and tangential fluid flow component is $Subscript[ϑ, t] = v$.

For horizontal boundaries, there is:

```wl
|        |
| ------ |
| ϑn = v |
| ϑt = u |
```

and

```wl
|                    |
| ------------------ |
| ∂ϑn / ∂n = ∂v / ∂y |
| ∂ϑt / ∂n = ∂u / ∂y |
```

[image]

At a horizontal boundary, marked in light red, a normal fluid flow component is $Subscript[ϑ, n] = v$ and tangential fluid flow component is $Subscript[ϑ, t] = u$.

In cases where the boundary is not axis-aligned, the components $Subscript[ϑ, n]$ and $Subscript[ϑ, t]$ need to be computed from $u$ and $v$.

Inflow Conditions

An inflow is also called an inlet. Both velocity components are prescribed:

```wl
|          |
| -------- |
| ϑn = ϑn0 |
| ϑt = ϑt0 |
```

Here $Subscript[ϑ, n]$ and $Subscript[ϑ, t]$ are the normal and tangential components of the velocity vector $**v**$. To specify the inflow condition, $Subsuperscript[ϑ, n, 0]$ and $Subsuperscript[ϑ, t, 0]$ are given. A typical inflow profile is:

```wl
|        |     |         |
| ------ | --- | ------- |
| **v**n |  =  | finflow |
| **v**t |  =  | 0       |
```

where $Subscript[f, inflow]$ is a function that specifies a fluid velocity in the normal direction.

One thing to keep in mind is that the inflow condition needs to be consistent with adjacent boundary conditions. For example, if the inflow profile at the boundary of the inlet is 0, then a non-slip wall boundary can be used. If the inflow profile is nonzero at the boundary of the inlet, then using an adjacent no-slip wall boundary condition will make it hard, if not impossible, for the solver to converge the model. In some cases, when the inflow condition at the boundary is nonzero, it may make sense to use a slip wall boundary condition instead.

An easy way to create a parabolic inflow profile is by using ``Fit``. Say there is an inflow from $1 <= x <= 2$ and a maximum inflow velocity of 1 [$Quantity[1, ("Meters"/"Seconds")]$] is given.

Create data points with $0$ velocity at the endpoints and 1 [$Quantity[1, ("Meters"/"Seconds")]$] at the middle:

```wl
In[1]:= fitData = {{1, 0}, {(1 + 2) / 2, 1}, {2, 0}};
```

Create a parabolic inflow profile:

```wl
In[3]:= parabolicInflowProfile = Fit[fitData, {1, x, x ^ 2}, x]

Out[3]= -8. + 12. x - 4. x^2
```

Visualize the parabolic inflow profile:

```wl
In[4]:= Plot[parabolicInflowProfile, {x, 1, 2}]

Out[4]= [image]
```

Outflow Conditions

An outflow is also called an outlet. Both velocity components remain constant with respect to the normal:

```wl
|              |
| ------------ |
| (∂ϑn/∂n) = 0 |
| (∂ϑt/∂n) = 0 |
```

Wall Conditions

A wall is any impermeable boundary containing a fluid. One condition at the wall boundary is then that the fluid cannot flow through the wall. A minimal requirement is that:

```wl
Subscript[ϑ, n] = 0
```

If only this minimal requirement is satisfied, the condition is called a slip condition. A second requirement is often that additionally there should be no fluid flow tangential to the wall boundary. In this second case, the condition is called a no-slip condition.

Pay attention that wall boundary conditions are consistent with inflow boundary conditions. For more details, see the section on [Inflow Conditions](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/LaminarFlow.en.md#308138130).

No-slip

A typical wall boundary condition is the no-slip condition. No-slip means that first, there is no fluid going through the wall and second, that the fluid adheres to the wall. The second condition models friction at the boundary. In sum, this means that:

```wl
|        |
| ------ |
| ϑn = 0 |
| ϑt = 0 |
```

Thus, the fluid flow velocity vector is set to $0$ at this boundary:

```wl
**v**	 = 	0
```

This can be set with a ``DirichletCondition``.

Specify a no-slip wall boundary:

```wl
In[92]:= DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, pred]

Out[92]= DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, pred]
```

No-slip example

This example shows the application of a no-slip wall boundary condition. A rectangular flow region with two semi-disk cutouts is considered. On the left there is an inflow profile, and on the right there is an outflow condition. The remainder of the domain is walls modeled as a no-slip boundary condition.

Set up, mesh and visualize the region:

```wl
In[26]:=
Ω = RegionDifference[Rectangle[{0, 0}, {10, 4}], RegionUnion[Disk[{5, 0}, 1], Disk[{5, 4}, 1]]];
mesh = ToElementMesh[Ω, MaxCellMeasure -> 0.02];
mesh["Wireframe"]

Out[28]= [image]
```

Set up the laminar flow equations:

```wl
In[96]:=
vars = {{u[x, y], v[x, y], p[x, y]}, {x, y}};
pars = <|"ReynoldsNumber" -> 100|>;
pde = FluidFlowPDEComponent[vars, pars] == {0, 0, 0};
```

Specify the inflow and outflow conditions:

```wl
In[99]:=
inflow = DirichletCondition[{u[x, y] == 3 / 32 (4 - y) y, v[x, y] == 0}, x == 0];
outflow = DirichletCondition[p[x, y] == 0, x == 10];
```

Specify the no-slip wall boundary conditions:

```wl
In[101]:= wall = DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, 0 < x < 10]

Out[101]= DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, 0 < x < 10]
```

Collect the boundary conditions:

```wl
In[102]:= bcs = {inflow, outflow, wall};
```

Solve the fluid flow PDE:

```wl
In[103]:= {uVel, vVel, pressure} = NDSolveValue[{pde, bcs}, vars[[1]], {x, y}∈mesh, Method -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}]

Out[103]=
{InterpolatingFunction[{{0., 10.000000000000142}, {0., 4.000000000000057}}, 
 {5, 4225, 0, {5857, 0}, {3, 0}, 0, 0, 0, 0, Indeterminate & , {}, {}, False}, 
 {NDSolve`FEM`ElementMesh[CompressedData["«114664»"], 
   {NDSolve`FEM`TriangleElement[Comp ...  5, 5, 5, 5, 13, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
      8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 
      4, 4, 4, 4, 4, 4, 4}]}]}, CompressedData["«15910»"], {Automatic}][x, y]}
```

Visualize the flow field:

```wl
In[104]:= VectorPlot[{uVel, vVel}, {x, y}∈mesh, AspectRatio -> Automatic]

Out[104]= [image]
```

Plot the $u$-velocity component at $x = 5$ :

```wl
In[105]:= Plot[uVel /. x -> 5, {y, 1, 3}]

Out[105]= [image]
```

Note that at $x = 5, y = 1$ and at $x = 5, y = 3$ the velocity is 0.

Plot the $v$-velocity component at $y = 1$ and $y = 3$ :

```wl
In[106]:= GraphicsRow[{Plot[vVel /. y -> 1, {x, 0, 10}, PlotRange -> All], Plot[vVel /. y -> 3, {x, 0, 10}, PlotRange -> All]}]

Out[106]= [image]
```

Note that at $x = 5, y = 1$ and at $x = 5, y = 3$ the velocity is 0.

Slip

In the slip, there is still no fluid going through the wall, but contrary to the no-slip condition, in the slip case there are no losses due to friction; this means that:

```wl
|              |
| ------------ |
| ϑn = 0       |
| ∂ϑt / ∂n = 0 |
```

In the axis-aligned case, the following holds for a upright wall:

```wl
u	 = 	0
```

Specify a slip wall boundary for an axis-aligned upright wall:

```wl
In[107]:= DirichletCondition[u[x, y] == 0, pred]

Out[107]= DirichletCondition[u[x, y] == 0, pred]
```

In the axis-aligned case, the following holds for a horizontal wall:

```wl
v	 = 	0
```

Specify a slip wall boundary for an axis-aligned horizontal wall:

```wl
In[108]:= DirichletCondition[v[x, y] == 0, pred]

Out[108]= DirichletCondition[v[x, y] == 0, pred]
```

Traction

Consider the Navier–Stokes momentum equation:

```wl
ρ((∂**v**/∂t) + **v**·∇**v**) + ∇·(-**τ** + p**I**) - **F**	 = 	0
```

The viscous stress tensor $**τ**$ is given as:

```wl
**τ**	 = 	2μ**Overscript[ϵ, .]** - (2/3)(∇·**v**)**I**
```

with the strain rate $Overscript[**ϵ**, .]$ :

```wl
**Overscript[ϵ, .]**	 = (1/2)(∇**v** + (∇**v**)^T)
```

Expanding all terms leads to:

```wl
ρ((∂**v**/∂t) + **v**·∇**v**) + ∇·(Underscript[-μ(∇**v** + (∇**v**)^T) + (2/3)μ(∇·**v**)**I** + p**I**, Underscript[︸,                         total stress tensor **σ**                          ]]) - **F**	 = 	0
```

where $σ = -τ + p**I**$ is the total stress. The total stress is the viscous stress $τ$ plus the pressure component.

To prescribe traction, it is necessary to specify

```wl
**n**·(τ - p**I****)** = **n**·σ = g
```

where $**g**$ is a traction value vector. With all terms substituted:

```wl
**n**·(μ**(****∇v + (∇v)^T****)** - (2/3)(∇·**v**)**I** - p**I**) = g
```

The total stress tensor $σ$ in component form is:

```wl
σ = (|                                                |                                                |                       |
| ---------------------------------------------- | ---------------------------------------------- | --------------------- |
| -2μ (∂u/∂x) - (2/3)((∂u/∂x) + (∂v/∂y) + …) + p | -μ((∂u/∂y) + (∂v/∂x))                          | -μ((∂u/∂z) + (∂w/∂x)) |
| -μ((∂u/∂y) + (∂v/∂x))                          | -2μ (∂v/∂y) - (2/3)((∂u/∂x) + (∂v/∂y) + …) + p | -μ((∂v/∂z) + (∂w/∂y)) |
| -μ((∂u/∂z) + (∂w/∂x))                          | -μ((∂v/∂z) + (∂w/∂y))                          | ⋱                     |)
```

The type of boundary conditions can be specified through ``NeumannValue``.

Traction example

The following example considers a flow with prescribed traction boundary conditions in an annulus region. On the outer edge, at $r = b$, there are the no-slip boundary conditions. These set the fluid velocity in the $u$ and $v$ directions to $0$, $**v** = **0**$.

On the inner boundary, at $r = a$, instead of defining a velocity, a traction is prescribed. The traction with the stress vector $**τ**$ is given as:

```wl
**n**·(τ - p**I****)**** = ****g** = Overscript[**r**,  ^ ]N(θ) + Overscript[**θ**,  ^ ]S(θ)
```

Here, $Overscript[**r**,  ^ ]$ is the radial unit vector given as $Overscript[**r**,  ^ ] = Overscript[**i**,  ^ ]cos(θ) + Overscript[**j**,  ^ ]sin(θ)$, and $Overscript[**θ**,  ^ ]$ is the tangent unit vector given as $Overscript[**θ**,  ^ ] = -Overscript[**i**,  ^ ]sin(θ) + Overscript[**j**,  ^ ]cos(θ)$. $Overscript[**i**,  ^ ]$ and $Overscript[**j**,  ^ ]$ are the Cartesian unit vectors in the $x$-axis and $y$-axis directions, respectively. For this example, the traction in the normal direction is set to $N(θ) = 0$ and the tangential traction to $S(θ) = cos(2θ)$. This means that there is no flow across the inner boundary but only tangentially to it.

Set up the variables and parameters:

```wl
In[4]:=
vars = {{u[x, y], v[x, y], p[x, y]}, {x, y}};
pars = <|"DynamicViscosity" -> 1, "MassDensity" -> 1, a -> 1 / 5, b -> 1|>;
```

Specify the geometry and a refined mesh:

```wl
In[6]:=
Ω = Annulus[{0, 0}, {a, b}] /. pars;
mesh = ToElementMesh[Ω, "RegionHoles" -> {0, 0}, MeshRefinementFunction -> Function[{vertices, area}, area > 0.000025 (0.1 + 20 Norm[Mean[vertices]])]];
```

Set up a laminar flow operator:

```wl
In[8]:= op = FluidFlowPDEComponent[vars, pars];
```

On the outer boundary there is the no-slip condition.

Set up a no-slip condition on the outer boundary:

```wl
In[9]:= ΓWall = DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, ElementMarker == 2]

Out[9]= DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, ElementMarker == 2]
```

On the inner surface, the tangent unit normal can be computed from the outward-pointing unit normal $**n**$ with a cross product.

Use a cross product to compute the unit tangent from the unit normal vector $**n**$:

```wl
In[10]:= Cross[{Subscript[n, x], Subscript[n, y]}]

Out[10]= {-Subscript[n, y], Subscript[n, x]}
```

To extract the first and second components of the cross product for the $u$ and $v$ directions, respectively, ``Indexed`` is used.

Set up the traction boundary conditions:

```wl
In[13]:=
Γu = NeumannValue[Indexed[Cross[NeumannBoundaryUnitNormal[{x, y}]], 1] * ((x / a) ^ 2 - (y / a) ^ 2), ElementMarker == 1];
Γv = NeumannValue[Indexed[Cross[NeumannBoundaryUnitNormal[{x, y}]], 2] * ((x / a) ^ 2 - (y / a) ^ 2), ElementMarker == 1];
```

Equivalently, the tangent unit normal can also be given as ${y / a, -x / a}$, because the unit normal vector on the inner surface is $**n** = {0, 0} - {x / a, y / a}$, and the unit tangent then is $**t** = {y / a, -x / a}$.

An equivalent set of traction boundary conditions:

```wl
In[15]:=
NeumannValue[(y / a) * ((x / a) ^ 2 - (y / a) ^ 2), ElementMarker == 1];
NeumannValue[-(x / a) * ((x / a) ^ 2 - (y / a) ^ 2), ElementMarker == 1];
```

Without specifying a boundary condition for the pressure, the pressure value will be floating, and ``NDSolve`` will give a warning that not enough boundary conditions are specified.

Solve the equations:

```wl
In[17]:= {ufun, vfun, pfun} = NDSolveValue[{op == {Γu, Γv, 0} /. pars, ΓWall}, {u, v, p}, {x, y}∈mesh, Method -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}, DependentVariables -> {u, v, p}]
```

NDSolveValue::femibcnd: No DirichletCondition or Robin-type NeumannValue was specified for {p}; the result may not be unique.

```wl
Out[17]= {InterpolatingFunction[{{-1., 1.}, {-1., 1.}}, <>], InterpolatingFunction[{{-1., 1.}, {-1., 1.}}, <>], InterpolatingFunction[«5»]}
```

Visualize the pressure solution and the velocity field:

```wl
In[18]:= Show[{DensityPlot[pfun[x, y], {x, y}∈mesh, ...], StreamPlot[{ufun[x, y], vfun[x, y]}, {x, y}∈mesh, ...]}, ...]

Out[18]= [image]
```

To verify that the traction boundary condition works, the radial and tangential components of the velocities can be computed and plotted at $r = a$ for all $θ$. For the radial velocity, a solution proportional to $-sin(2θ)$ is expected, and for the tangential velocity, a solution proportional to $-cos(2θ)$ is expected.

Define a function to compute the radial and tangential components of the velocities:

```wl
In[19]:=
uvR[x_, y_] := x / Sqrt[x ^ 2 + y ^ 2] * ufun[x, y] + y / Sqrt[x ^ 2 + y ^ 2] * vfun[x, y]
uvT[x_, y_] := -(y / Sqrt[x ^ 2 + y ^ 2]) * ufun[x, y] + x / Sqrt[x ^ 2 + y ^ 2] * vfun[x, y]
```

Plot the radial component of the velocities and a function proportional to $-sin(2θ)$ :

```wl
In[21]:= Plot[Evaluate[{uvR[a * Cos[θ], a * Sin[θ]], -Sin[2θ] * a / 20} /. pars], {θ, 0, 2π}, PlotLegends -> {"uvR", "∝ sin(2θ)"}]

Out[21]= [image]
```

Plot the tangential component of the velocities and a function proportional to $-cos(2θ)$ :

```wl
In[22]:= Plot[Evaluate[{uvT[a * Cos[θ], a * Sin[θ]], -Cos[2θ] * a / 5} /. pars], {θ, 0, 2π}, PlotLegends -> {"uvT", "∝ cos(2θ)"}]

Out[22]= [image]
```

The next step is to compute the normal stress and the shear at the inner boundary and check that they match the prescribed traction boundary conditions. The total stress tensor is given by:

```wl
σ = μ(∇v + (∇v)^T) - p**I**
```

The normal and tangential components of the traction vector can be computed by:

```wl
Subscript[F, n] = Subscript[n, i]Subscript[σ, ij]Subscript[n, j]
```

and

```wl
Subscript[F, t] = Subscript[t, i]Subscript[σ, ij]Subscript[n, j]
```

Create a function to compute the stress tensor $τ$ :

```wl
In[23]:=
gradV = Grad[{ufun[x, y], vfun[x, y]}, {x, y}];
sigma[x_, y_] := Evaluate[("DynamicViscosity" * (gradV + Transpose[gradV]) - pfun[x, y] * IdentityMatrix[2]) /. pars];
```

Compute the unit normal vector:

```wl
In[25]:= normal = ElementMeshBoundaryUnitNormal[mesh]

Out[25]= [image]
```

Define a function for the tangent vector:

```wl
In[26]:= tangent[x_, y_] := Cross[normal[x, y]]
```

Define a function to compute the normal stress:

```wl
In[27]:= normalStress[x_ ? NumberQ, y_ ? NumberQ] := normal[x, y].(sigma[x, y].normal[x, y]);
```

Define a function to compute the shear stress:

```wl
In[28]:= shearStress[x_ ? NumberQ, y_ ? NumberQ] := tangent[x, y].(sigma[x, y].normal[x, y]);
```

Visualize the computed normal and shear stress at $r = a$ :

```wl
In[29]:=
stresses = {normalStress[a * Cos[θ], a * Sin[θ]], shearStress[a * Cos[θ], a * Sin[θ]]} /. pars;
Plot[stresses, {θ, 0, 2 Pi}, ...]

Out[30]= [image]
```

Visualize the difference of the normal and shear stress at $r = a$ from the boundary conditions:

```wl
In[31]:= Plot[Evaluate[stresses - {0, Cos[2θ]}], {θ, 0, 2 Pi}, ...]

Out[31]= [image]
```

Initial Conditions

For time-dependent problems, initial conditions $Subscript[**v**, 0]$ need to be specified.

```wl
**v** = Subscript[**v**, 0](x, ...)
```

Initial conditions $Subscript[**v**, 0]$ need to satisfy the continuity equation.

If a nonzero initial velocity is specified, then the pressure field specified needs to be consistent with the velocity and vice versa. Obtaining consistent nonzero initial conditions can be tricky. And for this reason, very often, fluid flow models are started from a state of complete rest.

Convergence of Fluid Flow Models

Fluid flow PDEs can be difficult to get to converge. This section provides hints at what to try when ``NDSolve`` cannot find a solution to a fluid flow model. Combinations of the methods presented here are possible.

Initial Seeding

Typically, nonlinear solvers need an initial guess, a seed to get started. If this seed is too far off from the solution, the solver cannot converge. The radius in which the solver converges is called the radius of convergence. The closer the initial seed is to the final solution, the better.

Iterative stepping

Here the solution of a low Reynolds number flow field is used to seed the computation of a higher Reynolds number. For this, the parametric function is used.

Set up the lid-driven cavity example:

```wl
In[4]:=
vars = {{u[x, y], v[x, y], p[x, y]}, {x, y}};
pars = <|"ReynoldsNumber" -> re|>;
lid = DirichletCondition[{u[x, y] == 1, v[x, y] == 0}, y == 1];
wall = DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, y < 1];
referencePressure = DirichletCondition[p[x, y] == 0, x == 0 && y == 0];
bcs = {lid, wall, referencePressure};
pde = FluidFlowPDEComponent[vars, pars] == {0, 0, 0};
```

Set up a mesh:

```wl
In[11]:=
m1 = ToGradedMesh[Line[{{0}, {1}}], <|"Alignment" -> "BothEnds", "ElementCount" -> 75|>];
mesh = ElementMeshRegionProduct[m1, m1]

Out[12]= ElementMesh[{{0., 1.}, {0., 1.}}, {QuadElement[<5625>]}]
```

Specifying the dependent variables helps the solver to create an optimized system of equations.

Create a new parametric function, now using the constructed mesh with ``DependentVariables`` :

Compute the solution up to a Reynolds number of $30000$ in $20$ steps with an initial solution of $0$ :

```wl
In[14]:=
Clear[solution]
reMax = 30000;
nsteps = 30;
solution[0] = {0, 0, 0};
```

Iteratively solve for the high Reynolds number flow field:

```wl
In[18]:=
Monitor[Do[
	reNew = step * reMax / nsteps;
	solution[step] = pressure[reNew
	, "InitialSeeding" -> Thread[Equal[vars[[1]], solution[step - 1]]]
	];
	, {step, 1, nsteps, 1}], 
	step];
```

Vector plot of the solution at Reynolds number $30000$ :

```wl
In[19]:= VectorPlot[solution[nsteps][[{1, 2}]], {x, y}∈mesh, VectorPoints -> Tuples[Range[0.01, 0.99, 0.01], 2], Frame -> None, PlotLabel -> "Reynolds Number 30000"]

Out[19]= [image]
```

Equation Modification

This may sound absurd at first: Why would you modify the equations? Sometimes it is possible that a more stable or reduced set of equations can be used. The solution can then be used as a seed for the full equations.

Stokes equations

When no initial seed is specified, the solver assumes $0$ for all dependent variables. In the case of fluid dynamics, when an initial seed of $0$ is used for the velocity field and pressure is used, you are in effect solving Stokes-type equations.

References

*Simple Fields of Physics*; Backstrom, G.; 2006; GB Publishing; ISBN 91-975553-0-4

*Numerische Simulation in der Strömungsmechanik*; Griebel M. et al.; 1995; Vieweg; ISBN 3-528-06761-6

The Efficient Engineer; https://www.youtube.com/watch?v=VvDJyhYSJv8; retrieved 2022/12/05

"On the Use and Misuse of the Oberbeck–Boussinesq Approximation"; Barletta, A.; Celli, M; Rees, D.A.S; https://arxiv.org/abs/2202.10981

"The Validity of the Boussinesq Approximation for Liquids and Gases"; Gray, D.; Giorini, A.; *International Journal of Heat and Mass Transfer*, Volume 19, Issue 5, May 1976, pp. 545–551; https://doi.org/10.1016/0017-9310(76)90168-X

“High-Re Solutions for Incompressible Flow Using the Navier-Stokes Equations and a Multigrid Method"; Ghia, U., Ghia, K.N., Shin, C.T.; *Journal of Computational Physics*, vol. 48, pp. 387–411, 1982.; https://doi.org/10.1016/0021-9991(82)90058-4

"An Introduction to the Mechanics of Incompressible Fluids"; Deville, M. O.; 2022; Springer; ISBN 978-3-031-04682-7; https://doi.org/10.1007/978-3-031-04683-4

## Related Tech Notes

* [PDEModels Overview](https://reference.wolfram.com/language/PDEModels/tutorial/PDEModelsOverview.en.md)
* [Fluid Dynamics Model Verification Tests](https://reference.wolfram.com/language/PDEModels/tutorial/FluidDynamics/FluidDynamicsVerificationTests.en.md)