Some conceptual doubts about steady and transient solvers of the NS + energy equations

Question

Recently, I have been solving some transient and steady, flow and heat transfer problems in Mathematica. The transient problem is essentially the reciprocating (i.e., fully reversing) flow of a fluid over a thick heated block until the system reaches a cyclic steady-state (i.e., the system temperature oscillates around a mean when this point is reached). I received help regarding this in one of my earlier question on this site. Following is the code:

Needs["NDSolve`FEM`"]
{f = 0.5; L = 0.040, d = 0.003, e = 0.005, kf = 0.026499, ks = 16, 
 rho = 1.1492, rhos = 7860, mu = 18.92310^-6, cp = 1.006910^3, 
 cps = 502.4}; u0 = 3; nu = mu/rho; om = 2 Pi f;
tflow = 100;
t0 = .3;
NV = 2 f tflow;
nn = Round[NV [Pi]/(om t0)]
Ti = 307; q = 1000/Ti;
reg1 = ImplicitRegion[0 <= x <= L && 0 <= y <= d, {x, y}]; reg2 = 
 ImplicitRegion[0 <= x <= L && -e <= y <= d, {x, y}];
UX[0][x_, y_] := 0;
VY[0][x_, y_] := 0;
P[0][x_, y_] := 0;
Tfs[0][x_, y_] := 307/Ti; appro = 
 With[{k = 2. 10^6}, ArcTan[k #]/Pi + 1/2 &];
ade[y_] := (ks + (kf - ks) appro[y])
rde[y_] := (cps rhos + (cp rho - cps rhos) appro[y]);
Monitor[Do[{UX[i], VY[i], P[i]} = 
    NDSolveValue[{{Inactive[
            Div][({{-[Mu], 0}, {0, -[Mu]}}.Inactive[Grad][
              u[x, y], {x, y}]), {x, y}] + D[p[x, y], x] + 
          UX[i - 1][x, y]D[u[x, y], x] + 
          VY[i - 1][x, y]D[u[x, y], y] + (u[x, y] - UX[i - 1][x, y])/
           t0, Inactive[
        Div][({{-\[Mu], 0}, {0, -\[Mu]}}.Inactive[Grad][
          v[x, y], {x, y}]), {x, y}] + D[p[x, y], y] + 
      UX[i - 1][x, y]*D[v[x, y], x] + 
      VY[i - 1][x, y]*D[v[x, y], y] + (v[x, y] - VY[i - 1][x, y])/
       t0, D[u[x, y], x] + D[v[x, y], y]} == {0, 0, 0} /. \[Mu] -&gt;
     nu, {DirichletCondition[{u[x, y] == u0*Sin[om*i*t0], 
     v[x, y] == 0}, 
    x == L (1 - Sign[Sin[om*i*t0]])/2 &amp;&amp; 0 &lt; y &lt; d], 
   DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, 
    y == 0 || y == d]}, 
  DirichletCondition[p[x, y] == 0, 
   x == L (1 + Sign[Sin[om*i*t0]])/2 &amp;&amp; 0 &lt; y &lt; d]}, {u, v, 
  p}, {x, y} \[Element] reg1, 
 Method -&gt; {&quot;FiniteElement&quot;, 
   &quot;InterpolationOrder&quot; -&gt; {u -&gt; 2, v -&gt; 2, p -&gt; 1}, 
   &quot;MeshOptions&quot; -&gt; {&quot;MaxCellMeasure&quot; -&gt; 0.0000005}}];

ux = If[y <= 0, 0, UX[i][x, y]]; vy = If[y <= 0, 0, VY[i][x, y]];
   Tfs[i] = 
    NDSolveValue[{rde[
           y] ((uxD[T[x, y], x] + 
              vyD[T[x, y], y]) + (T[x, y] - Tfs[i - 1][x, y])/t0) - 
         Inactive[Div][
          ade[y]Inactive[Grad][T[x, y], {x, y}], {x, y}] == 
        NeumannValue[q, y == -e], 
       DirichletCondition[{T[x, y] == 1}, 
        x == L (1 - Sign[Sin[omi*t0]])/2 && 0 <= y <= d]}, 
      T, {x, y} [Element] reg2, 
      Method -> {"FiniteElement", "InterpolationOrder" -> {T -> 2}, 
        "MeshOptions" -> {"MaxCellMeasure" -> 0.0000001}}] // 
     Quiet;, {i, 1, nn}],ProgressIndicator[i,{1,nn}]]; // AbsoluteTiming
ListLinePlot[
 Table[{i t0, Tfs[i][0.5*L, -e/2]*Ti - 273.16}, {i, 0, nn}], 
 AxesLabel -> {"t(s)", "T(K)"}, PlotRange -> Full]

The above code runs the simulation for 100s flow-time for a flow oscillating with a frequency of 0.5Hz using the velocity profile $u = 3\sin(2*\pi*0.5*t)$. As can be seen in the above code the Do loop marches forward in time. Inside the loop the flow and temperature field is being solved at each iteration nn.

I modified the above code for a uni-directional flow, i.e., find a steady-state solution for a flow from left to right over the heated block. I removed the time-dependent terms inside the loop and the nn, instead of a time-stepper becomes an iterator to reach a converged solution. The code is:

Needs["NDSolve`FEM`"]
{f = 0.5; L = 0.040, d = 0.003, e = 0.005, kf = 0.026499, ks = 16, 
 rho = 1.1492, rhos = 7860, mu = 18.92310^-6, cp = 1.006910^3, 
 cps = 502.4}; u0 = 3; nu = mu/rho; om = 2 Pi f;
Ti = 307; q = 1000/Ti;
uavg = (1/(Pi/om)) Integrate[u0Sin[omt], {t, 0, Pi/om}] // N
nn = 10;
reg1 = ImplicitRegion[0 <= x <= L && 0 <= y <= d, {x, y}]; reg2 = 
 ImplicitRegion[0 <= x <= L && -e <= y <= d, {x, y}];
UX[0][x_, y_] := 0;
VY[0][x_, y_] := 0;
P[0][x_, y_] := 0;
Tfs[0][x_, y_] := 307/Ti; appro = 
 With[{k = 2. 10^6}, ArcTan[k #]/Pi + 1/2 &];
ade[y_] := (ks + (kf - ks) appro[y])
rde[y_] := (cps rhos + (cp rho - cps rhos) appro[y]);
Do[{UX[i], VY[i], P[i]} = 
    NDSolveValue[{{Inactive[
            Div][({{-[Mu], 0}, {0, -[Mu]}}.Inactive[Grad][
              u[x, y], {x, y}]), {x, y}] + D[p[x, y], x] + 
          UX[i - 1][x, y]D[u[x, y], x] + 
          VY[i - 1][x, y]D[u[x, y], y], 
         Inactive[
            Div][({{-[Mu], 0}, {0, -[Mu]}}.Inactive[Grad][
              v[x, y], {x, y}]), {x, y}] + D[p[x, y], y] + 
          UX[i - 1][x, y]D[v[x, y], x] + 
          VY[i - 1][x, y]D[v[x, y], y], 
         D[u[x, y], x] + D[v[x, y], y]} == {0, 0, 0} /. [Mu] -> 
        nu, {DirichletCondition[{u[x, y] == uavg, v[x, y] == 0}, 
        x == 0 && 0 < y < d], 
       DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, 
        y == 0 || y == d]}, 
      DirichletCondition[p[x, y] == 0, x == L && 0 < y < d]}, {u, v, 
      p}, {x, y} [Element] reg1, 
     Method -> {"FiniteElement", 
       "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}, 
       "MeshOptions" -> {"MaxCellMeasure" -> 0.0000005}}];
   ux = If[y <= 0, 0, UX[i][x, y]]; vy = If[y <= 0, 0, VY[i][x, y]];
   Tfs[i] = 
    NDSolveValue[{rde[y] ((uxD[T[x, y], x] + vyD[T[x, y], y])) - 
         Inactive[Div][
          ade[y]*Inactive[Grad][T[x, y], {x, y}], {x, y}] == 
        NeumannValue[q, y == -e], 
       DirichletCondition[{T[x, y] == 1}, x == 0 && 0 <= y <= d]}, 
      T, {x, y} [Element] reg2, 
      Method -> {"FiniteElement", "InterpolationOrder" -> {T -> 2}, 
        "MeshOptions" -> {"MaxCellMeasure" -> 0.0000001}}] // 
     Quiet;, {i, 1, nn}]; // AbsoluteTiming

My question is:

Why do we need an iterator to reach a steady-state solution?

In other words, the time-dependent code solved for the flow and temperature field at each time-step (but within each time-step there seemed to be no need for further iterations to reach a converged solution). Going by that logic, NDSolve should have been able to directly solve for the steady-state solution (i.e, without the time-dependent terms), but it cannot as I tried below:

{U, V, P} = 
 NDSolveValue[{{Inactive[
         Div][({{-\[Mu], 0}, {0, -\[Mu]}}.Inactive[Grad][
           u[x, y], {x, y}]), {x, y}] + D[p[x, y], x] + 
       u[x, y]*D[u[x, y], x] + v[x, y]*D[u[x, y], y], 
      Inactive[
         Div][({{-\[Mu], 0}, {0, -\[Mu]}}.Inactive[Grad][
           v[x, y], {x, y}]), {x, y}] + D[p[x, y], y] + 
       u[x, y]*D[v[x, y], x] + v[x, y]*D[v[x, y], y], 
      D[u[x, y], x] + D[v[x, y], y]} == {0, 0, 0} /. \[Mu] -> 
     nu, {DirichletCondition[{u[x, y] == uavg, v[x, y] == 0}, 
     x == 0 && 0 < y < d], 
    DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, 
     y == 0 || y == d]}, 
   DirichletCondition[p[x, y] == 0, x == L && 0 < y < d]}, {u, v, 
   p}, {x, y} \[Element] reg1, 
  Method -> {"FiniteElement", 
    "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}, 
    "MeshOptions" -> {"MaxCellMeasure" -> 0.0000005}

I would like to have a descriptive explanation of this difference between modelling a transient system and a steady-state solution? In essence, why within each time-step in the transient solver there is no need of iteration ? However, a steady-state solution requires one to iterate. I hope I could clearly explain my query.

Please, pay attention, that in my code with UX[i], VY[i], P[i] we solve on every step the linear system of equations, while in your code with {U, V, P} you try to solve nonlinear system. This is a big difference for NDSolve, since we use FEM. — Alex Trounev, Dec 10 '22 at 11:52
Thankyou for your comment Alex. I do understand that in my last code block, I try to directly solve a highly non-linear system with significant convective terms. However, when I modified your original code to solve for the steady-problem (by removing the time terms, i.e., the second code block), it worked fine. Hence, I wanted to understand this basic difference between modelling a steady-state and a time-dependent system. It seems my understanding is flawed. — Avrana, Dec 10 '22 at 12:55
In my code the false transient method is implemented, it is why numerical solution converges to steady state in a case of the time term elimination - see explanation on https://community.wolfram.com/groups/-/m/t/1433064?p_p_auth=4O8xgpur — Alex Trounev, Dec 10 '22 at 13:59
@AlexTrounev Thankyou for the answer. I went to the link and found some helpful notes. One of your codes there on square driven lid cavity used the method of false transients, which is helpful. I have one other request. Can you run the first code block (i.e., reciprocating flow) for these parameter values f=1, t0=0.15, tflow=700 (Other parameters remain same). I have seen that whenever I am giving a large tflow, the calculation stops midway with a kernel error, although I am running the calculation on a system with 16GB of memory. — Avrana, Dec 11 '22 at 11:29
Yes, you are right, there are some unstably solutions. In your case just put t0=0.3 to compute solution. Regardless your method at $t0\rightarrow \infty$, this is also nice method for some cases. — Alex Trounev, Dec 12 '22 at 04:54
@AlexTrounev Thankyou for verifying that the problem actually exists. Now, let me try to explain why I cannot keep t0=0.3. If my flow is oscillating/reciprocating at $1Hz$, that means the time period is $1s$, which implies that the flow changes direction every $0.5s$. Hence, to resolve this change in flow direction, I need to keep my step-size $<0.5s$. Keeping it t0=0.3, obviously satisfies this criterion but it means that I am only calculating at one time-step in each half-cycle, which is not a very good resolution. I will prefer at least $2$ time-steps, which t0=0.15 achieves. — Avrana, Dec 12 '22 at 05:55
In any case, if my flow oscillates at $f=2Hz$, I will need a t0=0.1 to properly resolve the flow reversal. The problem is, I need to run these calculations till a large tflow, until the system reaches a cyclic steady state, i.e., the temperature of the system (for ex. any point in the solid domain) oscillates around a mean value. And this tflow can get quite large if my heat input, i.e., q increases. — Avrana, Dec 12 '22 at 05:59
My modifications for $t0\rightarrow \infty$, was for a situation, where there is no flow oscillation at all. Instead, the flow just occurs from $x=0$ to $x=L$. And I was aiming for a steady-state temperature distribution. Do you think, there can be a way by which one can directly calculate the cyclic-steady state, when the flow velocity is oscillating, instead of marching forward through the entire transient period ? Should these be asked as a separate question ? — Avrana, Dec 12 '22 at 06:03
Ok, your code working with uavg=1.91. What about your question? — Alex Trounev, Dec 12 '22 at 09:27
I have posted a separate question here. Have a look. And thankyou for your help here. — Avrana, Dec 12 '22 at 11:03

Some conceptual doubts about steady and transient solvers of the NS + energy equations

0 Answers0