1. Introduction
Mesoscopic methods provide a powerful framework that bridges the gap between microscopic molecular dynamics and macroscopic continuum mechanics [1–3]. While microscopic models resolve the dynamics of particles and thus incur prohibitively high computational costs for large systems, macroscopic models describe the evolution of averaged quantities such as density, velocity, and pressure with limited ability to capture nonequilibrium effects. In contrast, mesoscopic kinetic models, typically based on the Boltzmann or Enskog equations, evolve distribution functions that retain essential kinetic information while maintaining computational efficiency. These models therefore provide a physically rigorous yet numerically tractable bridge between the microscopic and macroscopic descriptions of fluid behavior. Among various mesoscopic approaches, the lattice Boltzmann method (LBM) [2, 4] and the discrete Boltzmann method (DBM) [3, 5] have emerged as two representative schemes capable of efficiently simulating complex nonequilibrium and multiscale flow phenomena.
Since the advent of lattice-gas or cellular automaton methods, the LBM has generally been classified into two types: one for constructing coarse-grained physical models [2, 6, 7], and the other for numerically solving governing equations [4, 8, 9]. The latter dominates the literature, so LBM typically refers to this category and is often called the standard LBM. Nevertheless, LBM can also be employed for the construction of cross-scale physical models, commonly referred to as the DBM. These two types differ in their objectives and construction principles, operate in complementary dimensions, and are each well justified in their respective contexts [5]. Unlike the standard LBM, which primarily serves as a numerical tool, the DBM is grounded in nonequilibrium statistical physics and is capable of capturing both hydrodynamic and thermodynamic nonequilibrium effects. It ensures physical consistency by constructing discrete distribution functions based on kinetic theory and decouples the time step, spatial resolution, and discrete velocity sets, offering greater flexibility and accuracy in modeling complex fluid phenomena. These factors have enabled the DBM to be widely applied over the past decade to multiphase flows [10, 11], combustion [12, 13], plasma [14, 15], Riemann problems [16], sound wave propagation [17], and various fluid instabilities [18–23].
Conventionally, the one-dimensional (1D) formulation is used for 1D problems, the two-dimensional (2D) formulation can be employed for both 1D and 2D configurations, and the three-dimensional (3D) formulation is applicable to 1D–3D cases. A natural question then arises: can the 1D model also be extended to simulate higher-dimensional systems? The answer is affirmative. In this context, we propose a preliminary cross-dimensional DBM framework capable of simulating 1D problems and readily extending to 2D and 3D cases. To achieve this, the operator splitting method [24] is adopted, which decomposes a complex evolution problem into a sequence of simpler subproblems solved successively to approximate the full solution. Splitting schemes can be first-order, such as Godunov splitting [25], or second-order, such as Strang splitting [26], with higher-order errors analyzable via the Taylor-Lie series [27]. Owing to its simplicity and modularity, operator splitting has been applied in mesoscopic simulations [28–31]. Subsequently, we employ first-order operator splitting to decompose the 3D Boltzmann equation into a sequence of 1D evolution processes, each advancing the system along a single spatial direction. The remainder of this paper is organized as follows: Part 2 introduces the construction of the cross-dimensional DBM, Part 3 presents its verification through the Sod shock tube, Lax shock tube, 2D Riemann problem, translational motion, and sound wave tests, and Part 4 summarizes the main conclusions.
2. Discrete Boltzmann method
The Bhatnagar–Gross–Krook discrete Boltzmann equations take the form:
$\begin{align} \frac{{\partial {{f}_i}}}{{\partial t}} + {\boldsymbol{v}_i}\cdot{\nabla}{{f}_i} = -\frac{1}{\tau}\left({{f}_i}-{{f}}^{\,\,{\textrm{eq}}}_i\right) .\end{align}$
Here, $t$ denotes time, $\boldsymbol{v}_i$ represents the discrete velocity vector, and $f_i$ and $f_i^{\,\,{\mathrm{eq}}}$ are the discrete distribution function and its equilibrium counterpart, respectively. The relaxation time $\tau$ determines the rate at which $f_i$ approaches $f_i^{\,\,{\mathrm{eq}}}$.The principle of physical consistency requires the following relationship,
$\begin{align} \iint f^{\,\,{\textrm{eq}}} \Psi \mathrm{d}\boldsymbol{v}\mathrm{d}{\eta} = \sum\limits_{i}f^{\,\,{\textrm{eq}}}_{i}{\Psi}_i ,\end{align}$
where $\Psi = 1, \boldsymbol{v}, \boldsymbol{v}\cdot\boldsymbol{v}+{\eta}^2, \boldsymbol{v}\boldsymbol{v}, (\boldsymbol{v}\cdot\boldsymbol{v}+{\eta}^2)\boldsymbol{v}, \dots$, and ${\Psi}_i = 1, \boldsymbol{v}_i, \boldsymbol{v}_i\cdot\boldsymbol{v}_i+{\eta}_i^2, \boldsymbol{v}_i\boldsymbol{v}_i, (\boldsymbol{v}_i\cdot\boldsymbol{v}_i+{\eta}_i^2)\boldsymbol{v}_i, \dots$. Here, $\boldsymbol{v}$ and $\boldsymbol{v}_i$ are employed to describe the translational energy, while $\eta$ and $\eta_i$ are associated with vibrational and/or rotational energies. The values of $\boldsymbol{v}_i$ and $\eta_i$ are tunable to enhance the robustness and accuracy of the DBM. Specifically, the values of discrete velocities $\boldsymbol{v}_i$ are typically chosen around the flow velocity $|\boldsymbol{u}|$ and sound speed $v_\mathrm s = \sqrt{\gamma T}$, where $\gamma = {(D+I+2)}/{(D+I)}$, $T$ denotes the temperature, $D$ is the spatial dimension, and $I$ represents the number of extra degrees of freedom. Meanwhile, the parameters $\eta_i$ are centered near $\bar{\eta} = \sqrt{{IT}/{m}}$, with the molar mass $m = 1$, in accordance with the equipartition theorem of energy. The equilibrium distribution function is expressed by $\begin{align} {f}^{\,\,{\textrm{eq}}} = \rho\left(\dfrac{1}{2\pi T}\right)^{D/2}\left(\dfrac{1}{2\pi I T}\right)^{1/2}\mathrm{\exp}\left(-\dfrac{|\boldsymbol{v}-\boldsymbol{u}|^2}{2T}-\dfrac{{\eta}^2}{2IT}\right) ,\end{align}$
where $\rho$ is the density.In addition, the mass, momentum, and energy are associated with the first three groups of equation (2 ), where $f^{\,\,{\textrm{eq}}}$ and $f^{\,\,{\textrm{eq}}}_{i}$ can be replaced by $f$ and $f_{i}$, respectively. To be specific,
$\begin{align} &\qquad\qquad\quad \rho = \sum\limits_{i}f_{i} ,\end{align}$
$\begin{align} &\qquad\quad\quad \boldsymbol{J} = \sum\limits_{i}f_{i}{\boldsymbol{v}}_{i} = \rho \boldsymbol{u} ,\end{align}$
$\begin{align} & E = \frac{1}{2}\sum\limits_{i}f_{i}\left(|\boldsymbol{v}_{i}|^2+{\eta^2_i}\right) = \frac{1}{2}\left(D + I\right)\rho T +\frac{1}{2 \rho}{|\boldsymbol{J}|}^2 .\end{align}$
Accordingly, $\boldsymbol{u}$ and $T$ can be obtained from the conserved moments as $\begin{align} &\qquad \boldsymbol{u} = \frac{\boldsymbol{J}}{\rho} ,\end{align}$
$\begin{align} & T = \frac{2E - \rho|\boldsymbol{u}|^{2}}{\left(D+I\right)\rho} .\end{align}$
Furthermore, the Chapman–Enskog multiscale analysis demonstrates that the Euler equations can be recovered in the continuum limit if the first five elements of $\Psi$ and ${\Psi}_i$ are satisfied in equation (2 ). Specifically, the Euler equations read
$\begin{eqnarray} \frac{\partial \rho }{\partial t}+{\nabla}\cdot\left({\rho\boldsymbol{u}}\right) = 0 ,\end{eqnarray}$
$\begin{eqnarray} \frac{\partial \left({\rho\boldsymbol{u}}\right)}{\partial t}+{\nabla}\cdot\left({\rho\boldsymbol{u}}\boldsymbol{u}+p \boldsymbol{I}\right) = 0 ,\end{eqnarray}$
$\begin{eqnarray} \frac{\partial E }{\partial t}+{\nabla}\cdot\left({E \boldsymbol{u}}+p \boldsymbol{u}\right) = 0 ,\end{eqnarray}$
where $p = \rho T$ denotes the pressure, and $\boldsymbol{I}$ is an identity matrix.At the Euler level, the 1D, 2D, and 3D models require five, nine, and fourteen discrete velocities, respectively [32, 33]. The motivation of this study is to demonstrate the idea that a 1D model is capable of simulating not only 1D physical systems but also 2D or 3D ones. For this sake, we construct the discrete velocity model D1V5 incorporating additional degrees of freedom, as illustrated in figure 1. The discrete parameter sets are given by
$\begin{align} \left\{ \begin{aligned} \left(v_1, v_2, v_3, v_4, v_5\right) & = \left(0, v_a, -v_a, v_b, -v_b\right), \\ \left(\eta_1, \eta_2, \eta_3, \eta_4, \eta_5\right) & = \left(\eta_a, \eta_b, \eta_b, \eta_c, \eta_c\right). \end{aligned} \right.\end{align}$
The spatial derivatives are discretized using a second-order non-oscillatory, parameter-free dissipative scheme [34], while the temporal evolution is integrated with a first-order forward Euler scheme.
Figure 1. Discrete velocity set. |
Let us introduce the flowchart of the 1D DBM applied to a 3D fluid system, as shown in figure 2. At the initial time, the physical field is specified by the density $\rho_0$, flow velocity $\boldsymbol{u}_0 = (u_{x0},u_{y0},u_{z0})$, and temperature $T_0$. These quantities serve as the initial input to the flowchart’s initialization module. The evolution proceeds through three sequential steps as follows.
Figure 2. Flowchart of the cross-dimensional DBM. |
Step 1: Evolution in the $x$ direction (${{f}_{i}}\xrightarrow{x}f_{i}^{\,{*}}$)
First of all, the physical field evolves along the $x$-axis. The update procedure includes the following substeps: (i) apply the $x$-direction boundary conditions; (ii) evaluate the local equilibrium distribution $f_i^{\,\,{\textrm{eq}}} = f_i^{\,\,{\textrm{eq}}}(\rho, u_x, T)$, and assign $f_i^{\,\,{\textrm{eq}}}$ to $f_i$, i.e. $f_i^{\,\,{\textrm{eq}}} \rightarrow f_i$; (iii) advance the discrete distribution functions from $f_i$ to $f^{\,{*}}_i$ according to
$\begin{align} \frac{\partial f_i}{\partial t} + v_{i}\frac{\partial f_i}{\partial x} = -\frac{1}{\tau}\left(f_i - f_i^{\,\,{\textrm{eq}}}\right),\end{align}$
and (iv) update the physical variables from the discrete distribution functions $f^{\,{*}}_i$, i.e. $\rho = \rho(f^{\,{*}}_i)$, $u_{x} = u_{x}(f^{\,{*}}_i)$, and $T = T(f^{\,{*}}_i)$.Step 2: Evolution in the $y$ direction ($f_{i}^{\,{*}}\xrightarrow{y}f_{i}^{\,{**}}$)
Afterwards, the field evolves along the $y$-axis, following the same sequence: (i) impose the $y$-direction boundary conditions; (ii) with updated $\rho(f^{\,{*}}_i)$ and $T(f^{\,{*}}_i)$, compute $f_i^{\,{\textrm{eq}*}} = f_i^{\,\,{\textrm{eq}}}(\rho, u_y, T)$, and set $f_i^{\,{\textrm{eq}*}}$ to $f_i^*$, i.e. $f_i^{\,{\textrm{eq}*}} \rightarrow f_i^*$; (iii) update the discrete distribution functions from $f^{\,{*}}_i$ to $f^{\,{**}}_i$ through
$\begin{align} \frac{\partial f^{\,{*}}_i}{\partial t} + v_{i}\frac{\partial f^{\,{*}}_i}{\partial y} = -\frac{1}{\tau}\left(f^{\,{*}}_i - f_i^{\,{\textrm{eq}*}}\right),\end{align}$
and (iv) update the macroscopic quantities $\rho = \rho(f^{\,{**}}_i)$, $u_{y} = u_{y}(f^{\,{**}}_i)$, and $T = T(f^{\,{**}}_i)$.Step 3: Evolution in the $z$ direction (${{f}^{\,{**}}_{i}}\xrightarrow{z}f_{i}^{\,{***}}$)
Finally, the field evolves along the $z$-axis: (i) apply the $z$-direction boundary conditions; (ii) with renewed $\rho(f^{\,{**}}_i)$ and $T(f^{\,{**}}_i)$, compute $f_i^{\,{\textrm{eq}**}} = f_i^{\,\,{\textrm{eq}}}(\rho, u_z, T)$, and assign $f_i^{\,{\textrm{eq}**}}$ to $f^{\,{**}}_i$, i.e. $f_i^{\,{\textrm{eq}**}}\rightarrow f^{\,{**}}_i$; (iii) advance the discrete distribution functions from $f^{\,{**}}_i$ to $f^{\,{***}}_i$ following
$\begin{align} \frac{\partial f^{\,{**}}_i}{\partial t} + v_{i}\frac{\partial f^{\,{**}}_i}{\partial z} = -\frac{1}{\tau}\left(f^{\,{**}}_i - f_i^{\,{\textrm{eq}**}}\right),\end{align}$
and (iv) update the macroscopic variables $\rho = \rho(f^{\,{***}}_i)$, $u_{z} = u_{z}(f^{\,{***}}_i)$, and $T = T(f^{\,{***}}_i)$.At this stage, one full evolution cycle is completed, and the distribution functions have been advanced from their previous values $f_i(t)$ to the updated values $f_i(t+\Delta t) = f^{\,{***}}_i$. Subsequently, the physical field continues to evolve following the same three-step procedure. This operator splitting method provides a natural framework using a 1D model to simulate a cross-dimensional system, where 1D and 2D cases can be readily obtained by omitting the later steps.
3. Numerical validation
To validate the reliability of the cross-dimensional model, five classical benchmarks are examined: the Sod shock tube, Lax shock tube, 2D Riemann problem, translational motion, and sound wave.
3.1. Sod shock tube
First, the 1D Sod shock tube problem is simulated, with the initial conditions specified as
$\begin{align} \left(\rho_\mathrm{L}, u_\mathrm{L}, T_\mathrm{L}\right) = \left(1, 0, 1\right),\,\,\,\,\,\, \left(\rho_\mathrm{R}, u_\mathrm{R}, T_\mathrm{R}\right) = \left(0.125, 0, 0.8\right) ,\end{align}$
where the subscripts ‘$\mathrm{L}$’ and ‘$\mathrm{R}$’ denote the left and right states separated by the initial discontinuity located at $x_0 = 0.5$. The simulations are carried out on a grid with $N_x = 5000$ points, using a spatial step $\Delta x = 2 \times 10^{-4}$ and a time step $\Delta t = 5 \times 10^{-6}$. The specific heat ratio is given by $\gamma = 1.4$, and the discrete parameters are set as $(v_a, v_b, \eta_a, \eta_b, \eta_c) = (1, 5, 3.2, 0, 0)$. Free inflow and outflow boundary conditions are imposed at the left and right boundaries, respectively.Figures 3(a)–(d) illustrate the profiles of density, pressure, velocity, and temperature at $t = 0.25$. The Mach number of the shock wave is $\textrm{M}_a = 1.66$, calculated as the ratio of the shock propagation speed to the upstream sound speed. Therefore, the flow in this configuration is compressible, as the Mach number exceeds the commonly used threshold of 0.3. Furthermore, symbols represent the DBM results, while solid lines correspond to the Riemann solutions. The numerical solution accurately resolves the characteristic wave structures, including a left-propagating rarefaction wave, a contact discontinuity in the middle, and a right-propagating shock wave. Excellent agreement is observed between the numerical and Riemann solutions, demonstrating that the proposed model is capable of accurately simulating compressible flows.
Figure 3. Sod shock tube field at $t = 0.25$: (a) density, (b) pressure, (c) velocity, and (d) temperature. |
Moreover, a quantitative comparison of the computational efficiency of the D1V5 and D2V9 models [33] is performed. The computational time of the D1V5 is 9 s, while that of the D2V9 is 36.33 s. The computational time of the D2V9 is about 4 times that of the D1V5. This significant efficiency difference is attributed to the reduced number of discrete velocities and the lower spatial dimensionality of the evolution process: the D1V5 evolves the distribution function along a single spatial direction, whereas the D2V9 requires evolution in two spatial dimensions, resulting in a higher number of computational operations per iteration.
3.2. Lax shock tube
We next consider another 1D Riemann problem, the Lax shock tube, with the initial states given by
$\begin{align} \left(\rho_\mathrm{L}, u_\mathrm{L}, T_\mathrm{L}\right) & = \left(0.445, 0.698, 7.928\right),\nonumber\\ \left(\rho_\mathrm{R}, u_\mathrm{R}, T_\mathrm{R}\right)& = \left(0.5, 0, 1.142\right) .\end{align}$
All other conditions are the same as those used in the Sod shock tube problem. Figures 4(a)–(d) compare the DBM results with the Riemann solutions for density, pressure, velocity, and temperature at $t = 0.15$, showing excellent agreement. In the rarefaction region, the density, pressure, and temperature decrease smoothly from left to right, while the velocity increases continuously. In the vicinity of the material interface, the density and temperature exhibit opposite variations, whereas the pressure and flow velocity remain constant. Meanwhile, all variables exhibit sharp jumps across the shock front which travels rightwards with a supersonic speed.
Figure 4. Lax shock tube field at $t = 0.15$: (a) density, (b) pressure, (c) velocity, and (d) temperature. |
3.3. 2D Riemann problem
This subsection utilizes the 2D Riemann problem in gas dynamics to compare the D1V5 and D2V9 models [33]. The initial physical field is divided into four adjacent rectangular subdomains, each of which is assigned a uniform state for density, pressure, and flow velocity. Two representative configurations are selected for comparative analysis. For configuration I,
$\begin{align} \left(\rho,p,u_x,u_y\right) = \left\{ \begin{aligned} &\left(0.8,1,0,0\right), 0 \lt x\unicode{x2A7D}0.1, 0 \lt y\unicode{x2A7D}0.1,\\ &\left(1,1,0.05,0\right), 0.1 \lt x\unicode{x2A7D}0.2, 0 \lt y\unicode{x2A7D}0.1,\\ &\left(1,1,0,0.05\right), 0 \lt x\unicode{x2A7D}0.1, 0.1 \lt y\unicode{x2A7D}0.2,\\ &\left(0.5,0.6,0,0\right), 0.1 \lt x\unicode{x2A7D}0.2, 0.1 \lt y\unicode{x2A7D}0.2.\\ \end{aligned} \right.\end{align}$
For configuration II, $\begin{align} \left(\rho,p,u_x,u_y\right) = \left\{ \begin{aligned} &\left(1.1,1.1,0.05,0.05\right), 0 \lt x\unicode{x2A7D}0.1, 0 \lt y\unicode{x2A7D}0.1,\\ &\left(0.3,0.35,0.05,0\right), 0.1 \lt x\unicode{x2A7D}0.2, 0 \lt y\unicode{x2A7D}0.1,\\ &\left(0.3,0.35,0,0.05\right), 0 \lt x\unicode{x2A7D}0.1, 0.1 \lt y\unicode{x2A7D}0.2,\\ &\left(1.1,1.1,0,0\right), 0.1 \lt x\unicode{x2A7D}0.2, 0.1 \lt y\unicode{x2A7D}0.2.\\ \end{aligned} \right.\end{align}$
The simulation employs a grid of $N_x = N_y = 400$, with inflow/outflow boundary conditions applied in all directions. The grid and time steps are set to $\Delta x = \Delta y = 5 \times 10^{-4}$ and $\Delta t = 1 \times 10^{-4}$, respectively.Figures 5 and 6 display the density fields for configurations I and II, respectively. Each figure shows the snapshots obtained using the D1V5 and D2V9 models at time instants $t = 0$, $0.02$, and $0.04$. For configuration I, the discrete parameters are set as $(v_a, v_b, \eta_a, \eta_b, \eta_c) = (1, 0.9, 3.8, 0, 0)$ for D1V5, and $(v_a, v_b, v_c, \eta_a, \eta_b, \eta_c) = (1, 1.9, 1, 3.8, 0, 0)$ for D2V9; For configuration II, $(v_a, v_b, \eta_a, \eta_b, \eta_c) = (1, 3.2, 1.8, 0, 0)$ for D1V5, and $(v_a, v_b, v_c, \eta_a, \eta_b, \eta_c) = (1, 3.2, 1, 1.8, 0, 0)$ for D2V9. At the interfaces between the four subdomains, various waves are generated and interact with one another, leading to the formation of complex flow structures near the central region. Meanwhile, the complex patterns simulated by the D1V5 and D2V9 models are similar to each other.
Figure 5. Density contours of the 2D Riemann problem for configuration I at various time instants. Top row: D1V5; Bottom row: D2V9. |
Figure 6. Density contours of the 2D Riemann problem for configuration II at various time instants. Top row: D1V5; Bottom row: D2V9. |
Taking configuration I as an example, a quantitative comparison of the computational efficiency of the two models is performed. The computational time of the D1V5 is 29.33 s, while that of the D2V9 is 24 s. These results indicate that, for this configuration, the computational efficiency of the D2V9 is approximately 22% higher than that of the D1V5. This difference in efficiency results from the difference in algorithmic complexity between the two models.
3.4. Translational motion
Subsequently, Galilean invariance of the proposed model is examined through a translational motion test. A circular fluid region is initially placed at the center of the square computational domain, located at $[L_x/2, L_y/2]$, with a radius $R = L_x/4$, where $L_x = N_x \Delta x$ and $L_y = N_y \Delta y$ denote the domain width and length, respectively. The inner and outer densities are set to $\rho_{\mathrm{in}} = 1$ and $\rho_{\mathrm{out}} = 1.1$, respectively, with a smooth transition layer of thickness $W = L_x/50$ connecting the two regions. Mathematically, the initial density field is prescribed as
$\begin{align} \rho\left(y\right) = \frac{\rho_\mathrm{in}+\rho_\mathrm{out}}{2}-\frac{\rho_\mathrm{in}-\rho_\mathrm{out}}{2}\tanh\left(\frac{\sqrt{\left(x-\frac{L_x}{2}\right)^2+\left(y-\frac{L_y}{2}\right)^2}-R}{W}\right) ,\end{align}$
while the pressure is initially $p = 1$. The remaining parameters are specified as $N_x = N_y = 1000$, spatial steps $\Delta x = \Delta y = 2 \times 10^{-4}$, time step $\Delta t = 1 \times 10^{-5}$, specific heat ratio $\gamma = 1.5$, discrete parameters $(v_a, v_b, \eta_a, \eta_b, \eta_c) = (1, 5, 3.2, 0, 0)$. Periodic boundary conditions are applied in both spatial directions.A translational velocity is imposed in the diagonal direction, with $(u_x, u_y) = (0.5, 0.5)$. Snapshots of the density field at four representative times, $t = 0.1$, $0.2$, $0.3$, and $0.4$ from left to right, are shown in figure 7. In each snapshot, the density increases from bule to red, and the flow direction is indicated by white arrows. The system is expected to translate by a distance $L_x$ at $t = 0.4$, according to the theoretical relation $u_x = L_x/t$.
Figure 7. Diagonal motion of fluid at $t = 0.1$, $0.2$, $0.3$, and $0.4$. |
3.5. Sound wave
Finally, the classical sound wave is employed for verification. The initial fields in 1D, 2D, and 3D cases are set to a uniform density $\rho_0 = 1$, flow velocity $\boldsymbol{u}_0 = 0$, and temperature $T_0 = 0.5$. The computational domain is discretized using $N_{\alpha} = 1000$ grids in each coordinate direction, with a spatial resolution of $\Delta \alpha = 1 \times 10^{-4}$ and a time step of $\Delta t = 1 \times 10^{-5}$. The specific heat ratio is $\gamma = 5/3$, and the discrete parameters are set as $(v_a, v_b, \eta_a, \eta_b, \eta_c) = (1, 5, 3.2, 0, 0)$. Periodic boundary conditions are applied on all boundaries. A small perturbation is imposed at the center of each computational domain. In the 1D configuration, the perturbation develops into a planar wave that propagates bidirectionally along the tube. In 2D, the disturbance expands radially as a circular wave, while in 3D it evolves into a spherical wave propagating outward from the domain center.
Figures 8(a)–(c) present the sound wave pressure distributions in 1D, 2D, and 3D, respectively, at a time instant $t = 0.065$. Panel (a) shows the 1D pressure profile along the $x$-axis within the region $[-0.05, 0.05]$, where two peaks appear at approximately $x \approx \pm 0.04$ and the surrounding pressure remains close to the initial value of $0.5$. Panels (b) and (c) depict the 2D and 3D pressure fields over the domains $[-0.05, 0.05]^2$ and $[0, 0.05]^3$, respectively, where the perturbation has propagated outward from the center and already re-entered the computational domain. Figures 9(a)–(c) display the position of the sound wave over time in 1D, 2D, and 3D, respectively. The simulation results, obtained using the wavefront location as the reference, are consistent with the theoretical prediction $x = x_0 + v_\mathrm s t$, where the sound speed is given by $v_\mathrm s = \sqrt{\gamma T}$. These results demonstrate that the proposed model successfully captures the outward propagation of a small perturbation in 1D, 2D, and 3D simultaneously.
Figure 8. Pressure distribution of the sound wave at $t = 0.065$: (a) 1D, (b) 2D, and (c) 3D. |
Figure 9. Position of the sound wave versus time: (a) 1D, (b) 2D, and (c) 3D. |
4. Conclusions
In this work, we develop a 1D DBM capable of simulating not only 1D systems but also 2D and 3D configurations with an operator splitting method. Specifically, the model advances the physical fields sequentially along individual spatial directions, enabling a natural transition from one to higher-dimensional systems within a unified framework. Its performance is evaluated by the Sod shock tube, Lax shock tube, 2D Riemann problem, translational motion, and sound waves, and the results demonstrate the feasibility of the proposed approach for cross-dimensional simulations. One shortcoming of the method is the lack of the viscous, heat conductive, and other thermodynamic nonequilibrium effects when simulating 2D or 3D systems, thereby opening a promising avenue for future developments in mesoscopic fluid dynamics.