5
$\begingroup$

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]

enter image description here

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

enter image description here

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]
Improve this question
$\endgroup$
2
  • $\begingroup$ It is a good approach (+1), but my post was about very special case of turbulent boundary layer. $\endgroup$
    – Alex Trounev
    Commented Feb 6, 2023 at 12:08
  • $\begingroup$ @Alex, yes I start to understand that now. I am thinking about a more general approach for turbulence modeling. I'll think more about this but your post has already been helpful. $\endgroup$
    – user21
    Commented Feb 6, 2023 at 12:13

1 Answer 1

Reset to default
2
$\begingroup$

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

Figure 1

Improve this answer
$\endgroup$
1
  • $\begingroup$ (+1) Thanks for you time Alex. Please allow for a few days to study this. $\endgroup$
    – user21
    Commented Feb 6, 2023 at 8:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.