Fluid Dynamics

PDEModels Overview

Acoustics

Mass Transport

Heat Transfer

Solid Mechanics

Electromagnetics

Fluid Dynamics

Schrödinger

FluidFlowPDEComponent

Laminar Flow Monograph

FluidFlowPDEComponent

Compressible Navier-Stokes

Incompressible Navier-Stokes

Non-Newtonian Flow

Multiphysics

Navier-Stokes

Compressible Navier-Stokes

Viscous stress tensor τ

Strain rate (OverDot[ϵ])

Navier-Stokes

Compressible Navier-Stokes

Incompressible Navier-Stokes

Lid driven cavity

Ω=Rectangle[{0,0},{1,1}]

Lid driven cavity

vars={{u[x,y],v[x,y],p[x,y]},{x,y}};
pars=<|"ReynoldsNumber"->1000|>;

ℛℯ=((inertial forces)/(viscous forces))

pde=FluidFlowPDEComponent[vars,pars]=={0,0,0}

Lid driven cavity

bcs={
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]
};

Lid driven cavity

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

Lid driven cavity

StreamPlot[{uVel,vVel},{x,y}∈Ω]

InterpolationOrder

NDSolveValue[{pde,bcs},{u[x,y],v[x,y],p[x,y]},{x,y}∈Ω(*,Method->{"FiniteElement","InterpolationOrder"->{u->2,v->2,p->1}}*)];

({{ρ(u((∂u)/(∂x))+v((∂u)/(∂y)))-μ(((∂^2)u)/(∂(x^2)))-μ(((∂^2)u)/(∂(y^2)))+((∂p)/(∂x))-(Subscript[f, x]), =, 0}, {ρ(u((∂v)/(∂x))+v((∂v)/(∂y)))-μ(((∂^2)v)/(∂(x^2)))-μ(((∂^2)v)/(∂(y^2)))+((∂p)/(∂y))-(Subscript[f, y]), =, 0}, {((∂u)/(∂x))+((∂v)/(∂y)), =, 0}})

Method->{"FiniteElement","InterpolationOrder"->{u->2,v->2,p->1}};

Method->{
"PDEDiscretization"->{"MethodOfLines",
"SpatialDiscretization"->{"FiniteElement","InterpolationOrder"->{u->2,v->2,p->1}
}}};

NDSolveValue[{pde,bcs},{u[x,y],v[x,y],p[x,y]},{x,y}∈Ω,FluidFlowOptions[vars,pars,Ω]]

NDSolveValue[PDEModel[]]

Parameter specification

FluidFlowPDEComponent[vars, <|"ReynoldsNumber"->ℛℯ|>];

FluidFlowPDEComponent[vars, <|"Material"->(Entity["glycerol", Entity["Chemical","Glycerol"], "Entity[\"Chemical\", \"Glycerol\"]", "chemical"])|>];

FluidFlowPDEComponent[vars, <|"DynamicViscosity"->μ,"MassDensity"->ρ|>];

Boundary Conditions

Inflow

DirichletCondition[{u[x,y]==U,v[x,y]==V},pred]

Outflow

DirichletCondition[p[x,y]==0,pred]

Wall - no-slip

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

Wall - slip

DirichletCondition[u[x,y]==0,pred]

Boundary Conditions

({{ρ(((∂v)/(∂t))+v·∇v)+∇·(-τ+pI)-F, =, 0}})

Traction

n·(τ-pI)=n·σ=g

n·(μ(∇v+(((∇v))^T))-(2/3)(∇·v)I-pI)=g

Example in Laminar Flow monograph.

Time dependent

vars={{u[t,x,y],v[t,x,y],p[t,x,y]},t,{x,y}};

FluidFlowPDEComponent[vars,(Entity["glycerol", Entity["Chemical","Glycerol"], "Entity[\"Chemical\", \"Glycerol\"]", "chemical"])];

Time dependent

rules={length->22/10,height->41/100};
Ω=RegionDifference[Rectangle[{0,0},{length,height}],Disk[{1/5,1/5},1/20]]/.rules;
(mesh=ToElementMesh[Ω,"MaxCellMeasure"->0.0005])["Wireframe"]

Time dependent

vars={{u[t,x,y],v[t,x,y],p[t,x,y]},t,{x,y}};
pars=<|"MassDensity"->1,"DynamicViscosity"->10^-3|>;
op=FluidFlowPDEComponent[vars,pars];

Time dependent

RegionPlot[Ω,AspectRatio->Automatic]

bcs={
DirichletCondition[
{u[t,x,y]==sf[t]*6*y*(height-y)/height^2,v[t,x,y]==0},
x==0],
DirichletCondition[{u[t,x,y]==0,v[t,x,y]==0},0<x<length],
DirichletCondition[p[t,x,y]==0,x==length]
}/.rules;

ic={u[0,x,y]==0,v[0,x,y]==0,p[0,x,y]==0};

Time dependent

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]

Time dependent

AbsoluteTiming[{xVel,yVel,pressure}=NDSolveValue[{op=={0,0,0},bcs,ic},{u,v,p},{x,y}∈mesh,{t,0,10},];]

Time dependent

{minX,maxX}=MinMax[xVel["ValuesOnGrid"]];
AbsoluteTiming[frames=Table[Rasterize[ContourPlot[xVel[t,x,y],{x,y}∈mesh,PlotRange->All,AspectRatio->Automatic,ColorFunction->"TemperatureMap",Contours->Range[minX,maxX,(maxX-minX)/7],Axes->False,Frame->None],RasterSize->2*{360,68},ImageSize->{360,68}],{t,3,10,1/12}];];
ListAnimate[frames,SaveDefinitions->True]

Non-Newtonian Flow

({{ρ(((∂v)/(∂t))+v·∇v)+∇·(-τ+pI), =, F}})

Nonlinear stress-strain rate relation

Apparent viscosity

({{τ, =, 2(Subscript[μ, app])((OverDot[γ]))}})(OverDot[ϵ])

Shear rate

({{(OverDot[γ]), =, (Sqrt[2(OverDot[ϵ]):(OverDot[ϵ])])}})

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

Non-Newtonian Flow Models

Power Law

Carreou

Bingham-Papanstasiou

Herschel-Bulkley-Papanastasiou

....

Custom

Implementation code is available

Power law

({{(Subscript[μ, app]), =, ((m(((max((OverDot[γ]),(Subscript[(OverDot[γ]), min])))/(Subscript[(OverDot[γ]), ref]))))^(n-1))}})

"FluidDynamicsMaterialModel"-><|"PowerLaw"-><|"Exponent"->n,"MinimalShearRate"->(Subscript[(OverDot[γ]), min]),"ReferenceShearRate"->(Subscript[(OverDot[γ]), ref]),"PowerLawViscosity"->m|>|>;

Power law

r=RegionUnion[Rectangle[{0,-1/2},{5,1/2}],Rectangle[{5,-3/2},{30,3/2}]];
mesh=ToElementMesh[r,MaxCellMeasure->0.1];
mesh["Wireframe"]

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]
};

Power law

vars={{u[x,y],v[x,y],p[x,y]},{x,y}};
pars=<|"DynamicViscosity"->1,"MassDensity"->1|>;

parsPL=pars;
parsPL["FluidDynamicsMaterialModel"]=<|"PowerLaw"-><|"Exponent"->2,"PowerLawViscosity"->1/40|>|>;

pde=FluidFlowPDEComponent[vars,pars]=={0,0,0};
pdePL=FluidFlowPDEComponent[vars,parsPL]=={0,0,0};

Power law

uVel=NDSolveValue[{pde,bcs},u,{x,y}∈mesh,];
uVelPL=NDSolveValue[{pdePL,bcs},u,{x,y}∈mesh,];

Plot[{uVel[28,y*1.5],uVelPL[28,y*1.5]},{y,0,1},]

Energy Transport

Couple Fluid Flow with Heat Transfer

Boussinesq Approximation

Boussinesq Approximation

The density ρ is constant except in the buoyancy term

({{ρ(T), =, (Subscript[ρ, ref])(1-α(T-(Subscript[T, ref])))}})

α coefficient of thermal expansion

Energy Transport

Flow driven by heat difference

Ω=Rectangle[{0,0},{4,1}];
RegionPlot[Ω,AspectRatio->Automatic]

Ω=Rectangle[{0,0},{4,1}];
mesh=ToElementMesh[Ω,MaxCellMeasure->0.002];
mesh["Wireframe"]

varsFluid={{u[t,x,y],v[t,x,y],p[t,x,y]},t,{x,y}};
varsHeat={T[t,x,y],t,{x,y}};

Energy Transport

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

model=Flatten[{
FluidFlowPDEComponent[varsFluid,pars]-{0,boussinesq,0},HeatTransferPDEComponent[varsHeat,pars]
}];

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

Energy Transport

Ω=Rectangle[{0,0},{4,1}];
RegionPlot[Ω,AspectRatio->Automatic]

wall=DirichletCondition[{u[t,x,y]==0,v[t,x,y]==0},True];
reference=DirichletCondition[p[t,x,y]==0,x==0&&y==0];
temperatures={HeatTemperatureCondition[y==0,varsHeat,pars,<|"SurfaceTemperature"->(Subscript[T, hot])|>],HeatTemperatureCondition[y==1,varsHeat,pars,<|"SurfaceTemperature"->(Subscript[T, cold])|>]};
bcs={wall,reference,temperatures};

Energy Transport

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

ic={u[0,x,y]==0,v[0,x,y]==0,p[0,x,y]==0,T[0,x,y]==(Subscript[T, ref])};

tEnd=500;
pde={model=={0,0,0,0},bcs,ic}/.parameters;

Energy Transport

AbsoluteTiming[
{uVel,vVel,pressure,temperature}=NDSolveValue[pde,{u,v,p,T},{x,y}∈mesh,{t,0,tEnd},
];
]

Energy Transport

numberOfFrames=30;
bmesh=ToBoundaryMesh[uVel["ElementMesh"]];
frames=Table[Show[
bmesh["Wireframe"],
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]

Compressible Fluid Flow

vars={{u[x,y],v[x,y],p[x,y]},{x,y}};

FluidFlowPDEComponent[vars,<|"DynamicViscosity"->μ,"MassDensity"->ρ|>]//MatrixForm

FluidFlowPDEComponent[vars,<|"DynamicViscosity"->μ[x,y],"MassDensity"->ρ[x,y]|>]//MatrixForm

NDSolve rocks!

Summary:

2D/3D (in)compressible (non-)Newtonian Fluid Flow

ruebenko@wolfram.com

Thanks

Turbulence?

Init

Needs["NDSolve`FEM`"]
$HistoryLength=0;