Wolfram Language & System Documentation Center

FluidFlowPDEComponent

FluidFlowPDEComponent[vars,pars]

yields a flow PDE term with variables vars and parameters pars.

Details

FluidFlowPDEComponent models flow of viscous fluids subject to applied forces and constraints.
FluidFlowPDEComponent returns a sum of differential operators to be used as a part of partial differential equations:

FluidFlowPDEComponent models fluid flow phenomena with velocities , and in units of [] as dependent variables, as independent variables in units of [] and time variable in units of [].

FluidFlowPDEComponent creates PDE components for stationary, time-dependent, parametric analysis.
FluidFlowPDEComponent creates PDE components in two and three space dimensions.
Stationary variables vars are vars={{u[x₁,…,x_n],v[x₁,…,x_n],…,p[x₁,…,x_n]},{x₁,…,x_n}}.
Time-dependent or eigenmode variables vars are vars={{u[t,x₁,…,x_n],v[t,x₁,…,x_n],…,p[x₁,…,x_n]},t,{x₁,…,x_n}}.
The equations for different analysis types that FluidFlowPDEComponent generates depend on the form of vars.
FluidFlowPDEComponent creates a system of equations with the vector-valued Navier–Stokes equation combined with the continuity equation.
The time-dependent equilibrium equation of the fluid dynamics PDE FluidFlowPDEComponent with mass density [ $TemplateBox[{InterpretationBox[, 1], {"kg", , "/", , {"m", ^, 3}}, kilograms per meter cubed, {{(, "Kilograms", )}, /, {(, {"Meters", ^, 3}, )}}}, QuantityTF]$ ], fluid velocity [], time [], viscous stress tensor [], pressure [], identity matrix and body load vectors [ $TemplateBox[{InterpretationBox[, 1], {"N", , "/", , {"m", ^, 3}}, newtons per meter cubed, {{(, "Newtons", )}, /, {(, {"Meters", ^, 3}, )}}}, QuantityTF]$ ] is based on the Navier–Stokes equation and the continuity equation:

In compressible form, the viscous stress tensor is given as:

Here [] is the dynamic viscosity, and infinitesimal, small deformation strain rate measure [1/] is given as:

FluidFlowPDEComponent creates a PDE model for compressible or incompressible fluid flow, depending on the nature of the values of the mass density .
The compressible fluid dynamics model is given as:

For constant values of the mass density , the mass continuity equation simplifies to a volumetric continuity equation , and with that, the viscous stress tensor simplifies to .
The incompressible fluid dynamics model is given as:

The stationary equilibrium equations have .
The units of the Navier–Stokes model terms are a force density in [ $TemplateBox[{InterpretationBox[, 1], {"N", , "/", , {"m", ^, 3}}, newtons per meter cubed, {{(, "Newtons", )}, /, {(, {"Meters", ^, 3}, )}}}, QuantityTF]$ ].
The units of the mass continuity equation model terms are a force density in [ $TemplateBox[{InterpretationBox[, 1], {"kg", , "/(", , {"m", ^, 3}, , "s", , ")"}, kilograms per meter cubed second, {{(, "Kilograms", )}, /, {(, {{"Meters", ^, 3}, , "Seconds"}, )}}}, QuantityTF]$ ], and for the volumetric continuity [1/].
Laminar flow is typical for , where is the Reynolds number.
The Reynolds number is defined as , where [] is a characteristic length and the flow velocity.
The following parameters pars can be given:

parameter	default	symbol
"CoordinateChart"	-	none
"DynamicViscosity"	-	, dynamic viscosity []
"FluidLoad"	0	, body force density [ $TemplateBox[{InterpretationBox[, 1], {"N", , "/", , {"m", ^, 3}}, newtons per meter cubed, {{(, "Newtons", )}, /, {(, {"Meters", ^, 3}, )}}}, QuantityTF]$ ]
"FluidDynamicsMaterialModel"	"Newtonian"	none
"MassDensity"	-	, density in [ $TemplateBox[{InterpretationBox[, 1], {"kg", , "/", , {"m", ^, 3}}, kilograms per meter cubed, {{(, "Kilograms", )}, /, {(, {"Meters", ^, 3}, )}}}, QuantityTF]$ ]
"Material"	-	none
"ModelForm"	"Conservative"	none
"RegionSymmetry"	None	-
"ReynoldsNumber"	-

If a "Material" is specified, the material constants are extracted from the material data; otherwise, relevant material parameters need to be specified.
Instead of material parameters, a Reynolds number can be specified:

"RegionSymmetry"->"Axisymmetric" uses a truncated cylindrical coordinate system and assumes in the direction and pressure and velocity components and independent of .
The axis of symmetry is the axis.
A "CoordinateChart" uses the same specification as CoordinateChartData.
"CoordinateChart"->"Cylindrical" with pars assumes and can be used to model swirling flow. »
The default material model is a Newtonian flow model.
Alternate material models can be specified by setting the "FluidDynamicsMaterialModel" key in parameters pars.
The following non-Newtonian material models are available:
material model name

"PowerLaw"

"Carreau"

"Bingham-Papanastasiou"

"Herschel-Bulkley-Papanastasiou"

"Casson-Papanastasiou"

The apparent viscosity is a function of the shear rate .
Additional material model-specific parameters for a model with "ModelName" can be specified through "FluidDynamicsMaterialModel"-><|"ModelName"-><|...|>|>.
The "PowerLaw" model, a general-purpose model, implements .
The following parameters can be given for the "PowerLaw" model:

parameter	default	symbol
"PowerLawExponent"		, exponent
"MinimalShearRate"		, minimal shear rate
"ReferenceShearRate"		, reference shear rate
"PowerLawViscosity"		, power law viscosity

The generalized "Carreau" model, useful for polymer or blood flow, implements .
The following parameters can be given for the "Carreau" model:

parameter	default	symbol
"PowerLawExponent"		, exponent
"TransitionExponent"	2	, exponent
"InfiniteShearRateViscosity"		, viscosity at inifinite shear rate
"Lambda"		, relaxation time []
"ZeroShearRateViscosity"		, viscosity at zero shear rate

The "Carreau" model can also be used for a Cross model where with .
The "Bingham-Papanastasiou" model, useful for viscoplastic material, implements .
The following parameters can be given for the "Bingham-Papanastasiou" model:

parameter	default	symbol
"PlasticViscosity"		, plastic viscosity
"YieldStress"		, yield stress
"ShearRateFactor"		, shear rate factor

The "Herschel-Bulkley-Papanastasiou" model, a mixture of a "PowerLaw" and "Bingham-Papanastasiou" model, implements . The model uses a power law to compute the plastic viscosity of the Bingham–Papanastasiou model, and parameters for both models can be set.
The "Casson-Papanastasiou" model, useful for viscoplastic material or blood flow, implements .
The following parameters can be given for the "Casson-Papanastasiou" model:

parameter	default	symbol
"PlasticViscosity"		, plastic viscosity
"YieldStress"		, yield stress
"ShearRateFactor"		, shear rate factor

A custom apparent viscosity function fun can be specified as "FluidDynamicsMaterialModel"<|"Custom"<|"ApparentViscosityFunction"fun|>|>
A custom apparent viscosity function fun has a function signature fun[vars_,pars_,data__].
A custom viscous stress tensor function fun can be specified as "FluidDynamicsMaterialModelFunction"fun.
A custom viscous stress tensor function fun has a function signature fun[vars_,pars_,data__].
Non-isothermal flow can be modeled through the Boussinesq approximation.
FluidFlowPDEComponent uses "SIBase" units. The geometry has to be in the same units as the PDE.
If the FluidFlowPDEComponent depends on parameters that are specified in the association pars as …,keyp_i…,p_iv_i,…, the parameters are replaced with .

Examples

open all close all

Basic Examples (4)

Define a flow PDE model:

Wolfram Language code: FluidFlowPDEComponent[{{u[x, y], v[x, y], p[x, y]}, {x, y}}, Entity["Chemical", "Glycerol"]]//MatrixForm

Define a symbolic flow PDE:

Wolfram Language code:

FluidFlowPDEComponent[{{u[x, y], v[x, y], p[x, y]}, {x, y}}, <|"DynamicViscosity" -> μ, "MassDensity" -> ρ|>]//MatrixForm

Define a stationary flow PDE model with a Reynolds number of 100:

Wolfram Language code: FluidFlowPDEComponent[{{u[x, y], v[x, y], p[x, y]}, {x, y}}, <|"ReynoldsNumber" -> 100|>]//MatrixForm

Solve for the flow velocity and pressure in a driven cavity with a Reynolds number of 1000:

Wolfram Language code:

{uVel, vVel, pressure} = NDSolveValue[{FluidFlowPDEComponent[{{u[x, y], v[x, y], p[x, y]}, {x, y}}, <|"ReynoldsNumber" -> 1000|>] == {0, 0, 0}, DirichletCondition[{u[x, y] == 1, v[x, y] == 0}, y == 1], DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, y < 1], DirichletCondition[p[x, y] == 0, x == 0 && y == 0]}, {u[x, y], v[x, y], p[x, y]}, {x, y}∈Rectangle[{0, 0}, {1, 1}], Method -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}];

Visualize the velocity of the fluid:

Wolfram Language code: StreamPlot[{uVel, vVel}, {x, y}∈Rectangle[{0, 0}, {1, 1}]]

Scope (13)

2D (4)

Define a flow PDE model for a specific material:

Wolfram Language code:

FluidFlowPDEComponent[{{u[x, y], v[x, y], p[x, y]}, {x, y}}, <|"Material" -> Entity["Chemical", "Glycerol"]|>]//MatrixForm

Specify a flow PDE with a dynamic viscosity of and a mass density of :

Wolfram Language code:

FluidFlowPDEComponent[{{u[x, y], v[x, y], p[x, y]}, {x, y}}, <|"DynamicViscosity" -> Quantity[0.934, "Pascals"*"Seconds"], "MassDensity" -> Quantity[1.25, "Grams"/"Centimeters"^3]|>]//MatrixForm

Activate a flow PDE model for a specific material:

Wolfram Language code: Activate[FluidFlowPDEComponent[{{u[x, y], v[x, y], p[x, y]}, {x, y}}, <|"Material" -> Entity["Chemical", "Glycerol"]|>]]

Specify a symbolic stationary fluid dynamics PDE in two dimensions with a dynamic viscosity of and a mass density of :

Wolfram Language code: FluidFlowPDEComponent[{{u[x, y], v[x, y], p[x, y]}, {x, y}}, <|"DynamicViscosity" -> μ, "MassDensity" -> ρ|>]

2D Axisymmetric (1)

Specify a symbolic stationary axisymmetric fluid dynamics PDE with a dynamic viscosity of and a mass density of :

Wolfram Language code:

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

3D (1)

Specify a symbolic stationary fluid dynamics PDE in three dimensions with a dynamic viscosity of and a mass density of :

Wolfram Language code:

FluidFlowPDEComponent[{{u[x, y, z], v[x, y, z], w[x, y, z], p[x, y, z]}, {x, y, z}}, <|"DynamicViscosity" -> μ, "MassDensity" -> ρ|>]

Time Dependent (3)

Define a time-dependent flow PDE model:

Wolfram Language code: FluidFlowPDEComponent[{{u[t, x, y], v[t, x, y], p[t, x, y]}, t, {x, y}}, Entity["Chemical", "Glycerol"]]//MatrixForm

Specify a symbolic time-dependent fluid dynamics PDE in two dimensions with a dynamic viscosity of and a mass density of :

Wolfram Language code:

FluidFlowPDEComponent[{{u[t, x, y], v[t, x, y], p[t, x, y]}, t, {x, y}}, <|"DynamicViscosity" -> μ, "MassDensity" -> ρ|>]

Specify a symbolic time-dependent axisymmetric fluid dynamics PDE with a dynamic viscosity of and a mass density of :

Wolfram Language code:

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

Compressible Flow (1)

A compressible fluid dynamics PDE is generated if the mass density is a function of space or time :

Wolfram Language code: FluidFlowPDEComponent[{{u[x, y], v[x, y], p[x, y]}, {x, y}}, <|"DynamicViscosity" -> μ, "MassDensity" -> ρ[x, y]|>]

An incompressible fluid dynamics PDE is generated if the mass density is a constant:

Wolfram Language code: FluidFlowPDEComponent[{{u[x, y], v[x, y], p[x, y]}, {x, y}}, <|"DynamicViscosity" -> μ, "MassDensity" -> ρ|>]

Non-Newtonian Flow (3)

Specify a symbolic PDE for a non-Newtonian power-law model:

Wolfram Language code:

FluidFlowPDEComponent[{{u[x, y], v[x, y], p[x, y]}, {x, y}}, <|"MassDensity" -> ρ, "FluidDynamicsMaterialModel" -> <|"PowerLaw" -> <|"Exponent" -> n, "MinimalShearRate" -> Subscript[Overscript[γ, .], min], "ReferenceShearRate" -> Subscript[Overscript[γ, .], ref], "PowerLawViscosity" -> m|>|>|>]//MatrixForm

Custom fluid models can be obtained by specifying the apparent viscosity. Define the custom apparent viscosity function:

Wolfram Language code:

apparentViscosity[vars_, pars_, data_] := Module[{strainRate, shearRate, viscosity, exponent}, viscosity = pars["FluidDynamicsMaterialModel"]["Custom"]["viscosity"];
	exponent = pars["FluidDynamicsMaterialModel"]["Custom"]["exponent"];
	strainRate = "StrainRate" /. data;
	shearRate = Sqrt[2 * ArrayDot[strainRate, strainRate, 2]];
	viscosity * (shearRate) ^ (exponent - 1)]

Specify a symbolic non-Newtonian fluid dynamics PDE with the custom apparent viscosity:

Wolfram Language code:

FluidFlowPDEComponent[{{u[x, y], v[x, y], p[x, y]}, {x, y}}, <|"MassDensity" -> ρ, "FluidDynamicsMaterialModel" -> <|"Custom" -> <|"ApparentViscosityFunction" -> apparentViscosity, "viscosity" -> m, "exponent" -> n|>|>|>]//MatrixForm

Custom fluid models can be specified by specifying the viscous stress tensor. Define the custom viscous stress tensor function:

Wolfram Language code:

viscousStress[vars_, pars_, data_] := Module[{strainRate, shearRate, viscosity, exponent, apparentViscosity}, viscosity = pars["viscosity"];
	exponent = pars["exponent"];
	strainRate = "StrainRate" /. data;
	shearRate = Sqrt[2 * ArrayDot[strainRate, strainRate, 2]];
	apparentViscosity = viscosity * (shearRate) ^ (exponent - 1);
	2 * apparentViscosity * strainRate]

Specify a symbolic non-Newtonian fluid dynamics PDE with the custom viscous stress tensor:

Wolfram Language code:

FluidFlowPDEComponent[{{u[x, y], v[x, y], p[x, y]}, {x, y}}, <|"MassDensity" -> 1, "FluidDynamicsMaterialModelFunction" -> viscousStress, "viscosity" -> m, "exponent" -> n|>]//MatrixForm

Applications (8)

Stationary Analysis (1)

Solve for the velocity and pressure in a cavity, where the flow is driven at the top of a box:

Wolfram Language code:

{uVel, vVel, pressure} = NDSolveValue[{FluidFlowPDEComponent[{{u[x, y], v[x, y], p[x, y]}, {x, y}}, <|"ReynoldsNumber" -> 1000|>] == {0, 0, 0}, DirichletCondition[{u[x, y] == 1, v[x, y] == 0}, y == 1], DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, y < 1], DirichletCondition[p[x, y] == 0, x == 0 && y == 0]}, {u[x, y], v[x, y], p[x, y]}, {x, y}∈Rectangle[{0, 0}, {1, 1}], Method -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}];

Visualize the velocity of the fluid:

Wolfram Language code: StreamPlot[{uVel, vVel}, {x, y}∈Rectangle[{0, 0}, {1, 1}]]

Non-Newtonian Flow (1)

Compute the fluid flow of a non-Newtonian fluid in an opening channel by using the power law fluid flow model.

Set up the region:

Wolfram Language code: Ω = RegionUnion[Rectangle[{0, -1 / 2}, {5, 1 / 2}], Rectangle[{5, -3 / 2}, {30, 3 / 2}]];

Set up variables, flow parameters and the non-Newtonian power law for parameters for the exponent and power law viscosity:

Wolfram Language code:

vars = {{u[x, y], v[x, y], p[x, y]}, {x, y}};
pars = <|"MassDensity" -> 1, 
	"FluidDynamicsMaterialModel" -> <|"PowerLaw" -> <|"PowerLawExponent" -> 1 / 2, "PowerLawViscosity" -> 0.008838835|>|>
	|>;

Solve the PDE with an inflow profile given as {1/2,0} and with the outflow pressure set to 0. The walls are set up as no-slip walls:

Wolfram Language code:

{uVelocity, vVelocity, pressure} = NDSolveValue[{FluidFlowPDEComponent[vars, pars] == {0, 0, 0}, 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]}, {u, v, p}, {x, y}∈Ω, Method -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}]

Visualize the velocity in the region:

Wolfram Language code: ContourPlot[Sqrt[uVelocity[x, y] ^ 2 + vVelocity[x, y] ^ 2], {x, y}∈Ω, ...]

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

Wolfram Language code: Plot[uVelocity[28, y * 1.5], {y, 0, 1}]

Time-Dependent Analysis (1)

Solve a time-dependent driven cavity problem.

Create and visualize an auxiliary function to ramp up the flow speed on top of the box:

Wolfram Language code:

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]

Set up the PDE:

Wolfram Language code:

pde = FluidFlowPDEComponent[{{u[t, x, y], v[t, x, y], p[t, x, y]}, t, {x, y}}, <|"ReynoldsNumber" -> 1000|>] == {0, 0, 0};

Set up the boundary conditions:

Wolfram Language code:

bcs = {DirichletCondition[{u[t, x, y] == ramp[t], v[t, x, y] == 0}, y == 1], DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, y < 1], DirichletCondition[p[t, x, y] == 0, x == 0 && y == 0]};

Set up the initial conditions:

Wolfram Language code: ics = {u[0, x, y] == 0, v[0, x, y] == 0, p[0, x, y] == 0};

Monitor the solution process and measure the time it takes to solve the PDE:

Wolfram Language code:

AbsoluteTiming[Monitor[{uVelocity, vVelocity, pressure} = NDSolveValue[{pde, bcs, ics}, {u, v, p}, {x, y}∈Rectangle[{0, 0}, {1, 1}], {t, 0, 20}, 
	Method -> {
	"PDEDiscretization" -> {"MethodOfLines", 
	"SpatialDiscretization" -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}
	}}}, EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])];, monitor]]

Visualize rasterized frames of the velocity field for various times:

Wolfram Language code:

frames = StreamPlot[{uVelocity[#, x, y], vVelocity[#, x, y]}, {x, y}∈Rectangle[{0, 0}, {1, 1}], PlotLabel -> "Time: " <> ToString[#] <> " s", Frame -> False]& /@ Range[1, 20, 1];
frames = Rasterize[#, "Image", ImageResolution -> 120]& /@ frames;
ListAnimate[frames]

Axisymmetric Flow (2)

Compute the fluid flow in a narrowing cylindrical pipe driven by a parabolic inflow. Set up the axisymmetric domain:

Wolfram Language code:

domain = RegionUnion[RegionDifference[RegionDifference[RegionUnion[Rectangle[{0, 0}, {2, 7}], Rectangle[{0, 4}, {1, 10}]], Disk[{3 / 2, 7}, 1 / 2]], Disk[{2, 13 / 2}, 1 / 2]], Disk[{3 / 2, 6}, 1 / 2]]

Set up the PDE:

Wolfram Language code:

pde = FluidFlowPDEComponent[{{u[r, z], 0, w[r, z], p[r, z]}, {r, θ, z}}, <|"MassDensity" -> 1, "DynamicViscosity" -> 1, "RegionSymmetry" -> "Axisymmetric"|>] == {0, 0, 0};

Set up the boundary conditions:

Wolfram Language code:

bcs = {DirichletCondition[{u[r, z] == 0, w[r, z] == 0.25 * (4 - r ^ 2)}, z == 0], 
	DirichletCondition[{u[r, z] == 0, w[r, z] == 0}, r > 0 && z != 0 && z != 10], 
	DirichletCondition[u[r, z] == 0, r == 0], 
	DirichletCondition[u[r, z] == 0, z == 10], 
	DirichletCondition[p[r, z] == 0, r == 0 && z == 10]};

Solve the equations:

Wolfram Language code:

{rVel, zVel, pressure} = NDSolveValue[{pde, 
	bcs}, {u[r, z], w[r, z], p[r, z]}, {r, z}∈domain, 
	Method -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, w -> 2, p -> 1}}];

Visualize the velocity of the solution:

Wolfram Language code: VectorPlot[{rVel, zVel}, {r, z}∈domain]

Plot the pressure in 3D for part of the pipe:

Wolfram Language code:

rmf = RegionMember[domain];
Legended[RegionPlot3D[rmf[{Sqrt[x ^ 2 + y ^ 2], z}] && 0 <= PlanarAngle[{0, 0} -> {{0.02, 0}, {x, y}}] <= (4 π) / 3, {x, -2, 2}, {y, -2, 2}, {z, 0, 10}, ...], BarLegend[...]]

Compute the velocity profile for time-dependent axisymmetric flow driven by a pressure gradient. Set up the domain:

Wolfram Language code:

R = 3 / 2;L = 10;
domain = RegionDifference[Rectangle[{0, 0}, {R, L}], Disk[{R, L / 2}, R / 2]]

Set up the PDE:

Wolfram Language code:

pde = FluidFlowPDEComponent[{{u[t, r, z], 0, w[t, r, z], p[t, r, z]}, t, {r, θ, z}}, <|"DynamicViscosity" -> 1 / 10, "MassDensity" -> 1, "RegionSymmetry" -> "Axisymmetric"|>] == {0, 0, 0};

Set up the boundary conditions with a prescribed pressure gradient:

Wolfram Language code:

bcs = {DirichletCondition[u[t, r, z] == 0, z == 0], 
	DirichletCondition[u[t, r, z] == 0, r == 0], 
	DirichletCondition[{u[t, r, z] == 0, w[t, r, z] == 0}, r > 0 && z != 0 && z != L], 
	DirichletCondition[p[t, r, z] == Min[t, 1] * 3 / 2 * L, z == 0], 
	DirichletCondition[{p[t, r, z] == 0, u[t, r, z] == 0}, z == L]};

Set up the initial condition:

Wolfram Language code: ic = {u[0, r, z] == 0, w[0, r, z] == 0, p[0, r, z] == 0};

Solve the equations and measure the time it takes:

Wolfram Language code:

tEnd = 5;
AbsoluteTiming[{rVel, zVel, pressure} = NDSolveValue[{pde, bcs, ic}, {u, w, p}, {r, z}∈domain, {t, 0, tEnd}, Method -> {"PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, w -> 2, p -> 1}}}}];]

Visualize the rasterized velocity for various times:

Wolfram Language code:

frames = VectorPlot[{rVel[#, r, z], zVel[#, r, z]}, {r, z}∈domain, ...]& /@ Range[1, tEnd, 1];
frames = Rasterize[#, "Image", ImageResolution -> 120]& /@ frames;
ListAnimate[frames]

Plot the outflow velocity profile for various times:

Wolfram Language code: Plot[Evaluate[zVel[#, r, L]& /@ Range[1, tEnd, 1]], {r, 0, R}, ...]

Swirling Flow (3)

When a wall in an axisymmetric fluid flow rotates, it produces a swirling flow. Model a swirl flow in a tube where the inner wall rotates.

Set up a region:

Wolfram Language code:

ri = 1 / 2;ro = 3 / 2;h = 3;
Ω = Rectangle[{ri, 0}, {ro, h}]

Create the variables and flow parameters:

Wolfram Language code:

vars = {{u[r, z], v[r, z], w[r, z], p[r, z]}, {r, θ, z}};
pars = <|"DynamicViscosity" -> 3, "MassDensity" -> 2, "CoordinateChart" -> "Cylindrical"|>;

Set up the swirl flow operator:

Wolfram Language code: op = Activate[FluidFlowPDEComponent[vars, pars]];

At the inner wall, a nonzero angular velocity of is set up:

Wolfram Language code: stirrer = {DirichletCondition[{u[r, z] == 0, v[r, z] == r * 5π, w[r, z] == 0}, r == ri]};

The top and bottom have no velocity in the direction:

Wolfram Language code: topbottom = DirichletCondition[w[r, z] == 0, z == 0 || z == h];

The outer wall has a no-slip boundary condition:

Wolfram Language code: wall = DirichletCondition[{u[r, z] == 0, v[r, z] == 0, w[r, z] == 0}, r == ro];

A reference pressure is set up at the lower-left corner of the geometry:

Wolfram Language code: refPressure = DirichletCondition[p[r, z] == 0, r == ri && z == 0];

Solve the equation:

Wolfram Language code:

{ufun, vfun, wfun, pfun} = NDSolveValue[{op == {0, 0, 0, 0}, stirrer, topbottom, wall, refPressure}, {u, v, w, p}, {r, z}∈Ω, Method -> {"PDEDiscretization" -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, w -> 2, p -> 1}}}];

Construct the analytical solution of the azimuthal velocity :

Wolfram Language code:

analyticalSolution = With[{ω = 5π}, {a = -ω * ri ^ 2 / (ro ^ 2 - ri ^ 2), b = ω * ri ^ 2 * ro ^ 2 / (ro ^ 2 - ri ^ 2)}, 
	Function[r, a * r + b / r]]

Compare the analytical and the numerical solution:

Wolfram Language code: Plot[{analyticalSolution[r] - vfun[r, h / 2]}, {r, ri, ro}, PlotRange -> All]

Construct the velocity magnitude:

Wolfram Language code: velMag = FunctionInterpolation[Sqrt[ufun[r, z] ^ 2 + vfun[r, z] ^ 2 + wfun[r, z] ^ 2], {r, z}∈Ω];

Visualize the velocity magnitude:

Wolfram Language code: Show[DensityPlot[velMag[r, z], {r, z}∈Ω, ...], ContourPlot[vfun[r, z], {r, z}∈Ω, ...]]

Model a three-blade mixer in a beaker. Define the geometry:

Wolfram Language code:

Ω = 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}}];
RegionPlot[Ω]

Create the fluid flow operator in cylindrical coordinates:

Wolfram Language code:

op = Activate[FluidFlowPDEComponent[{{u[r, z], v[r, z], w[r, z], p[r, z]}, {r, θ, z}}, <|"DynamicViscosity" -> 10 ^ -3, "MassDensity" -> 10 ^ 3, "CoordinateChart" -> "Cylindrical"|>]];

Set a rotation speed at the shaft:

Wolfram Language code:

rotatingShaft = DirichletCondition[{u[r, z] == 0, v[r, z] == r * π, w[r, z] == 0}, 0.015 ≤ z < 0.05 && 0.001 ≤ r < 0.02]

Set wall boundary conditions:

Wolfram Language code: wall = DirichletCondition[{u[r, z] == 0, v[r, z] == 0, w[r, z] == 0}, r == 0.02 || z == 0];

On the axis of symmetry, the and velocities are set to 0:

Wolfram Language code: 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 direction:

Wolfram Language code: top = DirichletCondition[w[r, z] == 0, z == 0.05];

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

Wolfram Language code: refPressure = DirichletCondition[p[r, z] == 0, r == 0.02 && z == 0.05];

Solve the equations:

Wolfram Language code:

{ufun, vfun, wfun} = NDSolveValue[{op == {0, 0, 0, 0}, rotatingShaft, wall, symmetryAxis, top, refPressure}, {u, v, w}, {r, z}∈Ω, Method -> {"PDEDiscretization" -> {"FiniteElement"
	, "InterpolationOrder" -> {u -> 2, v -> 2, w -> 2, p -> 1}
	}}];

Compute the velocity magnitude:

Wolfram Language code: velMag = FunctionInterpolation[Sqrt[ufun[r, z] ^ 2 + vfun[r, z] ^ 2 + wfun[r, z] ^ 2], {r, z}∈Ω];

Visualize the result:

Wolfram Language code:

Show[HighlightMesh[Polygon[{{0.01, 0.017}, {0.01, 0.015}, {0.001, 0.015}, {0., 0.015}, {0., 0.017}, {0., 0.019}, 
  {0., 0.020999999999999998}, {0., 0.023}, {0., 0.025}, {0., 0.05}, {0.001, 0.05}, {0.001, 0.025}, 
  {0.01, 0.025}, {0.01, 0.023}, {0.001, 0.023}, {0.001, 0.020999999999999998}, 
  {0.006666666666666666, 0.020999999999999998}, {0.006666666666666666, 0.019}, {0.001, 0.019}, 
  {0.001, 0.017}}], {...}], 
	DensityPlot[velMag[r, z], {r, z}∈Ω, ...], ContourPlot[vfun[r, z], {r, z}∈Ω, ...], Rule[...]]

Simulate a Taylor–Couette flow with various rotating speeds of the inner cylinder. Set up a region:

Wolfram Language code: Ω = Rectangle[{1, -1.5}, {1.5, 1.5}];

Setup a swirling flow PDE with the rotational speed as parameter:

Wolfram Language code:

parafun = ParametricNDSolveValue[{Activate[FluidFlowPDEComponent[{{u[r, z], v[r, z], w[r, z], p[r, z]}, {r, θ, z}}, <|"DynamicViscosity" -> 1, "MassDensity" -> 1, "CoordinateChart" -> "Cylindrical"|>]] == {0, 0, 0, 0}, 
	DirichletCondition[{u[r, z] == 0, v[r, z] == r * ω, w[r, z] == 0}, r == 1], DirichletCondition[{u[r, z] == 0, v[r, z] == 0, w[r, z] == 0}, r > 1], DirichletCondition[p[r, z] == 0, r == 1 && z == -1.5]}, {u, v, w, p}, {r, z}∈Ω, ω, Method -> {"PDEDiscretization" -> {"FiniteElement"
	, "InterpolationOrder" -> {u -> 2, v -> 2, w -> 2, p -> 1}}}]

Compute the solution for :

Wolfram Language code: {ufun, vfun, wfun, pfun} = parafun[50];

Verify that the numerical solution matches the analytical solution for this rotational speed:

Wolfram Language code:

Plot[Evaluate[vfun[r, 0] - (-ω * ri ^ 2 / (ro ^ 2 - ri ^ 2) * r + ω * ro ^ 2 / (ro ^ 2 - ri ^ 2) / r /. {ri -> 1, ro -> 3 / 2, ω -> 50})], {r, 1, 1.5}]

Visualize the velocity magnitude:

Wolfram Language code:

velMag = FunctionInterpolation[Sqrt[ufun[r, z] ^ 2 + vfun[r, z] ^ 2 + wfun[r, z] ^ 2], {r, z}∈Ω];
Show[DensityPlot[velMag[r, z], {r, z}∈Ω, ...], ContourPlot[vfun[r, z], {r, z}∈Ω, ...], ...]

Compute and visualize the velocity magnitude for a rotational speed of :

Wolfram Language code:

{ufun, vfun, wfun, pfun} = parafun[150];
velMag = FunctionInterpolation[Sqrt[ufun[r, z] ^ 2 + vfun[r, z] ^ 2 + wfun[r, z] ^ 2], {r, z}∈Ω];
Show[DensityPlot[velMag[r, z], {r, z}∈Ω, ...], ...]

Note that the flow starts to show toroidal vortices. To increase the rotational velocity to , an initial seeding with the previous solution is used to solve the now highly nonlinear equations. Compute and visualize the velocity magnitude for a rotational speed of :

Wolfram Language code:

{ufun, vfun, wfun, pfun} = parafun[250, InitialSeeding -> {u[r, z] == ufun[r, z], v[r, z] == vfun[r, z], w[r, z] == wfun[r, z], p[r, z] == pfun[r, z]}];
velMag = FunctionInterpolation[Sqrt[ufun[r, z] ^ 2 + vfun[r, z] ^ 2 + wfun[r, z] ^ 2], {r, z}∈Ω];
Show[
	DensityPlot[velMag[r, z], {r, z}∈Ω, ...], ContourPlot[vfun[r, z], {r, z}∈Ω, ...], ...]

Top

More Learning

Tech Support

Wolfram Solutions

Wolfram Solutions For Education

Get Started

Grow Your Skills

Work with Us

Educational Programs for Adults

Educational Programs for Youth

Read

FluidFlowPDEComponent

Details

Examples

Basic Examples (4)

Scope (13)

2D (4)

2D Axisymmetric (1)

3D (1)

Time Dependent (3)

Compressible Flow (1)

Non-Newtonian Flow (3)

Applications (8)

Stationary Analysis (1)

Non-Newtonian Flow (1)

Time-Dependent Analysis (1)

Axisymmetric Flow (2)

Swirling Flow (3)

Text

CMS

APA

BibTeX

BibLaTeX

FluidFlowPDEComponent

Details

Examples

Basic Examples (4)

Scope (13)

2D (4)

2D Axisymmetric (1)

3D (1)

Time Dependent (3)

Compressible Flow (1)

Non-Newtonian Flow (3)

Applications (8)

Stationary Analysis (1)

Non-Newtonian Flow (1)

Time-Dependent Analysis (1)

Axisymmetric Flow (2)

Swirling Flow (3)

See Also

Tech Notes

Related Guides

History

Text

CMS

APA

BibTeX

BibLaTeX