Computational Approaches to Compressible Micropolar Fluid Flow in Moving Parallel Plate Configurations

Abstract

In this paper, we consider the unsteady flow of a compressible micropolar fluid between two moving, thermally isolated parallel plates. The fluid is characterized as viscous and thermally conductive, with polytropic thermodynamic properties. Although the mathematical model is inherently three-dimensional, we assume that the variables depend on only a single spatial dimension, reducing the problem to a one-dimensional formulation. The non-homogeneous boundary conditions representing the movement of the plates lead to moving domain boundaries. The model is formulated in mass Lagrangian coordinates, which leads to a time-invariant domain. This work focuses on numerical simulations of the fluid flow for different configurations. Two computational approaches are used and compared. The first is based on the finite difference method and the second is based on the Faedo–Galerkin method. To apply the Faedo–Galerkin method, the boundary conditions must first be homogenized and the model equations reformulated. On the other hand, in the finite difference method, the non-homogeneous boundary conditions are implemented directly, which reduces the computational complexity of the numerical scheme. In the performed numerical experiments, it was observed that, for the same accuracy, the Faedo–Galerkin method was approximately 40 times more computationally expensive compared to the finite difference method. However, on a dense numerical grid, the finite difference method required a very small time step, which could lead to an accumulation of round-off errors. On the other hand, the Faedo–Galerkin method showed the convergence of the solutions as the number of expansion terms increased, despite the higher computational cost. Comparisons of the obtained results show good agreement between the two approaches, which confirms the consistency and validity of the numerical solutions.

1. Introduction

Micropolar fluid theory extends classical fluid mechanics by incorporating the microstructure and microrotation of fluid particles, making it suitable for modeling complex fluids such as polymers, blood, liquid crystals, and certain lubricants [1]. In the case of compressible fluids, considered in this work, the application can refer to the flow of gaseous lubricants, explosive mixtures, aerosols, and similar substances. These fluids exhibit non-Newtonian behavior, characterized by microrotational effects that classical fluid models cannot capture. The micropolar fluid model may also be more adequate for modeling flows at the micro- and nanoscale, as shown, for example, in [2], where the numerical model based on micropolar fluid theory predicts the experimental data better than the Navier–Stokes model.

The study of micropolar fluid flow between two moving parallel plates is motivated by the need to better understand and predict the behavior of these complex fluids in practical applications. The model under consideration takes into account microrotations, coupling stresses and microrotational viscosity, allowing the accurate modeling of fluids with significant microstructural characteristics. The model also includes the energy equation, which accounts for heat conduction within the fluid. The inclusion of the energy equation enables thermodynamic monitoring, which could be crucial in applications where heat transfer impacts fluid behavior.

The development of the micropolar fluid model began with the work of Eringen [3], while the theoretical analysis of compressible viscous heat-conducting micropolar fluid flow models involving temperature began with various one-dimensional problems studied by Mujaković, such as in [4,5]. The numerical approximations of different one-dimensional problems were performed in [6,7,8] using the finite difference numerical method. In these works, it was proved that the proposed scheme converges and also that a global solution to the considered problems, obtained as the limit of the approximate solutions, does exist. In most of the mentioned works, the case with homogeneous boundary conditions was considered. The problems with non-homogeneous boundary conditions have been studied, for example, in [5]. In [8], the finite difference method was successfully applied and used as a tool to prove the existence of a solution for the free boundary problem.

The flow of a compressible micropolar fluid between two parallel plates was derived and investigated in [9]. In that work, a shear flow model was developed with boundary conditions that modeled a fixed lower plate in the plane $x=0$ and an upper plate moving irrotationally at a constant distance in the plane $x=h$. The numerical simulations of the flow were performed using numerical approximations based on the Faedo–Galerkin method. Similar problems, but for classical fluids, were studied, for example, in [10,11,12].

In this paper, a similar problem is considered, but with the key assumption that the distance between the plates is not fixed but varies over time as the plates move. This type of problem is motivated by the non-homogeneous 1D problem for classical fluids described in [13], where it is called the piston problem. Such a model is obtained by assuming non-homogeneous boundary conditions for the velocity, which leads to moving boundaries and a domain that changes with time. These boundaries can represent two moving pistons, where the piston velocities coincide with the gas velocities on both sides. A similar model with two moving pistons for a classical fluid is studied, for example, in [14], where the control problem is considered. In that work, however, the piston motion is modeled using an ordinary differential equation expressing Newton’s second law. Since, in this case, the motion of the pistons is not known in advance, this type of problem is a free piston problem. Here, however, we consider the classical piston problem, where the flow is generated by the known motion of the piston defined by the given boundary velocities.

There are not many studies that deal with problems involving non-homogeneous boundary conditions and moving domain boundaries, and those that do focus mainly on classical fluids. For example, in [15], the existence of global solutions for similar one-dimensional problems with classical fluids in time-decreasing non-rectangular domains is considered, while in [16], a free piston problem with non-homogeneous boundary conditions is studied. In contrast to these works, the present study is, to our knowledge, the first to investigate a shear flow problem with a micropolar fluid and non-homogeneous boundary conditions defining moving domain boundaries.

This study contributes to the understanding of the dynamics of micropolar fluids by addressing a previously unexplored problem: the flow of compressible micropolar fluids between two moving plates that define a dynamically evolving model domain. The novelty of this work lies in the use of mass Lagrangian coordinates, which reformulate the variable-domain problem into a fixed-domain model, allowing the application of efficient numerical methods. The second contribution focuses on the definition of robust numerical methods to solve the problem, namely the finite difference method and the Faedo–Galerkin method. To ensure that the methods are both stable and convergent, great attention is paid to the discretization process. In the case of the finite difference method, a staggered grid is used to preserve the physical relationships between the approximate variables to prevent potential instabilities. On the other hand, for the Faedo–Galerkin method, an appropriate substitution is introduced to reformulate the problem into a problem with homogeneous boundary conditions so that the standard Faedo–Galerkin approach can be applied. Finally, the performance of both methods was validated by comparing their results in selected test experiments with moving plates. The strong agreement and convergence of the results show that the proposed methods are stable and applicable to problems of this type.

To the best of our knowledge, this study is the first to address the shear flow of micropolar fluids with non-homogeneous boundary conditions in a time-dependent domain. By extending classical fluid mechanics to account for microrotations, couple stresses, and heat conduction, this work provides a novel framework for the simulation of complex fluids in dynamically changing geometries. These findings could be relevant for applications involving dynamically changing domains, such as lubricant flow in moving mechanical systems, microfluidics in expanding or contracting channels, or fluid transport in moving micro-/nanostructures.

2. Mathematical Model

We consider a mathematical model for the flow of the compressible viscous micropolar fluid between two parallel plates that was developed and described in [9] and which has the following formulation:

$${\partial }_t\rho +u{\partial }_x\rho +\rho {\partial }_{x}u=0,$$

(1)

$$\rho ({\partial }_{t}u+u{\partial }_{x}u)=-{\partial }_x\left(R\rho \theta \right)+(\lambda +2\mu ){\partial }_{xx}u,$$

(2)

$$\rho ({\partial }_t{v}_2+u{\partial }_x{v}_2)=(\mu +{\mu }_r){\partial }_{xx}{v}_2-2{\mu }_r{\partial }_x{\omega }_3,$$

(3)

$$\rho ({\partial }_t{v}_3+u{\partial }_x{v}_3)=(\mu +{\mu }_r){\partial }_{xx}{v}_3+2{\mu }_r{\partial }_x{\omega }_2,$$

(4)

$$j_I\rho ({\partial }_{t}w+u{\partial }_{x}w)=(c_0+2c_d){\partial }_{xx}w-4{\mu }_{r}w,$$

(5)

$$j_I\rho ({\partial }_t{\omega }_2+u{\partial }_x{\omega }_2)=(c_d+c_a){\partial }_{xx}{\omega }_2-4{\mu }_r{\omega }_2-2{\mu }_r{\partial }_x{v}_3,$$

(6)

$$j_I\rho ({\partial }_t{\omega }_3+u{\partial }_x{\omega }_3)=(c_d+c_a){\partial }_{xx}{\omega }_3-4{\mu }_r{\omega }_3+2{\mu }_r{\partial }_x{v}_2,$$

(7)

$$\begin{array}{c}\hfill c_v\rho ({\partial }_t\theta +u{\partial }_x\theta )=k_{\theta }{\partial }_{xx}\theta -R\rho \theta {\partial }_{x}u+(\lambda +2\mu ){\left({\partial }_{x}u\right)}^2+(\mu +{\mu }_r)\left({\left({\partial }_x{v}_2\right)}^2+{\left({\partial }_x{v}_3\right)}^2\right)\\ \hfill +4{\mu }_r{\left|\mathit{\omega }\right|}^2+4{\mu }_r({\omega }_2{\partial }_x{v}_3-{\omega }_3{\partial }_x{v}_2)+(c_0+2c_d){\left({\partial }_{x}w\right)}^2+(c_d+c_a)\left({\left({\partial }_x{\omega }_2\right)}^2+{\left({\partial }_x{\omega }_3\right)}^2\right)\end{array}$$

(8)

in the domain $\Omega \times ]0,T[$, for $T>0$ and $\Omega \subset \mathbf{R}$.

Here, $\rho$, $\mathbf{v}=(u,v_1,v_2)$, $\mathit{\omega }=(w,{\omega }_1,{\omega }_2)$, and $\theta$ are variables in the system and denote the mass density, the velocity vector, the microrotation velocity vector, and the absolute temperature, respectively. Equations (1)–(8) are formulated in the Eulerian framework and represent the conservation laws of mass, momentum, angular momentum, and energy. To obtain the system in the presented form, it is assumed that the flow is such that all functions entering into the equations are independent of the horizontal spatial coordinates y and z and depend only on x (see Figure 1), i.e.,

$$\rho (\mathbf{x},t)=\rho (x,t),\theta (\mathbf{x},t)=\theta (x,t),$$

(9)

$$\mathbf{v}(\mathbf{x},t)=\mathbf{v}(x,t)=\left(u(x,t),v_2(x,t),v_3(x,t)\right),$$

(10)

$$\mathit{\omega }(\mathbf{x},t)=\mathit{\omega }(x,t)=\left(w(x,t),{\omega }_2(x,t),{\omega }_3(x,t)\right).$$

(11)

The model (1)–(8) incorporates the constitutive equations for a micropolar continuum, which define the stress tensor and the couple stress tensor. These relations take into account the interactions between stresses, microrotations, and deformation rates, as detailed in [9]. The internal energy density E is assumed to be proportional to the temperature $\theta$, specifically $E=c_v\theta$, where $c_v>0$ is the specific heat at a constant volume. Accordingly, the corresponding energy in Equation (8) is expressed in terms of $\theta$. Additionally, Equation (8) incorporates Fourier’s law, which states that the heat flux density is given by $-k_{\theta }{\partial }_x\theta$, where $k_{\theta }>0$ is the thermal conductivity. It is also assumed that the fluid is perfect and polytropic, which leads to the pressure being defined as $p=R\rho \theta$, where $R>0$ is the specific gas constant. The other symbols in the previous expressions are as follows: $j_I$ stands for the microinertial density; $\lambda$ and $\mu$ denote the viscosity coefficients; and ${\mu }_r$, $c_0$, $c_d$, and $c_a$ are the microviscosity coefficients. These coefficients must satisfy the Clausius–Duhamel inequalities to ensure that the principles of thermodynamics are maintained (see, for example, [9]).

To formulate the specific type of flow, the initial and boundary conditions should be added to the system (1)–(8). Let us assume that the gas occupies the domain $]0,h[$ at the initial time $t=0$ and that there are moving pistons at both boundaries. The position of the left piston at the time t is given by $h_L\left(t\right)$ and the position of the right piston is given by $h_R\left(t\right)$. We assume, furthermore, that $h_R\left(t\right)-h_L\left(t\right)\ge \delta >0$, which ensures that the distance between the pistons remains positive at all times. It is also assumed that no gas passes through the pistons and that the velocities of the pistons on both sides coincide with the fluid velocities, i.e.,

$$u(h_L\left(t\right),t)=u_L\left(t\right)={\displaystyle \frac{d{h}_L}{dt}},u(h_R\left(t\right),t)=u_R\left(t\right)={\displaystyle \frac{d{h}_R}{dt}}.$$

(12)

It follows that the domain has moving boundaries and, at the time t, is equal to $[h_L\left(t\right),h_R\left(t\right)]$. We assume that the velocity vectors at the boundaries are defined by

$$\mathbf{v}(h_L\left(t\right),t)={\mathbf{v}}_L\left(t\right)=(u_L\left(t\right),v_{2L}\left(t\right),v_{3L}\left(t\right)),\mathbf{v}(h_R\left(t\right),t)={\mathbf{v}}_R\left(t\right)=(u_R\left(t\right),v_{2R}\left(t\right),v_{3R}\left(t\right)),$$

(13)

where the first components of the velocity vectors govern the dynamics of the moving boundaries. If the velocity is prescribed, the positions of the boundaries can be determined by solving (12), with the initial conditions $h_L\left(0\right)=0$ and $h_R\left(0\right)=h$. The described flow configuration is schematically presented in Figure 1. We assume the following homogenous boundary conditions for the microrotation field and the heat flux:

$$\mathit{\omega }(h_L\left(t\right),t)=\mathit{\omega }(h_R\left(t\right),t)=\mathbf{0},{\partial }_x\theta (h_L\left(t\right),t)={\partial }_x\theta (h_R\left(t\right),t)=0.$$

(14)

For the initial conditions, we take

$$\rho (x,0)={\rho }_0\left(x\right),\theta (x,0)={\theta }_0\left(x\right)$$

(15)

$$\mathbf{v}(x,0)=\left(u_0\left(x\right),v_{20}\left(x\right),v_{30}\left(x\right)\right),\mathit{\omega }(x,0)=\left(w_0\left(x\right),{\omega }_{20}\left(x\right),{\omega }_{30}\left(x\right)\right),$$

(16)

for $x\in [0,h],$ where the functions on the right-hand sides are smooth enough. To define the problem appropriately, we assume that the compatibility conditions between the initial conditions on the boundary and the boundary conditions are satisfied. We also assume that the initial density and the initial temperature are strictly positive and bounded, i.e.,

$$0<m\le {\rho }_0\left(x\right)\le M,m\le {\theta }_0\left(x\right)\le M,$$

for some $m,M>0$.

As mentioned earlier, the developed model is formulated in Eulerian coordinates in a domain characterized by moving boundaries. In order to analyze such models effectively, it is common to transform them into Lagrangian coordinates. This transformation simplifies the handling of moving boundaries and improves the accuracy of numerical approximations.

Specifically, we transform the problem described by Equations (1)–(8) into mass Lagrangian coordinates. This is achieved by introducing a new variable, y, which is defined as follows:

$$y={\displaystyle \frac{1}{L}}\Psi (x,t),t\ge 0,x\in [h_L\left(t\right),h_R\left(t\right)],$$

(17)

where

$$\Psi (x,t)={\int }_{h_L\left(t\right)}^{x(y,t)}\rho (s,t)ds\mathrm{and}L=\Psi (h_R\left(t\right),0).$$

From (1) and (13), it follows that, for $t>0$,

$${\int }_{h_L\left(t\right)}^{h_R\left(t\right)}\rho (x,t)dx={\int }_0^h{\rho }_0\left(x\right)=L,$$

which means that the total mass of the gas is conserved. Note that if ${\rho }_0\left(x\right)>0,x\in [0,h]$, which is physically justified because it ensures that no vacuum occurs in the domain, then for any $t\ge 0$, $\Psi (x,·)$ is one-to-one from $[h_L\left(t\right),h_R\left(t\right)]$ onto $[0,1]$, so that ${\Psi }^{-1}(y,·)$ is well defined. By introducing new coordinates, the system becomes parabolic. Also, note that the domain in the new coordinates becomes fixed and equal to $[0,1]$, which is spatially dimensionless, which can be particularly advantageous for modeling at the micro- and nanoscales. The variables $\rho$, $\mathbf{v}$, $\mathit{\omega }$, and $\theta$ in the new coordinates become

$$\widehat{\rho }(y,t)=\rho ({\Psi }^{-1}(Ly,t),t),\widehat{\theta }(y,t)=\theta ({\Psi }^{-1}(Ly,t),t),$$

(18)

$$\widehat{\mathbf{v}}(y,t)=\mathbf{v}({\Psi }^{-1}(Ly,t),t),\widehat{\mathit{\omega }}(y,t)=\mathit{\omega }({\Psi }^{-1}(Ly,t),t),$$

(19)

for $y\in \left[0,1\right]$. The mass Lagrangian coordinates are connected to the Eulerian coordinates by

$$x(y,t)=x_0\left(y\right)+{\int }_0^t\widehat{u}(y,\tau )d\tau ,x_0\left(y\right)={\Psi }^{-1}(Ly,0).$$

(20)

Written in new coordinates, the system (1)–(8), after omitting the “hat” notation, becomes the same as in [9] and equal to

$${\partial }_t\rho +L^{-1}{\rho }^2{\partial }_{y}u=0,$$

(21)

$${\partial }_{t}u=-L^{-1}R{\partial }_y\left(\rho \theta \right)+L^{-2}(\lambda +2\mu ){\partial }_y\left(\rho {\partial }_{y}u\right),$$

(22)

$${\partial }_t\mathbf{u}=L^{-2}(\mu +{\mu }_r){\partial }_y\left(\rho {\partial }_y\mathbf{u}\right)+2L^{-1}{\mu }_r\nabla \times \mathbf{w},$$

(23)

$$j_I{\partial }_{t}w=L^{-2}(c_0+2c_d){\partial }_y\left(\rho {\partial }_{y}w\right)-4{\mu }_r{\displaystyle \frac{w}{\rho }},$$

(24)

$$j_I{\partial }_t\mathbf{w}=L^{-2}(c_d+c_a){\partial }_y\left(\rho {\partial }_y\mathbf{w}\right)-{\displaystyle \frac{4{\mu }_r}{\rho }}\mathbf{w}+2L^{-1}{\mu }_r\nabla \times \mathbf{u},$$

(25)

$$\begin{array}{c}\hfill c_v{\partial }_t\theta =L^{-2}{k}_{\theta }{\partial }_y\left(\rho {\partial }_y\theta \right)-L^{-1}R\rho \theta {\partial }_{y}u+L^{-2}(\lambda +2\mu )\rho {\left({\partial }_{y}u\right)}^2\\ \hfill +L^{-2}(\mu +{\mu }_r)\rho {\left|{\partial }_y\mathbf{u}\right|}^2+4{\mu }_r{\displaystyle \frac{w^2}{\rho }}+4{\mu }_r{\displaystyle \frac{{\left|\mathbf{w}\right|}^2}{\rho }}-4L^{-1}{\mu }_r(\nabla \times \mathbf{u})·\mathbf{w}+\\ \hfill L^{-2}(c_0+2c_d)\rho {\left({\partial }_{y}w\right)}^2+L^{-2}(c_d+c_a)\rho {\left|{\partial }_y\mathbf{w}\right|}^2,\end{array}$$

(26)

for $(y,t)\in ]0,1[\times ]0,T[$, where we denote $\mathbf{u}=(0,v_2,v_3),\mathbf{w}=(0,{\omega }_2,{\omega }_3).$ Observe that the vectors $\mathbf{u}$ and $\mathbf{w}$ are obtained from the velocity vector $\mathbf{v}$ and the microrotation velocity vector $\mathit{\omega }$ by setting their first components to zero. The boundary conditions for $t\in \left[0,T\right]$ now read as

$$u(0,t)=u_L\left(t\right),u(1,t)=u_R\left(t\right),$$

(27)

$$\mathbf{u}(0,t)={\mathbf{u}}_L\left(t\right)=(0,v_{2L}\left(t\right),v_{3L}\left(t\right)),\mathbf{u}(1,t)={\mathbf{u}}_R\left(t\right)=(0,v_{2R}\left(t\right),v_{3R}\left(t\right)),$$

(28)

$$w(0,t)=w(1,t)=0,\mathbf{w}(0,t)=\mathbf{w}(1,t)=\mathbf{0},{\partial }_y\theta (0,t)={\partial }_y\theta (1,t)=0.$$

(29)

The initial conditions are derived by transforming the initial functions, defined in (15)–(16), into the new coordinate system, that is,

$$\rho (y,0)={\rho }_0\left({\Psi }^{-1}(Ly,0)\right),\theta (y,0)={\theta }_0\left({\Psi }^{-1}(Ly,0)\right),$$

(30)

$$u(y,0)=u_0\left({\Psi }^{-1}(Ly,0)\right),\mathbf{u}(y,0)=\left(0,v_{20}\left({\Psi }^{-1}(Ly,0)\right),v_{30}\left({\Psi }^{-1}(Ly,0)\right)\right),$$

(31)

and

$$w(y,0)=w_0\left({\Psi }^{-1}(Ly,0)\right),\mathbf{w}(y,0)=\left(0,{\omega }_{20}\left({\Psi }^{-1}(Ly,0)\right),{\omega }_{30}\left({\Psi }^{-1}(Ly,0)\right)\right).$$

(32)

Note that Equation (21) could also be written in the form

$${\partial }_t\left(\frac{1}{\rho }\right)=L^{-1}{\partial }_{y}u,$$

(33)

which, after integration over $[0,1]\times [0,t]$, and using boundary conditions for u, gives

$${\displaystyle {\int }_0^1{\displaystyle \frac{dy}{\rho (y,t)}}={\int }_0^1{\displaystyle \frac{dy}{{\rho }_0\left(x\right)}}+\frac{1}{L}{\int }_0^t(u_R\left(\tau \right)-u_L\left(\tau \right))d\tau =A\left(t\right).}$$

(34)

To define the problem appropriately, we assume that there exists a constant, $\delta >0$, such that $A\left(t\right)\ge \delta$ for $t\in [0,T]$. This condition is physically justified as it means that the distance between the moving domain boundaries is greater than zero. Under this assumption, the local existence of the solution for a 1D problem with non-homogeneous boundary conditions for velocities was proved in [5]. Therefore, we expect that if this condition is satisfied, and the initial and boundary functions are smooth enough, the generalized solution to our problem exists. We are interested here in the numerical simulations of the flow modeled by the described model.

In the continuation of this paper, we provide numerical approximate solutions to this problem using two different numerical methods and show that they provide compatible solutions to the selected experiments. The proposed approximate solutions can be the basis for proving the existence of a solution to this problem. The main difference between the considered problem and the problem presented in [9] is the use of the non-homogeneous boundary conditions modeling the moving domain boundaries. The movement of the upper plate modeled in [9] was only in the $yz$ direction, so the domain was fixed.

3. Finite Difference Method

We first propose a finite difference scheme based on the semi-discrete approach, which means that the spatial discretization is performed first, while the problem remains continuous in time. A uniform staggered grid is used for the spatial discretization. The convergence of the numerical scheme of this type for similar problems was proved in [6,8]. The detailed description of the scheme follows. arrows.meta

Let h denote an increment in y such that $h={\displaystyle \frac{1}{N}}$ for the chosen number of discretization points, $N\in {\mathrm{Z}}^{+}$. The staggered grid points (Figure 2) are denoted by $y_k=kh$, $k\in \{0,1,\dots ,N\}$, and $y_j=jh$, $j\in \left\{{\displaystyle \frac{1}{2}},\dots ,N-{\displaystyle \frac{1}{2}}\right\}$.

We construct the following time-dependent functions

$${\rho }_j\left(t\right),{\theta }_j\left(t\right),j=\frac{1}{2},\dots ,N-\frac{1}{2},$$

(35)

$$u_k\left(t\right),{\mathbf{u}}_k\left(t\right),w_k\left(t\right),{\mathbf{w}}_k\left(t\right)k=0,1,\dots ,N,$$

(36)

which represent the discrete approximation of the solution at the corresponding grid points. Note that $\rho$ and $\theta$ are approximated at the points indexed by j, that is, at the cell centers, while the velocity and the microrotation velocity components are approximated at the points indexed by k, which could be considered the cell edges. The spatial partial derivatives are approximated by using simple centered finite differences of the form

$${\displaystyle {\partial }_{y}g(y_l,t)\approx \delta g_l=\frac{g_{l+\frac{1}{2}}-g_{l-\frac{1}{2}}}{h},}$$

(37)

for $l=j$ or $l=k$. Note that if the function is approximated at the cell edge, its gradient is approximated at the cell center and vice versa. The defined arrangement of the variables on the staggered grid ensures that the velocity components, which physically depend on the pressure gradients, must be approximated at the cell edges, which is consistent with the position of the evaluation of the pressure gradient with (37), since the pressure is given by $p=R\rho \theta$, where $\rho$ and $\theta$ are approximated at the cell centers. Similarly, the density and temperature are influenced by the velocity gradients and are calculated at the same locations. This leads to a more accurate and stable coupling between pressure, velocity, and temperature fields and makes the finite difference approximations consistent with the physics being represented.

The functions (35)–(36) are determined as the solutions of the ordinary differential equation system obtained by introducing spatial discretization to the system (21)–(26), where the spatial derivatives are approximated by (37). The following system is obtained for the functions ${\rho }_j\left(t\right)$, $u_k\left(t\right)$, ${\mathbf{u}}_k\left(t\right),w_k\left(t\right),{\mathbf{w}}_k\left(t\right),{\theta }_j\left(t\right),$$j={\displaystyle \frac{1}{2}},\dots ,N-{\displaystyle \frac{1}{2}}$, and $k=1,\dots ,N-1$:

$${\dot{\rho }}_j=-L^{-1}{\rho }_j^2\delta u_j,$$

(38)

$${\dot{u}}_k=-L^{-1}R\delta {\left(\rho \theta \right)}_k+L^{-2}(\lambda +2\mu )\delta {\left(\rho \delta u\right)}_k,$$

(39)

$${\dot{\mathbf{u}}}_k=L^{-2}(\mu +{\mu }_r)\delta {\left(\rho \delta \mathbf{u}\right)}_k+2L^{-1}{\mu }_r\nabla \times {\mathbf{w}}_k,$$

(40)

$$j_I{\dot{w}}_k=L^{-2}(c_0+2c_d)\delta {\left(\rho \delta w\right)}_k-4{\mu }_r{\displaystyle \frac{w_k}{{\rho }_k}},$$

(41)

$$j_I{\dot{\mathbf{w}}}_k=L^{-2}(c_d+c_a)\delta {\left(\rho \delta \mathbf{w}\right)}_k-{\displaystyle \frac{4{\mu }_r}{{\rho }_k}}{\mathbf{w}}_k+2L^{-1}{\mu }_r\nabla \times {\mathbf{u}}_k,$$

(42)

$$\begin{array}{c}\hfill c_v{\dot{\theta }}_j=L^{-2}{k}_{\theta }\delta {\left(\rho \delta \theta \right)}_j-L^{-1}R{\rho }_j{\theta }_j\delta u_j+L^{-2}(\lambda +2\mu )\rho {\left({\partial }_{y}u\right)}^2\\ \hfill +L^{-2}(\mu +{\mu }_r)\rho {\left|{\partial }_y{\mathbf{u}}_j\right|}^2+4{\mu }_r{\displaystyle \frac{w_j^2}{{\rho }_j}}+4{\mu }_r{\displaystyle \frac{|{\mathbf{w}}_j{|}^2}{{\rho }_j}}-4L^{-1}{\mu }_r(\nabla \times {\mathbf{u}}_j)·{\mathbf{w}}_j+\\ \hfill L^{-2}(c_0+2c_d){\rho }_j{\left(\delta w_j\right)}^2+L^{-2}(c_d+c_a)\rho {\left|\delta {\mathbf{w}}_j\right|}^2,\end{array}$$

(43)

The values of the functions $u_k\left(t\right),{\mathbf{u}}_k\left(t\right),w_k\left(t\right)$, and ${\mathbf{w}}_k\left(t\right),$ as well as $\delta {\theta }_k\left(t\right)$, for $k=0,N$, that is, at the domain boundaries, are obtained from the boundary conditions. The initial conditions for the approximate functions are determined as the mean values of the initial functions over the corresponding cells.

In the described numerical method, the key point is the use of staggered grid discretization, where approximations for velocity and microrotation velocity are considered at the cell edges, and those for density and temperature are considered at the cell centers. Note that, in this way, the velocities and heat fluxes are located directly at the cell edges, leading to a natural and accurate application of the boundary conditions (27)–(29) in the numerical scheme. In this way, the stability of the solution is improved by reducing the potential for spurious modes and numerical artifacts from the boundaries. In [6,8], it was proved that the numerical methods defined in this way converge in 1D cases. Due to the similarity of the problem, we expect that convergence will also be achieved for the system under consideration.

Finally, the ODE system (38)–(43) obtained after the described spatial discretization can be solved numerically by using a standard numerical scheme. In this work, the explicit strong-stability-preserving Runge–Kutta method is used [17]. The steps of the finite difference scheme are summarized in Algorithm 1 below.

Algorithm 1 The finite difference method.

Require: Spatial domain $[0,h]$, $T>0$ defining temporal domain $[0,T]$; Physical parameters in (1)–(8); Initial conditions: ${\rho }_0$, ${\mathbf{v}}_0$, ${\mathit{\omega }}_0$, ${\theta }_0$ in (15)–(16); Boundary conditions: defined by (13)–(14): Number of discretization points N, time step for temporal discretization $\Delta t$ Ensure: Approximate solutions at time moments $t^m=m\Delta t\le T$, $m=1,2,\dots$: ${\rho }_j\left(t^m\right),{\theta }_j\left(t^m\right),u_k\left(t^m\right),{\mathbf{u}}_k\left(t^m\right),w_k\left(t^m\right),{\mathbf{w}}_k\left(t^m\right),$ for $j={\displaystyle \frac{1}{2}},\dots ,N-{\displaystyle \frac{1}{2}}$, $k=0,\dots ,N$. 1: Let initially: Determine the transformation map (17) at the time moment $t=0$, and transform the initial conditions using (30)–(32). Apply the initial conditions to determine the initial values of the approximate functions: ${\rho }_j\left(0\right)={\rho }_0\left(y_j\right),{\theta }_j\left(0\right)={\theta }_0\left(y_j\right),j={\displaystyle \frac{1}{2}},\dots ,N-{\displaystyle \frac{1}{2}}$ $u_k\left(0\right)=u_0\left(y_k\right),{\mathbf{u}}_k\left(0\right)=(0,v_{20}\left(y_k\right),v_{30}\left(y_k\right))$ $w_k\left(0\right)=w_0\left(y_k\right),{\mathbf{w}}_k\left(0\right)=(0,{\omega }_{20}\left(y_k\right),{\omega }_{30}\left(y_k\right)),k=1,\dots ,N-1$ Apply the boundary conditions at the initial time moment to determine: $u_k\left(0\right),{\mathbf{u}}_k\left(0\right),w_k\left(0\right),{\mathbf{w}}_k\left(0\right),\delta {\theta }_k\left(0\right),k=0,N$ 2: For $t^m=t^1,t^2,\dots$: Apply a step of the Runge-Kutta method (or an alternative algorithm) to the ODE system (38)–(43) to compute the approximate solution values at time $t^m$. These values are determined based on the approximations from previous time steps and correspond to the solution in the “inner” cells: $j={\displaystyle \frac{1}{2}},\dots ,N-{\displaystyle \frac{1}{2}}$ and $k=1,\dots ,N-1$, Apply the boundary conditions at time moment $t^m$ to determine: $u_k\left(t^m\right),{\mathbf{u}}_k\left(t^m\right),w_k\left(t^m\right),{\mathbf{w}}_k\left(t^m\right),\delta {\theta }_k\left(t^m\right),k=0,N.$

4. Faedo–Galerkin Approximations

The second approach for the numerical approximation of the proposed system is based on the Faedo–Galerkin method. This method is usually used to prove the local existence of the solution to a problem of this type by proving that the numerical solutions converge to the solution of the problem. In this work, the proposed scheme is used to analyze the behavior of the solution in the chosen experiments as well as for comparison with the numerical solutions obtained by using the finite difference method.

In order to apply the Faedo–Galerkin method, it is necessary to first homogenize the boundary conditions, similarly to [5]. For this purpose, we define the following functions:

$$U(y,t)=u(y,t)-g(y,t),\mathbf{U}(y,t)=\mathbf{u}(y,t)-\mathbf{h}(y,t),$$

(44)

where

$$\begin{array}{ccc}\hfill g(y,t)& =& {\displaystyle {\displaystyle \frac{u_R\left(t\right)-u_L\left(t\right)}{A\left(t\right)}}{\int }_0^y{\displaystyle \frac{ds}{\rho (s,t)}}+u_L\left(t\right),}\hfill \\ \hfill \mathbf{h}(y,t)& =& {\displaystyle {\displaystyle \frac{{\mathbf{u}}_R\left(t\right)-{\mathbf{u}}_L\left(t\right)}{A\left(t\right)}}{\int }_0^y{\displaystyle \frac{ds}{\rho (s,t)}}+{\mathbf{u}}_L\left(t\right),}\hfill \end{array}$$

(45)

and $A\left(t\right)$ is defined by (34). After introducing the new variables (44) into the boundary conditions, the new functions should satisfy the homogeneous boundary conditions

$$U(0,t)=U(1,t)=0,\mathbf{U}(0,t)=\mathbf{U}(1,t)=\mathbf{0},$$

and the initial conditions

$$U(y,0)=u_0\left(y\right)-{\displaystyle \frac{u_R\left(0\right)-u_L\left(0\right)}{A\left(0\right)}}{\int }_0^y{\displaystyle \frac{dy}{{\rho }_0\left(y\right)}}-u_L\left(0\right),$$

(46)

$$\mathbf{U}(y,0)={\mathbf{u}}_0\left(y\right)-{\displaystyle \frac{{\mathbf{u}}_R\left(t\right)-{\mathbf{u}}_L\left(t\right)}{A\left(0\right)}}{\int }_0^y{\displaystyle \frac{dy}{{\rho }_0\left(y\right)}}-{\mathbf{u}}_L\left(0\right),$$

(47)

while the system becomes

$${\partial }_t\rho +{\displaystyle \frac{1}{L}}{\rho }^2{\partial }_{y}U+{\displaystyle \frac{u_R\left(t\right)-u_L\left(t\right)}{LA\left(t\right)}}\rho =0,$$

(48)

$${\partial }_{t}U=-L^{-1}R{\partial }_y\left(\rho \theta \right)+L^{-2}(\lambda +2\mu ){\partial }_y\left(\rho {\partial }_{y}U\right)-{\partial }_{t}g,$$

(49)

$${\partial }_t\mathbf{U}=L^{-2}(\mu +{\mu }_r){\partial }_y\left(\rho {\partial }_y\mathbf{U}\right)+2L^{-1}{\mu }_r\nabla \times \mathbf{w}-{\partial }_t\mathbf{h},$$

(50)

$$j_I{\partial }_{t}w=L^{-2}(c_0+2c_d){\partial }_y\left(\rho {\partial }_{y}w\right)-4{\mu }_r{\displaystyle \frac{w}{\rho }},$$

(51)

$$j_I{\partial }_t\mathbf{w}=L^{-2}(c_d+c_a){\partial }_y\left(\rho {\partial }_y\mathbf{w}\right)-{\displaystyle \frac{4{\mu }_r}{\rho }}\mathbf{w}+2L^{-1}{\mu }_r\nabla \times \mathbf{U}+2L^{-1}{\mu }_r\nabla \times \mathbf{h},$$

(52)

$$\begin{array}{c}\hfill c_v{\partial }_t\theta =L^{-2}{k}_{\theta }{\partial }_y\left(\rho {\partial }_y\theta \right)-L^{-1}R\rho \theta ({\partial }_{y}U+{\partial }_{y}g)+L^{-2}(\lambda +2\mu )\rho {({\partial }_{y}U+{\partial }_{y}g)}^2\\ \hfill +L^{-2}(\mu +{\mu }_r)\rho |{\partial }_y{\mathbf{V}|}^2+2L^{-2}(\mu +{\mu }_r)\rho {\partial }_y\mathbf{V}{\partial }_y\mathbf{h}+L^{-2}(\mu +{\mu }_r)\rho {\left|{\partial }_y\mathbf{h}\right|}^2\\ \hfill +4{\mu }_r{\displaystyle \frac{w^2}{\rho }}+4{\mu }_r{\displaystyle \frac{{\left|\mathbf{w}\right|}^2}{\rho }}-4L^{-1}{\mu }_r(\nabla \times \mathbf{U})·\mathbf{w}\\ \hfill -4L^{-1}{\mu }_r(\nabla \times \mathbf{h})·\mathbf{w}+L^{-2}(c_0+2c_d)\rho {\left({\partial }_{y}w\right)}^2+L^{-2}(c_d+c_a)\rho {\left|{\partial }_y\mathbf{w}\right|}^2,\end{array}$$

(53)

We will find a numerical solution to the problem (48)–(53) by using the Faedo–Galerkin method, as described below. For a chosen $n\in \mathbf{N}$, we denote by

$$\left({\rho }^n,U^n,{\mathbf{U}}^n,w^n,{\mathbf{w}}^n,{\theta }^n\right),n\in \mathbf{N}$$

(54)

the Faedo–Galerkin approximations of the functions $\rho$, u, $\mathbf{U}$, w, $\mathbf{w}$, and $\theta$, where

$$U^n(y,t)=\sum _{i=1}^n{U}_i^n\left(t\right)sin\left(\pi iy\right),{\mathbf{U}}^n(y,t)=\sum _{i=1}^n{\mathbf{V}}_i^n\left(t\right)sin\left(\pi iy\right),$$

(55)

$$w^n(y,t)=\sum _{i=1}^n{w}_i^n\left(t\right)sin\left(\pi iy\right),{\mathbf{w}}^n(y,t)=\sum _{i=1}^n{\mathbf{w}}_i^n\left(t\right)sin\left(\pi iy\right),$$

(56)

$${\theta }^n(y,t)=\sum _{k=0}^n{\theta }_k^n\left(t\right)cos\left(\pi ky\right).$$

(57)

The coefficients of the approximate functions are obtained as a solution to the ordinary differential equation system obtained in accordance with the Faedo–Galerkin method applied to (50)–(53). The system reads as

$${\dot{U}}_i^n\left(t\right)=2{\int }_0^1\left[-L^{-1}R{\partial }_y\left({\rho }^n{\theta }^n\right)+L^{-2}(\lambda +2\mu ){\partial }_y\left({\rho }^n{\partial }_y{u}^n\right)-{\partial }_t{g}^n\right]sin\left(\pi iy\right)dy,$$

(58)

$${\dot{\mathbf{U}}}_i^n\left(t\right)=2{\int }_0^1\left[L^{-2}(\mu +{\mu }_r){\partial }_y\left({\rho }^n{\partial }_y{\mathbf{U}}^n\right)+2L^{-1}{\mu }_r\nabla \times {\mathbf{w}}^n-{\partial }_t{\mathbf{h}}^n\right]sin\left(\pi iy\right)dy,$$

(59)

$${\dot{w}}_i^n\left(t\right)=2{\int }_0^1\left[{\displaystyle \frac{c_0+2c_d}{L^2{j}_I}}{\partial }_y\left({\rho }^n{\partial }_y{w}^n\right)-4{\displaystyle \frac{{\mu }_r}{j_I}}{\displaystyle \frac{w^n}{{\rho }^n}}\right]sin\left(\pi iy\right)dy,$$

(60)

$${\dot{\mathbf{w}}}_i^n\left(t\right)=2{\int }_0^1\left[{\displaystyle \frac{c_d+c_a}{L^2{j}_I}}{\partial }_y\left({\rho }^n{\partial }_y{\mathbf{w}}^n\right)-{\displaystyle \frac{4{\mu }_r}{j_I}}{\displaystyle \frac{{\mathbf{w}}^n}{{\rho }^n}}+2{\displaystyle \frac{{\mu }_r}{L{j}_I}}\left(\nabla \times {\mathbf{U}}^n+\nabla \times {\mathbf{h}}^n\right)\right]sin\left(\pi iy\right)dy,$$

(61)

$$\begin{array}{c}\hfill {\dot{\theta }}_k^n\left(t\right)={\lambda }_k{\int }_0^1\left[{\displaystyle \frac{k}{L^2{c}_v}}{\partial }_y\left({\rho }^n{\partial }_y{\theta }^n\right)-{\displaystyle \frac{R}{L{c}_v}}{\rho }^n{\theta }^n{\partial }_y{U}^n-L^{-1}R\rho \theta {\partial }_y{g}^n\right.\\ \hfill +{\displaystyle \frac{\lambda +2\mu }{L^2{c}_v}}{\rho }^n{\left({\partial }_y{U}^n\right)}^2+2L^{-2}(\lambda +2\mu )\rho {\partial }_{y}U{\partial }_y{g}^n+L^{-2}(\lambda +2\mu )\rho {\left({\partial }_y{g}^n\right)}^2\\ \hfill +{\displaystyle \frac{\mu +{\mu }_r}{L^2{c}_v}}{\rho }^n|{\partial }_y{\mathbf{U}}^n{|}^2+2{\displaystyle \frac{\mu +{\mu }_r}{L^2{c}_v}}{\rho }^n{\partial }_y{\mathbf{U}}^n{\partial }_y{\mathbf{h}}^n+{\displaystyle \frac{\mu +{\mu }_r}{L^2{c}_v}}{\rho }^n{\left|{\partial }_y{\mathbf{h}}^n\right|}^2\\ \hfill +4{\displaystyle \frac{{\mu }_r}{c_v}}{\displaystyle \frac{{\left(w^n\right)}^2}{{\rho }^n}}+4{\displaystyle \frac{{\mu }_r}{c_v}}{\displaystyle \frac{|{\mathbf{w}}^n{|}^2}{{\rho }^n}}-4{\displaystyle \frac{{\mu }_r}{L{c}_v}}(\nabla \times {\mathbf{V}}^n)·{\mathbf{w}}^n-4{\displaystyle \frac{{\mu }_r}{L{c}_v}}(\nabla \times {\mathbf{h}}^n)·{\mathbf{w}}^n\\ \hfill \left.+{\displaystyle \frac{c_0+2c_d}{L^2{c}_v}}{\rho }^n{\left({\partial }_y{w}^n\right)}^2+{\displaystyle \frac{c_d+c_a}{L^2{c}_v}}{\rho }^n{\left|{\partial }_y{\mathbf{w}}^n\right|}^2\right]cos\left(\pi ky\right)dy,\end{array}$$

(62)

$${\dot{z}}_m^n\left(t\right)=u_m^n\left(t\right)\lambda \left(t\right),$$

(63)

for $i,m=1,\dots ,n$, and $k=0,1,\dots ,n$, where ${\displaystyle \lambda \left(t\right)=exp\left\{-\frac{1}{L}{\int }_0^t\frac{u_R\left(\tau \right)-u_L\left(\tau \right)}{A\left(\tau \right)}d\tau \right\}}$ and

$$\begin{array}{ccc}\hfill g^n(y,t)& =& {\displaystyle \frac{u_R\left(t\right)-u_L\left(t\right)}{A\left(t\right)}{\int }_0^y{\displaystyle \frac{ds}{{\rho }^n(s,t)}}+u_L\left(t\right),}\hfill \\ \hfill {\mathbf{h}}^n(y,t)& =& {\displaystyle \frac{{\mathbf{u}}_R\left(t\right)-{\mathbf{u}}_L\left(t\right)}{A\left(t\right)}{\int }_0^y{\displaystyle \frac{ds}{{\rho }^n(s,t)}}+{\mathbf{u}}_L\left(t\right).}\hfill \end{array}$$

(64)

We find the approximation ${\rho }^n$ of the function $\rho$ as the solution to the initial problem

$${\partial }_t{\rho }^n+{\displaystyle \frac{1}{L}}{\left({\rho }^n\right)}^2{\partial }_y{U}^n+{\displaystyle \frac{u_R\left(t\right)-u_L\left(t\right)}{LA\left(t\right)}}{\rho }^n=0,{\rho }^n(y,0)={\rho }_0\left(y\right),$$

(65)

which can be written in the form

$${\rho }^n(y,t)={\displaystyle \frac{{\rho }_0\left(y\right)\lambda \left(t\right)}{{\displaystyle 1+\frac{{\rho }_0\left(y\right)}{L}{\int }_0^t{\partial }_y{U}^n(y,\tau )\lambda \left(\tau \right)d\tau }}}={\displaystyle \frac{{\rho }_0\left(y\right)\lambda \left(t\right)}{{\displaystyle 1+\frac{{\rho }_0\left(y\right)}{L}\sum _{i=1}^{n}i\pi cos\left(\pi ix\right)z_i^n\left(t\right)}}}.$$

(66)

The initial conditions for Equations (58)–(62) are obtained as the Fourier coefficients of the initial functions and $z_m^n\left(0\right)=0.$ The functions on the right-hand side of the approximating system (58)–(63) satisfy the conditions of the Cauchy–Picard theorem; thus, the existence of its solution is ensured.

The numerical solution of the system (58)–(63) is determined with the same method as in the finite difference scheme, that is, with the strong-stability-preserving explicit Runge–Kutta method. The right-hand sides of Equations (58)–(62) are approximated numerically by using the Gauss–Legendre quadrature formulas of the order 10, similarly to [9]. In all expressions in which the spatial derivatives of the term ${\rho }^n$ occur, partial integration is performed first and then the integrals are approximated accordingly. It can be shown that the partial derivatives of the terms $g^n$ and ${\mathbf{h}}^n$, which are defined by (64), can be determined as

$${\partial }_y{g}^n(y,t)={\displaystyle \frac{u_R\left(t\right)-u_L\left(t\right)}{A\left(t\right)}}{\displaystyle \frac{1}{{\rho }^n(y,t)}},{\partial }_y{\mathbf{h}}^n(y,t)={\displaystyle \frac{{\mathbf{u}}_R\left(t\right)-{\mathbf{u}}_L\left(t\right)}{A\left(t\right)}}{\displaystyle \frac{1}{{\rho }^n(y,t)}}$$

(67)

and

$$\begin{array}{ccc}\hfill {\partial }_t{g}^n(y,t)& =& {\displaystyle \frac{{\dot{u}}_R\left(t\right)-{\dot{u}}_L\left(t\right)}{A\left(t\right)}}{\int }_0^y{\displaystyle \frac{dy}{{\rho }^n(y,t)}}+{\displaystyle \frac{{\dot{u}}_R\left(t\right)-{\dot{u}}_L\left(t\right)}{LA\left(t\right)}}{U}^n(y,t)+{\dot{u}}_L\left(t\right),\hfill \end{array}$$

(68)

$$\begin{array}{ccc}\hfill {\partial }_t{\mathbf{h}}^n(y,t)& =& {\displaystyle \frac{{\dot{\mathbf{u}}}_R\left(t\right)-{\dot{\mathbf{u}}}_L\left(t\right)}{A\left(t\right)}}{\int }_0^y{\displaystyle \frac{dy}{{\rho }^n(y,t)}}+{\displaystyle \frac{{\mathbf{u}}_R\left(t\right)-{\mathbf{u}}_L\left(t\right)}{LA\left(t\right)}}{U}^n(y,t)+{\dot{\mathbf{u}}}_L\left(t\right).\hfill \end{array}$$

(69)

Similarly to the finite difference method, we provide here Algorithm 2 summarizing the steps of the Faedo–Galerkin method.

Algorithm 2 The Faedo–Galerkin method.

Require: Spatial domain $[0,h]$, $T>0$ defining temporal domain $[0,T]$; Physical parameters in (1)–(8); Initial conditions: ${\rho }_0$, ${\mathbf{v}}_0$, ${\mathit{\omega }}_0$, ${\theta }_0$ in (15)–(16); Boundary conditions: defined by (13)–(14): Number of terms in the Faedo-Galerkin approximations n, time step for temporal discretization $\Delta t$ Ensure: Approximate solutions at time moments $t^m=m\Delta t\le T$, $m=1,2,\dots$: ${\rho }^n(x,t^m),{\theta }^n(x,t^m),u^n(x,t^m),{\mathbf{u}}^n(x,t^m),w^n(x,t^m),{\mathbf{w}}^n(x,t^m).$ 1: Let initially: Determine the transformation map (17) at the time moment $t=0$, and transform the initial conditions using (30)–(32). Determine the initial conditions for the new variables $U^n,{\mathbf{U}}^n$ using (46), (47). Determine the Fourier coefficients of the initial functions to obtain the initial values for the coefficients of FG approximate functions: $U_i^n\left(0\right),{\mathbf{U}}_i^n\left(0\right)$, $w_i^n\left(0\right),{\mathbf{w}}_i^n\left(0\right)$ and ${\theta }_k^n\left(0\right)$, $i=1,\dots ,n$, $k=0,1,\dots ,N$. 2: For $t^m=t^1,t^2,\dots$: Apply the partial integration to the integrals on RHS of system (58)–(62), take into account the boundary conditions if possible, and perform the numerical integration to approximate the remaining integrals. Apply a step of the Runge-Kutta method (or an alternative algorithm) to the ODE system (58)–(62) to compute the coefficients of the approximate functions at time $t^m$: $U_i^n\left(t^m\right),{\mathbf{U}}_i^n\left(t^m\right)$, $w_i^n\left(t^m\right),{\mathbf{w}}_i^n\left(t^m\right)$, ${\theta }_k^n\left(t^m\right)$, and $z_i^n\left(t^m\right)$, $i=1,\dots ,n$, $k=0,1,\dots ,n$ Use (66), (55)–(57) to obtain the Faedo-Galerkin approximations: ${\rho }^n(x,t^m),U^n(x,t^m),{\mathbf{U}}^n(x,t^m),w^n(x,t^m),{\mathbf{w}}^n(x,t^m),{\theta }^n(x,t^m)$. Apply (68) to determine $g^n(x,t^m)$ and ${\mathbf{h}}^n(x,t^m)$ and finally: $u^n(x,t^m)=U^n(x,t^m){+}^n(x,t^m)$ and ${\mathbf{u}}^n(x,t^m)={\mathbf{U}}^n(x,t^m)+{\mathbf{h}}^n(x,t^m)$.

5. Numerical Experiments

In this section, we compare both methods’ performance in the selected experiments. Since the methods are entirely different, the comparison of the obtained solutions will, to some extent, justify the validity of both methods.

In both considered problems, the physical parameters of the system were set to $L=1\left[{\displaystyle \frac{kg}{m^2}}\right]$; $c_d=c_0=c_v={10}^{-3}\left[{\displaystyle \frac{kg}{ms}}\right]$; $R=1\left[{\displaystyle \frac{J}{kgK}}\right]$; $j_I=1\left[m^2\right]$; $c_a=0\left[{\displaystyle \frac{kg}{ms}}\right]$; $\lambda =-2·{10}^{-3}\left[{\displaystyle \frac{kg}{ms}}\right]$; ${\mu }_r={10}^{-3}\left[{\displaystyle \frac{kg}{ms}}\right]$; $\mu =3·{10}^{-3}\left[{\displaystyle \frac{kg}{ms}}\right]$; and $k_{\theta }=0.024\left[K\right]$. Since the parameters used in the model, such as the microrotation coefficients, are still under investigation and their values for compressible fluids are currently unknown, we lacked sufficient information to simulate the real flow scenario. The parameters were chosen in such a way that the changes occurring within a relatively short time interval were clearly visible.

In the presented examples, the type of flow described in Section 2 was considered. We constructed approximate solutions to the problem (21)–(29), formulated in mass Lagrangian coordinates, and most of the results we present in the same coordinate system.

5.1. Example 1

We first considered the numerical solutions of the problem (21)–(29) defined with the following initial functions already given in mass Lagrangian coordinates:

$$\begin{array}{cc}\hfill {\rho }_0\left(y\right)& =1,\hfill \\ \hfill u_0\left(y\right)& =sin\left(\pi y\right),{\mathbf{u}}_0\left(y\right)=sin\left(\pi y\right)(0,1,1),\hfill \\ \hfill w_0\left(y\right)& =sin\left(2\pi y\right),{\mathbf{w}}_0\left(y\right)=sin\left(2\pi y\right)(0,1,1),\hfill \\ \hfill {\theta }_0\left(y\right)& =2+cos\left(\pi y\right),\hfill \end{array}$$

The boundary conditions were defined by (27)–(29) for three different cases, where the functions defining the boundary conditions modeling the motion of the plates were given with the following functions:

• Case A 

$$\begin{array}{cc}\hfill u_L\left(t\right)& =-0.5sin\left(\pi t\right),u_R\left(t\right)=0.5sin\left(\pi t\right),\hfill \\ \hfill {\mathbf{u}}_L\left(t\right)& =\mathbf{0},{\mathbf{u}}_R\left(t\right)=0.5sin\left(\pi t\right)(0,1,1);\hfill \end{array}$$

(70)

• Case B 

$$u_L\left(t\right)=0,u_R\left(t\right)=\left\{\begin{array}{cc}0.5t,\hfill & t\le 1\hfill \\ 0.5,\hfill & t>1\hfill \end{array}\right.,{\mathbf{u}}_L\left(t\right)=\mathbf{0},{\mathbf{u}}_R\left(t\right)=\mathbf{0};$$

(71)

• Case C 

$$u_L\left(t\right)=0,u_R\left(t\right)=0,{\mathbf{u}}_L\left(t\right)=\mathbf{0},{\mathbf{u}}_R\left(t\right)=\mathbf{0}.$$

(72)

In Case A, both plates moved; in Case B, the upper plate only moved in the x-direction while the lower plate remained fixed; and in Case C, both plates were fixed.

5.1.1. Case A

In the calculations, the finite difference numerical method with $N=64$ spatial cells and the Faedo–Galerkin method with $n=15$ expansion terms were used. The time step used for temporal discretization was adjusted so that the stability properties were satisfied and we set it here to $\Delta t={10}^{-3}$ for both methods.

The aim is to analyze the solutions in the presence of non-homogeneous boundary conditions and to compare the proposed computational approaches. For clearer visualization and comparison, the solutions in Figure 3, Figure 4, Figure 5 and Figure 6 are shown in mass Lagrangian coordinates, that is, they are displayed in a fixed domain. The numerical solutions obtained with both methods for the selected moments in time are compared. It can be observed that the solutions generally agree well. However, in the case of the Faedo–Galerkin method, minor oscillations occurred in the solution, which is particularly evident in Figure 6. Of course, the intensity of these oscillations decreased as more terms were included in the Faedo–Galerkin approximation.

To better illustrate the influence of non-homogeneous boundary conditions on the solutions, Figure 7 shows the temporal evolution of the solution using contour plots in Eulerian coordinates. Since the simulation was performed in mass Lagrangian coordinates, it was necessary to track the evolution of the domain points by applying Equation (20), that is, by solving the corresponding differential equation ${\displaystyle \frac{dx(y,t)}{dt}}=u(y,t)$. This visualization clearly shows both the change in the domain over time and the corresponding changes in the solution values.

5.1.2. Case B

Since the novelty of this study lies in the investigation of the flow of a micropolar fluid in a time-varying domain, we investigated the effect of boundary motion by considering the problem with non-homogeneous boundary conditions only in the x-direction. These conditions were defined by the functions given in (71). The parameters of the numerical methods were same as those in Case A.

The evolution of the problem domain can be observed in Figure 8, which shows the trajectories of the discretization points $y_k$ used in the finite difference method. The trajectories corresponding to $y_0=0$ and $y_N=1$ define, in fact, the boundaries of the domain at each moment in time. The solutions obtained for this case are shown in Figure 9 and Figure 10, where the numerical results for $\rho$ and u from both methods are compared for the selected moments in time. For clarity, the solutions are again shown in mass Lagrangian coordinates in a fixed domain. Since the boundary conditions were set such that for $t>1$ the velocity of the upper plate remained constant at $0.5$, while the lower plate remained fixed, it was physically justified to expect the solutions to converge to a steady state. It can be observed that the density decreased and converged to zero, which was expected as the mass was conserved, while the domain expanded with time. Moreover, since the velocity was 0 for the bottom plate and $0.5$ for the upper plate, the velocity u clearly converged to a linear profile within the domain, as can be seen from the solution at $t=100$.

Regarding the behavior of the numerical methods, it is visible in Figure 9 and Figure 10 that the Faedo–Galerkin method exhibited oscillations in regions with large solution gradients, while the finite difference method provided more stable results.

5.1.3. Case C

To analyze the convergence properties of the proposed numerical methods, we considered a problem with homogeneous boundary conditions to eliminate the influence of other boundary conditions on the behavior of the numerical solutions. The procedure is described below. Since the analytical solution was not available, we determined a reference solution using the Faedo–Galerkin method with a large number of terms, specifically $t=20$, in the expansion, for $t\in [0,1]$.

We examined first the convergence of the Faedo–Galerkin method with respect to the number of terms in the Fourier expansion. Specifically, we computed the $L_2$-norm of the difference between the reference solution and the approximate solutions obtained with varying numbers of terms in the expansion. The results, presented in Figure 11, demonstrate a clear convergence trend toward the reference solution. Although this is not a strict proof of the convergence rate, the trend is obvious.

A similar approach was employed to assess the convergence properties of the finite difference method. We used the same reference solution as before, namely the one obtained using the Faedo–Galerkin method with $n=20$. We then computed the $L_2$-norm of the difference between this reference solution and the solutions obtained using the finite difference method on grids with different numbers of discretization points, N. In this case, the deviations were measured only at the discretization points. The analysis was performed for N = 8, 16, 32, 64, 128, and 256 discretization points, where the time step was chosen to be sufficiently small to ensure numerical stability. The results are shown in Figure 12 on logarithmic scale. They clearly illustrate that the solutions converged as the number of the discretization points increased.

Furthermore, during the simulations, we observed that for large values of N in the finite difference method, the required time step had to be significantly reduced to ensure scheme stability. This, in turn, led to an accumulation of round-off errors, which ultimately affected the accuracy of the solutions in such cases.

In

Table 1. CPU times obtained with different computational parameters in Example 1.C.

Numerical Scheme	CPU Time (s)
finite difference, $\Delta t={10}^{-3}$, $N=32$	$0.44$
finite difference, $\Delta t={10}^{-3}$, $N=64$	$0.53$
finite difference, $\Delta t={10}^{-3}$, $N=128$	$0.95$
Faedo–Galerkin, $\Delta t={10}^{-3}$, $n=5$	$5.53$
Faedo–Galerkin, $\Delta t={10}^{-3}$, $n=10$	$10.46$
Faedo–Galerkin, $\Delta t={10}^{-3}$, $n=20$	$21.23$

, the computation times for both methods obtained with different parameters, for the final time $t_{max}=$ 10 s, using the time step $\Delta t={10}^{-3}$ are presented. By observing the numerical results, we concluded that in most of the considered cases, comparable accuracy was achieved with $N=64$ discretization points in the finite difference method and $n=20$ terms in the Faedo–Galerkin method. A comparison of the CPU times for these parameter settings revealed that the Faedo–Galerkin method was approximately 40 times more computationally expensive than the finite difference method. Overall, the finite difference method proved to be significantly more efficient than the Faedo–Galerkin method.

5.2. Example 2

In the second test example, it was assumed that the initial density belonged to $H^1\left((0,1)\right)$, while the other initial functions were smoothly differentiable:

$$\begin{array}{cc}\hfill {\rho }_0\left(y\right)& =\left|y^2-{\displaystyle \frac{1}{4}}\right|+1,\hfill \\ \hfill u_0\left(y\right)& =0,{\mathbf{u}}_0\left(y\right)=\mathbf{0},\hfill \\ \hfill w_0\left(y\right)& =4y^2(1-y^2),{\mathbf{w}}_0\left(y\right)=4y^2(1-y^2)(0,1,1),\hfill \\ \hfill {\theta }_0\left(y\right)& =0.1.\hfill \end{array}$$

The non-homogeneous boundary conditions (28) modeling moving plates were defined by (70), as in Example 1.

In this example, the chosen initial function for the density had lower smoothness and the Fourier expansions of the initial functions contained an infinite number of terms, in contrast to the Fourier expansion of the initial conditions in Example 1, which only had one term in the expansion. The aim of this example was to investigate the influence of the number of expansion terms on the computational results and to compare the methods in the case of less smooth initial conditions. In the finite difference scheme, a spatial discretization with $N=64$ cells was used, while in the Faedo–Galerkin method, we employed $n=20$ terms in the approximate functions. As before, the time step was set to $\Delta t={10}^{-3}$.

First, as in Case A in Example 1, we compared both methods at different moments in time, as shown in Figure 13, Figure 14, Figure 15 and Figure 16. The numerical results of both computational approaches show reasonable agreement. However, discrepancies occur near the points where the solution lacks smoothness. For example, deviations in the density can be observed near $y=0.5$ and in regions with steep gradients, such as at the right boundary of the velocity component $v_2$.

We also investigated how the number of Fourier series terms in the Faedo–Galerkin method affected the solution near such problematic points. The results at $t=7.5$ s are shown in Figure 17 and Figure 18. As observed in the previous figures, the presence of steep gradients or non-smooth solutions leads to oscillations in the Faedo–Galerkin method. This effect is particularly noticeable when fewer expansion terms were used. However, as the number of terms increased, the numerical solutions obtained with the Faedo–Galerkin method progressively approached those obtained by the finite difference method.

By increasing the density of the numerical grid in the finite difference method, or the number of terms, n, in the Faedo–Galerkin method, we expected that the approximate solutions would converge to the exact solutions of the observed system.

6. Conclusions

This study investigated compressible micropolar fluid flow between two moving parallel plates. The mathematical model was derived from a general three-dimensional micropolar fluid model, simplified under the assumption that the flow depends on only one spatial variable. The work contributes to the study of micropolar fluid dynamics in two main directions. The first contribution is the development of a formulation for the micropolar flow between two plates, whose motion is described by non-homogeneous boundary conditions for the velocity, leading to a time-dependent domain. By using mass Lagrangian coordinates, we transformed this variable-domain problem into a fixed-domain problem, which simplifies the numerical treatment.

Second, we introduced and analyzed robust computational strategies to solve this problem, in particular the finite difference method and the Faedo–Galerkin method. Great attention was paid to the discretization process to ensure stability and convergence. In the finite difference method, staggered grids were implemented to maintain the physical relationships between the approximate variables, preventing, in this way, possible instabilities. The Faedo–Galerkin method, on the other hand, was adapted by introducing a suitable substitution to reformulate the problem into one with homogeneous boundary conditions, making it compatible with the standard approach. Both methods were tested on selected numerical experiments, demonstrating good agreement in the results, which confirmed the validity of the proposed computational approaches. The finite difference method proved to be significantly more computationally efficient, while the Faedo–Galerkin method showed convergence as the number of expansion terms increased, albeit with higher computational costs. Furthermore, we analyzed the effects of smooth and non-smooth initial conditions on the numerical solutions, showing the tendency of the Faedo–Galerkin method to generate oscillations in regions where the solution was not smooth and in regions with steep gradients.

Several extensions of this work are possible. First, a deeper theoretical analysis of convergence and stability would strengthen the mathematical foundation of the numerical methods. Second, the numerical framework could be extended to account for additional physical effects, such as different gas state equations and more complex boundary conditions. Other possible directions include the exploration of adaptive mesh refinement techniques for the finite difference method, which could help to reduce the computational costs while maintaining accuracy, the use of implicit numerical schemes for temporal discretization in the finite difference case, and the application of finite element techniques. This study provides a basis for further research and a deeper understanding of the behavior of micropolar fluids in dynamic environments.