# Spalart-Allmaras turbulence model

In this post Alex gives an implementation of the Spalart-Allmaras turbulence model [1, 2]. The example produces reasonable results, as far as I can tell. However, the implementation Alex uses deviates from the paper and standard references in several ways. I would like to write a SA turbulence model that is closer to the original definition.

So we have:

saEQN = D[nut[t, x, y], t] + u[t, x, y]*D[nut[t, x, y], x] ==
cb1*(1 - ft2)*STilda*
nut[t, x,
y] - (cw1 fw - (cb1/\[Kappa]^2)* ft2)*(nut[t, x, y]/d)^2 + (1/
sigma (D[(nu + nut[t, x, y]) D[nut[t, x, y], x], x] +
D[(nu + nut[t, x, y]) D[nut[t, x, y], y], y]) +
cb2/sigma*(D[nut[t, x, y], x]^2 + D[nut[t, x, y], y]^2));


Couple this to the fluid equation:

mut = rho*nut[t, x, y]*fv1;
fluidEQN =
D[u[t, x, y], t] + u[t, x, y]*D[u[t, x, y], x] + px ==
D[(mu + mut) D[u[t, x, y], y], y] + D[u[t, x, y], x, x];


Here are the model definitions:

sigma = 2/3; kappa = .41; cb1 = 0.1355; cb2 = 0.622; eps = 10^-6; d = Sqrt[y^2 + eps^2]; cw1 = cb1/kappa^2 + (1 + cb2)/sigma;
cw2 = 0.3; cw3 = 2; cv1 = 7.1; ct1 = 1; ct2 = 2; ct3 = 1.2; ct4 = 0.5;
mu = 1.711*10^-5;
rho = 1;
nu = mu/rho;

omega = 1/2 Sqrt[(D[u[t, x, y], y] - D[u[t, x, y], x])^2];
chi = nut[t, x, y]/nu;
fv1 = chi^3/(chi^3 + cv1^3);
fv2 = 1 - chi/(1 + chi*fv1);
Omega = Sqrt[2*omega*omega];
STilda = Omega + nut[t, x, y]/(kappa^2*d^2)*fv2;
g = r + cw2 (r^6 - r);
r = Min[nut[t, x, y]/(STilda*kappa^2*d^2), 10];
fw = g*((1 + cw3^6)/(g^6 + cw3^6))^(1/6);
ft2 = ct3 Exp[-ct4*chi^2];


Set up the coupled equations:

eq = {saEQN, fluidEQN};


Set up the boundary conditions. These are updated boundary conditions Alex provided to avoid a message of NDSolve.

bcUpdate = {
nut[t, x, 0] == 0, u[t, x, 0] == 0,
nut[t, x, L] == 0.1, u[t, x, L] == 1,
nut[t, 0, y] == 0.1 y/L, u[t, 0, y] == y/L,
Derivative[0, 1, 0][nut][t, L, y] == 0,
Derivative[0, 1, 0][u][t, L, y] == 0};


A side question: My understanding is that for the SA turbulence model typically nut is set to 0 at the boundaries, why is that not the case here?

Initial conditions:

ic = {nut[0, x, y] == 0.1 y/L, u[0, x, y] == y/L};


Geometry, time and pressure:

L = 10^4;
t0 = 15;
px = 0;


Solve the equations:

Monitor[{nutTGP, uTGP} =
NDSolveValue[{eq, ic, bcUpdate}, {nut, u}, {t, 0, t0}, {x, 0,
L}, {y, 0, L},
EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])
];, monitor]


NDSolve give a warning message:

General::munfl: Exp[-87139.3] is too small to represent as a normalized machine number; precision may be lost.

I think we can ignore it. When I look at the result

Plot3D[nutTGP[t0, x, y], {x, 0, L}, {y, 0, L}, Mesh -> None,
ColorFunction -> "Rainbow", AxesLabel -> Automatic]


I was expecting something more like the result from Alex's plot:

I can not see what I am doing wrong. Perhaps I need to scale the equations differently? Does anyone spot where my implementation goes south?

Complete code for my implementation:

ClearAll["Global*"]
(* SA eqn *)
saEQN = D[nut[t, x, y], t] + u[t, x, y]*D[nut[t, x, y], x] ==
cb1*(1 - ft2)*STilda*
nut[t, x,
y] - (cw1 fw - (cb1/\[Kappa]^2)* ft2)*(nut[t, x, y]/d)^2 + (1/
sigma (D[(nu + nut[t, x, y]) D[nut[t, x, y], x], x] +
D[(nu + nut[t, x, y]) D[nut[t, x, y], y], y]) +
cb2/sigma*(D[nut[t, x, y], x]^2 + D[nut[t, x, y], y]^2));
(* fluid eqns *)
mut = rho*nut[t, x, y]*fv1;
fluidEQN =
D[u[t, x, y], t] + u[t, x, y]*D[u[t, x, y], x] + px ==
D[(mu + mut) D[u[t, x, y], y], y] + D[u[t, x, y], x, x];

(* model parameters *)
sigma = 2/3; kappa = .41; cb1 = 0.1355; cb2 = 0.622; eps = 10^-6; d =
Sqrt[y^2 + eps^2]; cw1 = cb1/kappa^2 + (1 + cb2)/sigma;
cw2 = 0.3; cw3 = 2; cv1 = 7.1; ct1 = 1; ct2 = 2; ct3 = 1.2; ct4 = \
0.5;
mu = 1.711*10^-5;
rho = 1;
nu = mu/rho;
omega = 1/2 Sqrt[(D[u[t, x, y], y] - D[u[t, x, y], x])^2];
chi = nut[t, x, y]/nu;
fv1 = chi^3/(chi^3 + cv1^3);
fv2 = 1 - chi/(1 + chi*fv1);
Omega = Sqrt[2*omega*omega];
STilda = Omega + nut[t, x, y]/(kappa^2*d^2)*fv2;
g = r + cw2 (r^6 - r);
r = Min[nut[t, x, y]/(STilda*kappa^2*d^2), 10];
fw = g*((1 + cw3^6)/(g^6 + cw3^6))^(1/6);
ft2 = ct3 Exp[-ct4*chi^2];

(* eqns, bcs, ics *)
eq = {saEQN, fluidEQN};
bcUpdate = {
nut[t, x, 0] == 0, u[t, x, 0] == 0,
nut[t, x, L] == 0.1, u[t, x, L] == 1,
nut[t, 0, y] == 0.1 y/L, u[t, 0, y] == y/L,
Derivative[0, 1, 0][nut][t, L, y] == 0,
Derivative[0, 1, 0][u][t, L, y] == 0};
ic = {nut[0, x, y] == 0.1 y/L, u[0, x, y] == y/L};
(* geometry, time, pressure *)
L = 10^4;
t0 = 15;
px = 0;

(* solving *)
Monitor[{nutTGP, uTGP} =
NDSolveValue[{eq, ic, bcUpdate}, {nut, u}, {t, 0, t0}, {x, 0,
L}, {y, 0, L},
EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])
];, monitor]

(* plot *)
Plot3D[nutTGP[t0, x, y], {x, 0, L}, {y, 0, L}, Mesh -> None,
ColorFunction -> "Rainbow", AxesLabel -> Automatic]


Complete code for my Alex' implementation:

ClearAll["Global*"]
L = 10^4;
t0 = 15;
px = 0;
sigma = 2/3; kappa = .41; cb1 = .1355; cb2 = .622; eps = 10^-6; d =
Sqrt[y^2 + eps^2]; cw1 =
cb1/kappa^2 + (1 + cb2)/
sigma; cw2 = .3; cw3 = 2; cv1 = 7.1; ct1 = 1; ct2 = 2; ct3 = 1.2; \
ct4 = .5;
mu = 1.711 10^-5;
omega = Sqrt[(D[u[t, x, y], y] - D[u[t, x, y], x])^2];
chi = nu[t, x, y]/mu;
S = omega + (1 - chi/(1 + chi^4/(cv1^3 + chi^3))) nu[t, x,
y]/(kappa d)^2;
r = nu[t, x, y]/S/(kappa d)^2;
fw = (r +
cw2 (r^6 - r)) ((1 + cw3^6/(cw3^6 + (r + cw2 (r^6 - r))^6)))^(1/
6);
ft2 = ct3 Exp[-ct4 chi^2];
nut = nu[t, x, y] chi^3/(cv1^3 + chi^3);
eq = {mu*(D[nu[t, x, y], t] +
u[t, x, y] D[nu[t, x, y], x]) == (cb1 (1 - ft2) S nu[t, x,
y] - (cw1 fw -
cb1/kappa^2 ft2) (nu[t, x, y]/d)^2 + (1/
sigma (D[(nut + mu) D[nu[t, x, y], x], x] +
D[(nut + mu) D[nu[t, x, y], y], y]) +
cb2/sigma  (D[nu[t, x, y], x]^2 + D[nu[t, x, y], y]^2))),
D[u[t, x, y], t] + u[t, x, y] D[u[t, x, y], x] + px ==
D[u[t, x, y], x, x] + D[(nut/mu + 1) D[u[t, x, y], y], y]};
bcUpdate = {nu[t, x, 0] == 0, u[t, x, 0] == 0, u[t, x, L] == 1,
nu[t, x, L] == 0.1, nu[t, 0, y] == 0.1 y/L, u[t, 0, y] == y/L,
Derivative[0, 1, 0][nu][t, L, y] == 0,
Derivative[0, 1, 0][u][t, L, y] == 0};
ic = {nu[0, x, y] == 0.1 y/L, u[0, x, y] == y/L};
Monitor[{nuTPG1, UTPG1} =
NDSolveValue[{eq, ic, bcUpdate}, {nu, u}, {t, 0, t0}, {x, 0,
L}, {y, 0, L}
, EvaluationMonitor :> (monitor =
Row[{"t = ", CForm[t]}])];, monitor]
Plot3D[nuTPG1[t0, x, y], {x, 0, L}, {y, 0, L}, Mesh -> None,
ColorFunction -> "Rainbow", AxesLabel -> Automatic]

• It is a good approach (+1), but my post was about very special case of turbulent boundary layer. Commented Feb 6, 2023 at 12:08

There are several remarks concerning Spalart-Allmaras turbulence model (SA) for aerodynamic applications, published in AIAA Paper 92 - 0439, 1992. First, to compare SA with our turbulence model we need to rescale SA so that the space scale $$l$$ and velocity scale $$u_0$$ connect with dynamic viscosity $$\mu$$ as follows $$l u_0=\mu /\rho$$. Due to this rescaling we have equations for the flow velocity and turbulent viscosity in the form

saEQN = mu (D[nut[t, x, y], t] + u[t, x, y]*D[nut[t, x, y], x]) ==
cb1*(1 - ft2)*STilda*
nut[t, x,
y] - (cw1 fw - (cb1/\[Kappa]^2)*ft2)*(nut[t, x, y]/d)^2 + (1/
sigma (D[(mu + mut) D[nut[t, x, y], x], x] +
D[(mu + mut) D[nut[t, x, y], y], y]) +
cb2/sigma*(D[nut[t, x, y], x]^2 + D[nut[t, x, y], y]^2));

fluidEQN =
D[u[t, x, y], t] + u[t, x, y]*D[u[t, x, y], x] + px ==
D[(1 + mut/mu) D[u[t, x, y], y], y] + D[u[t, x, y], x, x];



where mut = rho*nut[t, x, y]*fv1;. Second, in the wall region we use omega = Sqrt[(D[u[t, x, y], y] - D[u[t, x, y], x])^2]; and STilda = omega + nut[t, x, y]/(kappa^2*d^2)*fv2;. Third, we use the homogenous Neuman bc at x==L to avoid message from NDSolve. Finally, the code is given by

ClearAll["Global*"]

sigma = 2/3; kappa = .41; cb1 = 0.1355; cb2 = 0.622; eps = 10^-6; d =
Sqrt[y^2 + eps^2]; cw1 = cb1/kappa^2 + (1 + cb2)/sigma;
cw2 = 0.3; cw3 = 2; cv1 = 7.1; ct1 = 1; ct2 = 2; ct3 = 1.2; ct4 = \
0.5;
mu = 1.711*10^-5;
rho = 1;
nu = mu/rho;

omega = Sqrt[(D[u[t, x, y], y] - D[u[t, x, y], x])^2];
chi = nut[t, x, y]/nu;
fv1 = chi^3/(chi^3 + cv1^3);
fv2 = 1 - chi/(1 + chi*fv1);
STilda = omega + nut[t, x, y]/(kappa^2*d^2)*fv2;
g = r + cw2 (r^6 - r);
r = nut[t, x, y]/(STilda*kappa^2*d^2);
fw = g*((1 + cw3^6)/(g^6 + cw3^6))^(1/6);
ft2 = ct3 Exp[-ct4*chi^2];

saEQN = mu (D[nut[t, x, y], t] + u[t, x, y]*D[nut[t, x, y], x]) ==
cb1*(1 - ft2)*STilda*
nut[t, x,
y] - (cw1 fw - (cb1/\[Kappa]^2)*ft2)*(nut[t, x, y]/d)^2 + (1/
sigma (D[(mu + mut) D[nut[t, x, y], x], x] +
D[(mu + mut) D[nut[t, x, y], y], y]) +
cb2/sigma*(D[nut[t, x, y], x]^2 + D[nut[t, x, y], y]^2));
mut = rho*nut[t, x, y]*fv1;
fluidEQN =
D[u[t, x, y], t] + u[t, x, y]*D[u[t, x, y], x] + px ==
D[(1 + mut/mu) D[u[t, x, y], y], y] + D[u[t, x, y], x, x];

eq = {saEQN, fluidEQN}; bcUpdate = {nut[t, x, 0] == 0,
u[t, x, 0] == 0, nut[t, x, L] == 0.1, u[t, x, L] == 1,
nut[t, 0, y] == 0.1 y/L, u[t, 0, y] == y/L,
Derivative[0, 1, 0][nut][t, L, y] == 0,
Derivative[0, 1, 0][u][t, L, y] == 0}; ic = {nut[0, x, y] ==
0.1 y/L, u[0, x, y] == y/L}; L = 10^4;
t0 = 15;
px = 0;
Monitor[{nutTGP, uTGP} =
NDSolveValue[{eq, ic, bcUpdate}, {nut, u}, {t, 0, t0}, {x, 0,
L}, {y, 0, L},
EvaluationMonitor :> (monitor =
Row[{"t = ", CForm[t]}])];, monitor]


Visualization

{Plot3D[nutTGP[t0, x, y], {x, 0, L}, {y, 0, L}, Mesh -> None,
ColorFunction -> "Rainbow", AxesLabel -> Automatic],
Plot3D[uTGP[t0, x, y], {x, 0, L}, {y, 0, L}, Mesh -> None,
ColorFunction -> "Rainbow", AxesLabel -> Automatic]}
`

• (+1) Thanks for you time Alex. Please allow for a few days to study this. Commented Feb 6, 2023 at 8:47