Potential Singularity of the 3D Euler Equations in the Interior Domain

Abstract

Whether the 3D incompressible Euler equations can develop a finite time singularity from smooth initial data is one of the most challenging problems in nonlinear PDEs. In this paper, we present some new numerical evidence that the 3D axisymmetric incompressible Euler equations with smooth initial data of finite energy develop a potential finite time singularity at the origin. This potential singularity is different from the blow-up scenario revealed by Luo and Hou (111:12968–12973, 2014) and (12:1722–1776, 2014), which occurs on the boundary. Our initial condition has a simple form and shares several attractive features of a more sophisticated initial condition constructed by Hou and Huang in (arXiv:2102.06663, 2021) and (435:133257, 2022). One important difference between these two blow-up scenarios is that the solution for our initial data has a one-scale structure instead of a two-scale structure reported in Hou and Huang (arXiv:2102.06663, 2021) and (435:133257, 2022). More importantly, the solution seems to develop nearly self-similar scaling properties that are compatible with those of the 3D Navier–Stokes equations. We will present numerical evidence that the 3D Euler equations seem to develop a potential finite time singularity. Moreover, the nearly self-similar profile seems to be very stable to the small perturbation of the initial data.

Appendix A: Construction of the Adaptive Mesh

In this appendix, we describe our adaptive mesh strategy to study the singularity formation near the origin $(r,z) = (0,0)$. We will use the method described in Appendix B of [20] to construct our adaptive mesh maps $r=r(\rho )$ and $z=z(\eta )$. We will discretize the equations in the transformed variables $(\rho ,\eta )$ with $n_1$ grid points along the z direction and $n_2$ grid points along the r-direction.

The adaptive mesh strategy described in [20] was inspired by the adaptive mesh strategy introduced in [32]. On the other hand, the adaptive mesh strategy developed in [32] is simpler due to the fact that the singularity is located at a fixed stagnation point on the boundary $(r,z)=(1,0)$, $\omega _1 >0$ for $z>0$ and has a bell-shaped structure near the singularity. In our case, we have a traveling wave solution that approaches the origin with $\omega _1$ changing sign in the most singular region. The singular region has a more complicated shape since the potentially singular solution produces a strong shearing flow traveling downstream. Thus, we need to introduce a moving frame and design our adaptive mesh map $r(\rho )$ and $z(\eta )$ to resolve the solution in different regions. Our adaptive mesh strategy does not require that the solution has a bell-shaped structure in the most singular region. To construct $r(\rho )$, we use the distance dr between the location at which $u_1$ achieves its maximum and the location at which $u_{1r}$ achieves its maximum to define the boundary $r_i$ for different singular regions (phases). Similarly, we construct $z(\eta )$ by using the distance dz between the location at which $\omega _1$ achieves its maximum and the location at which $\omega _{1z}$ achieves its maximum to define the boundary $z_j$ for different phases.

1.1 A.1. The Adaptive (Moving) Mesh Algorithm

To effectively and accurately compute the potential blow-up, we have designed a special meshing strategy that is dynamically adaptive to the singular structure of the solution. The adaptive mesh covering the half-period computational domain $\mathcal {D}_1 = \{(r,z):0\le r\le 1,0\le z\le 1/2\}$ is characterized by a pair of analytic mesh mapping functions

$$\begin{aligned}r = r(\rho ),\quad \rho \in [0,1];\quad z = z(\eta ),\quad \eta \in [0,1].\end{aligned}$$

These mesh mapping functions are both monotonically increasing and infinitely differentiable on [0, 1], and satisfy $r(0) = 0,\; r(1) = 1,\; z(0) = 0,\; z(1) = 1/2.$ In particular, we construct these mapping functions by carefully designing their Jacobians/densities $r_\rho = r'(\rho ),\quad z_\eta = z'(\eta )$ using analytic functions that are even functions at 0. The even symmetries ensure that the resulting mesh can be smoothly extended to the full-period cylinder $\{(r,z):0\le r\le 1,-1/2\le z\le 1/2\}$. The density functions contain a small number of parameters, which are dynamically adjusted to the solution. Once the mesh mapping functions are constructed, the computational domain is covered with a tensor-product mesh:

$$\begin{aligned} \mathcal {G} = \{(r_i,z_j): 0\le i\le n_2,\ 0\le j\le n_1\}, \end{aligned}$$

(A.1)

where $r_i^h = r(ih_\rho ),\; h_\rho = 1/n_2;\; z_j^h = z(jh_\eta ),\; h_\eta = 1/n_1.$ The precise definition and construction of the mesh mapping functions are described in Appendix B of [20].

Figure 30 gives an example of the densities $r_\rho ,z_\eta$ (in log scale) we use in the computation. We design the densities $r_\rho ,z_\eta$ to have three phases:

• Phase 1 covers the inner profile of the smaller scale near the sharp front;

• Phase 2 covers the outer profile of the larger scale of the solution;

• Phase 3 covers the (far-field) solution away from the symmetry axis $r=0$.

We add a phase 0 in the density $r_\rho$ to cover the region near $r=0$ and also add a phase 0 in the density $z_\eta$ to cover the region near $z=0$ in the late stage. In our computation, the number (percentage) of mesh points in each phase are fixed, but the physical location of each phase will change in time, dynamically adaptive to the structure of the solution. Between every two neighboring phases, there is also a smooth transition region that occupies a fixed percentage of mesh points.

Fig. 30

The mapping densities $r_\rho$ (left) and $z_\eta$ (right) with phase numbers labeled. This figure is for illustration only. The parameters do not reflect the adaptive mesh used in our computation

1.2 A.2. Adaptive Mesh for the 3D Euler Equations

We use three different adaptive mesh strategies for three different time periods. The first time period corresponds to the time interval between $t=0$ and $t_1=0.002231338$ with $\Vert \omega (t_1)\Vert _{L^\infty }/\Vert \omega (0)\Vert _{L^\infty } \approx 46.54325$ for the $1536\times 1536$ grid and the number of time steps equal to 45000. The second time period corresponds to the time interval between $t_1=0.002231338$ and $t_2 = 0.002264353$ with $\Vert \omega (t_2)\Vert _{L^\infty }/\Vert \omega (0)\Vert _{L^\infty } \approx 295.39986$ for the $1536\times 1536$ grid and the number of time steps equal to 60000. The third time period is for $t \ge t_2$.

In the first time period, since we use a very smooth initial condition whose support covers the whole domain, we use the following parameters $r_1=0.001,\; r_2=0.05, \; r_3=0.2$ and $s_{\rho _1}=0.001$, $s_{\rho _2}=0.5$, $s_{\rho _3}=0.85$ to construct the mapping $r=r(\rho )$ using a four-phase map. Similarly, we use the following parameters $z_1=0.1,\; z_2=0.25$ and $s_{\eta _1}=0.5$, $s_{\eta _2}=0.85$ to construct the mapping $z=z(\eta )$ using a three-phase map. We then update the mesh $z=z(\eta )$ dynamically using $z_1=2 z(I)$ and $z_2=10 z(I)$ with $s_{\eta _1}=0.6$, $s_{\eta _2}=0.9$ when $I < 0.25n_1$, but keep $r=r(\rho )$ unchanged during this early stage. Here I is the grid point index along the z-direction at which $\omega _1$ achieves its maximum.

In the second time period, we use the following parameters $s_{\rho _1}=0.05$, $s_{\rho _2}=0.6$, $s_{\rho _3}=0.9$, $r_2 = r(J) + 2dr$, $r_1=\max ((s_{\rho _1}/s_{\rho _2})r_2,r(J_r) - 5dr)$, and $r_3 = \max (3r(J),(r_2 - r_1)(s_{\rho _3}-s_{\rho _2})/(s_{\rho _2}-s_{\rho _1}) + r_2)$, where J is the grid index at which $u_1$ achieves its maximum along the r-direction, $J_r$ is the grid index at which $u_{1,r}$ achieves its maximum along the r-direction, and $dr = r(J) - r(J_r)$. We update the mapping $r(\rho )$ dynamically when $J_r < 0.2n_2$. The adaptive mesh map for $z(\eta )$ in the second time period remains the same as in the first time period.

In the third time period, we need to allocate more grid points to resolve the sharp front. We use the following parameters $s_{\rho _1}=0.05$, $s_{\rho _2}=0.65$, $s_{\rho _3}=0.9$, $r_2 = r(J) + 10dr$, $r_1=\max ((s_{\rho _1}/s_{\rho _2})r_2,r(J_r) - 3dr)$, and $r_3 = \max (2.3r(J),(r_2 - r_1)(s_{\rho _3}-s_{\rho _2})/(s_{\rho _2}-s_{\rho _1}) + r_2)$. To construct the mesh map $z(\eta$, we use the following parameters $s_{\eta _1}=0.05$, $s_{\eta _2}=0.65$, $s_{\eta _3}=0.9$, $z_2 = z(I_w) + 2dz$, $z_1=\max ((s_{\eta _1}/s_{\eta _2})z_2,z(I_{wz}) - 16dz)$, and $z_3 = \max (2.3z(I_w),(z_2 - z_1)(s_{\eta _3}-s_{\eta _2})/(s_{\eta _2}-s_{\eta _1}) + z_2)$, where $I_w$ is the grid index at which $\omega _1$ achieves its maximum along the z-direction, $I_{wz}$ is the grid index at which $\omega _{1,z}$ achieves its maximum along the z-direction, and $dz = z(I_w) - r(I_{wz})$. We will update $r(\rho )$ dynamically when $J_r < 0.2 n_2$ and update $z(\eta )$ when $I_z < 0.23 n_1$.

The three time periods for different grids will be defined similarly through the number of time steps by using a linear scaling relationship. For example, for the $1024\times 1024$ grid, the first time period will be between $t=0$ and the time step $30000=45000 \cdot (1024/1536)$. The second time period will be between time steps 30000 and $40000 = 60000\cdot (1024/1536)$. The third time period will be beyond time step 40000. The three time periods for other grids are defined similarly.

As we mentioned in the Introduction, we use a second-order Runge-Kutta method to discretize the Euler and Navier–Stokes equations in time with an adaptive time stepping strategy. We discretize the Euler and Navier–Stokes equations in the transformed domain $(\eta ,\rho )$ using a uniform mesh with $h_1=1/n_1$ and $h_2=1/n_2$. We choose our adaptive time step as follows:

$$\begin{aligned} k= & {} \min (k_1,k_2), \;\;k_1 = \min (\;0.2\min (h1,h_2)/\text{umax },\;10^{-3}\Vert u_1\Vert _{L^\infty }^{-1},\;2.5\cdot 10^{-7}),\\ k_2= & {} 0.1\min (\min (h_1z_\eta )^2/\nu ,\min (h_2r_\rho )^2/\nu ), \end{aligned}$$

where $\text{ umax } = \max (\Vert u^r/r_\rho \Vert _{L^\infty },\Vert u^z/z_\eta \Vert _{L^\infty })$ is the maximum velocity in the transformed domain.