Infinite expresssion 1/0 encountered when use NDSolve for 3D axisymmetric Navier-Stokes (Euler) equations

Question

The PDEs we are interested in solving using NDSolve is the vorticity-stream formulation of the 3D axisymmetric Navier-Stokes (Euler) equations (Ref.1 :T. Y. Hou, Potential singularity of the 3D Euler equations in the interior domain, arXiv:2107.05870): \begin{equation}\label{NSvorticityStreamD2} u_{t}+(-xf_{z}) u_{x}+(2f+xf_{x}) u_{z} =\nu \left(u_{xx} +u_{zz}+\tfrac{3}{x}u_{x}\right) +2uf_{z}\\ % w_{t}+(-xf_{z}) w_{x}+(2f+xf_{x}) w_{z} =\nu \left(w_{xx}+\omega_{zz}+\tfrac{3}{x}w_{x}\right)+2u u_{z},\\ % -w=\,\,\,\,\left(f_{xx} +f_{zz}+\tfrac{3}{x}f_{x} \tag{1}\right). \end{equation}

The initial conditions are \begin{equation}\label{InitialConditions} u(0,x,z)=\frac{12000(1-x^2)^{18}\sin(2\pi z)}{1+12.5\sin^2(\pi z)},\\ w(0,x,z)=f(0,x,z)=0 \tag{2}. \end{equation}

The boundary conditions are \begin{equation}\label{BoundaryConditionD2} u(t,\pm1,z)=f(t,\pm1,z)=w(t,\pm1,z)=0.\\ u(t,x,\pm \tfrac12)=f(t,x,\pm \tfrac12)=w(t,x,\pm \tfrac12)=0. \end{equation}

When I used the following NDSolve command

NDSolve[{PDE, IC, BC}, {u, w, f}, {t, 0, 1/10}, {x, -1, 1}, {z, -1/2, 1/2}]

I got the error message:

Infinite expression 1/0 encountered. 

I will add Mathematica code tonight. Here is the Mathematica code:

PDE = {-2*u[t, x, z]*Derivative[0, 0, 1][f][t, x, z] + 
    Derivative[0, 0, 1][u][t, x, z]*
      (2*f[t, x, z] + x*Derivative[0, 1, 0][f][t, x, 
             z]) - x*Derivative[0, 0, 1][f][t, x, z]*
      Derivative[0, 1, 0][u][t, x, z] + 
    Derivative[1, 0, 0][u][t, x, z] == 0, 
-2*u[t, x, z]*Derivative[0, 0, 1][u][t, x, z] + 
    Derivative[0, 0, 1][w][t, x, z]*
      (2*f[t, x, z] + x*Derivative[0, 1, 0][f][t, x, 
             z]) - x*Derivative[0, 0, 1][f][t, x, z]*
      Derivative[0, 1, 0][w][t, x, z] + 
    Derivative[1, 0, 0][w][t, x, z] == 0, 
w[t, x, z] + Derivative[0, 0, 2][f][t, x, z] + 
    (3*Derivative[0, 1, 0][f][t, x, z])/x + 
    Derivative[0, 2, 0][f][t, x, z] == 0}
IC = {w[0, x, z] == 0, f[0, x, z] == 0, 
 u[0, x, z] == 12000(1 - x^2)^18
     (Sin[2Piz]/(1 + (25/2)Sin[Piz]^2))}
BC = {w[t, x, 1/2] == 0, f[t, x, 1/2] == 0, 
 u[t, x, 1/2] == 0, w[t, x, -2^(-1)] == 0, 
 f[t, x, -2^(-1)] == 0, u[t, x, -2^(-1)] == 0, 
 w[t, 1, z] == 0, f[t, 1, z] == 0, u[t, 1, z] == 0, 
 w[t, -1, z] == 0, f[t, -1, z] == 0, 
 u[t, -1, z] == 0}
NDSolve[{PDE, IC, BC}, {u, w, f}, {t, 0, 1/10}, {x, -1, 1}, {z, -1/2, 1/2}]

From inspection, we know that if $u(0,x,z),w(0,x,z),f(0,x,z)$ are even in $x$, then the solution $u(t,x,z),w(t,x,z),f(t,x,z)$ will remain to be even in $x$. So the coordinate singular term $(1/x)u_x$ is nonsingular as $x\to 0$. NDSolve is not capable of figuring out this fact.

Any suggestions?

Best regards- Mike

Answer 1

We can solve this problem using implicit-explicit method and linear FEM as follows. First, we add artificial viscosity of about $10^{-6}$ so that $\nu\ne 0$. Note that artificial viscosity also used in the paper even Thomas Hou supposed $\nu=0$. Second, we use time step of about $dt=0.00005675$, while in the paper $dt\le 10^{-6}$. Third, our mesh is of 4050 quad elements only, while in the paper there is adaptive mesh by $1536\times 1536$ elements. Nevertheless we can reproduce singularity similar as in the paper.

Needs["NDSolve`FEM`"];
mesh = ToElementMesh[Rectangle[{-1, -1/2}, {1, 1/2}], 
  AccuracyGoal -> 5, PrecisionGoal -> 5, "MaxCellMeasure" -> 0.0005] ;
U[0][x_, z_] := 
 12000(1 - x^2)^18(Sin[2Piz]/(1 + (25/2)Sin[Piz]^2)); 
W[0][x_, z_] := 0; F[0][x_, z_] := 0;
nu = 10^-6; dt = 
 0.00227/40; T = {0.002271815, 0.002274596, 0.002276480};
Do[{U[i], W[i], F[i]} = 
   NDSolveValue[{(u[x, z] - U[i - 1][x, z])/dt - 
       xDerivative[0, 1][F[i - 1]][x, z]Derivative[1, 0][u][x, z] + 
       Derivative[0, 1][u][x, 
         z](2F[i - 1][x, z] + xDerivative[1, 0][F[i - 1]][x, z]) - 
       2u[x, z]Derivative[0, 1][F[i - 1]][x, z] - 
       nu (Derivative[0, 2][u][x, z] + (3Derivative[1, 0][u][x, z])/
           x + Derivative[2, 0][u][x, z]) == 
      0, -2U[i - 1][x, z]Derivative[0, 1][u][x, z] + 
       Derivative[0, 1][w][x, 
         z](2F[i - 1][x, z] + xDerivative[1, 0][F[i - 1]][x, z]) - 
       xDerivative[0, 1][F[i - 1]][x, z]*
        Derivative[1, 0][w][x, z] + (w[x, z] - W[i - 1][x, z])/dt - 
       nu (Derivative[0, 2][w][x, z] + (3Derivative[1, 0][w][x, z])/
           x + Derivative[2, 0][w][x, z]) == 0, 
     w[x, z] + (Derivative[0, 2][f][x, 
          z] + (3Derivative[1, 0][f][x, z])/x + 
         Derivative[2, 0][f][x, z]) == 0, 
     DirichletCondition[{w[x, z] == 0, f[x, z] == 0, u[x, z] == 0}, 
      True]}, {u, w, f}, Element[{x, z}, mesh], 
    Method -> {"FiniteElement", 
      "InterpolationOrder" -> {u -> 2, w -> 2, f -> 2}}];, {i, 1, 40}]

Visualization

Table[Plot3D[U[i][x, z], Element[{x, z}, mesh], PlotRange -> All, 
  Mesh -> None, ColorFunction -> "Rainbow", PlotLabel -> 1. dt i], {i,
   0, 40, 10}]
Table[Plot3D[W[i][x, z], Element[{x, z}, mesh], PlotRange -> All, 
  Mesh -> None, ColorFunction -> "Rainbow", PlotLabel -> 1. dt i], {i,
   10, 40, 10}]

The detailed picture of singularity

{Plot3D[U[40][x, z], {x, 0, .04}, {z, 0, .04}, PlotRange -> All, 
  Mesh -> None, ColorFunction -> "Rainbow"], 
 Plot3D[W[40][x, z], {x, 0, .04}, {z, 0, .04}, PlotRange -> All, 
  Mesh -> None, ColorFunction -> "Rainbow"]}