Temporal stability of multiple similarity solutions for porous channel flows with expanding or contracting walls

In this paper, the temporal stability of multiple similarity solutions (flow patterns) for the incompressible laminar fluid flow along a uniformly porous channel with expanding or contracting walls is analyzed. This work extends the recent results of similarity perturbations of [1] by examining the temporal stability with perturbations of general form (including similarity and non-similarity forms). Based on the linear stability theory, two-dimensional eigenvalue problems associated with the flow equations are formulated and numerically solved by a finite difference method on staggered grids. The linear stability analysis reveals that the stability of the solutions is same with that under perturbations of a similarity form within the range of wall expansion ratio α (−5 ≤ α ≤ 3 as in [1]). Further, it is found that the expansion ratio α has a great influence on the stability of type I flows: in the case of wall contraction (α < 0), the stability region of the cross-flow Reynolds number (R) increases as the contraction ratio (|α|) increases; in the case of wall expansion and 0 < α ≤ 1, the stability region increases as the expansion ratio (α) increases; in the case of 1 ≤ α ≤ 3, type Corresponding author Email address: p.lin@dundee.ac.uk (Ping Lin ) Preprint submitted to PHYSICS OF FLUIDS July 24, 2021 T hi s is th e au th or ’s p ee r re vi ew ed , a cc ep te d m an us cr ip t. H ow ev er , t he o nl in e ve rs io n of r ec or d w ill b e di ffe re nt fr om th is v er si on o nc e it ha s be en c op ye di te d an d ty pe se t. P L E A S E C IT E T H IS A R T IC L E A S D O I: 1 0 .1 0 6 3 /5 .0 0 5 1 8 4 6 I flows are stable for all R where they exist. The flows of other types (types II and III with −5 ≤ α ≤ 3 and type IV with α = 3) are always unstable. As a nonlinear stability analysis or a validation of the linear stability analysis, the original nonlinear two-dimensional time dependent problem with an initial perturbation of general form over those flow patterns is solved directly. It is found that the stability with the non-linear analysis is consistent to the linear stability analysis.


Introduction
The laminar flow in a porous channel with expanding or contracting walls has attracted much attention due to its wide applications in engineering and biomedicine, including transpiration cooling, phase sublimation, propellant burning, filtration, and blood transport in organisms. For example, the sublimation process of carbon dioxide, during which the walls expanded ( [2]); propellant burning in a rocket motor with regressing walls ( [3]); and fluid transport produced by expansion and contraction of a blood vessel ( [4]).
The earliest investigations of steady flows across permeable and stationary walls can be traced back to Berman [5]. In his study, the laminar, twodimensional flow of a viscous incompressible fluid in a porous channel with uniform injection (or suction) was considered. By assuming that the transverse velocity component was independent of the streamwise coordinate, the Navier-Stokes equations were reduced to a nonlinear ordinary differential equation with appropriate boundary conditions. Then Berman obtained an asymptotic expression for a small Reynolds number R by a perturbation method. A number of studies of porous channel flow followed. For example, Terrill [6] extended Berman's small R case and obtained series solutions for large R (for large suction), and Proudman [7] investigated the case of large R using an integral approach. Using the method of averages, Morduchow [8] obtained an analytical 2 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0051846
solution for the entire injection range. Yuan [9] provided a perturbation solution for high injection case, and later, Terrill [10] modified the work of Yuan and provided a more accurate solution.
The earliest studies for moving walls can be traced back to Brady and Acrivos [11]. In their study, an exact solution to the Navier-Stokes equations for the flow in a channel with an accelerating surface velocity was presented. Along similar lines, Dauenhauer and Majdalani [2] obtained a self-similar solution for a porous channel flow with expanding or contracting walls. They assumed that the wall expansion ratio α was a constant and reduced the Navier-Stokes equations to a boundary value problem of a fourth-order nonlinear ordinary differential equation that could be solved by a shooting method. In a later study, asymptotic solutions for this problem were presented by Majdalani et al. [4] for small R and by Majdalani and Zhou [12] for moderate-to-large R.
Zhou and Majdalani [3] also provided an analytical solution for slab rocket motors with regressing walls. Recently, Xu et al. [13] investigated multiple solutions of the case for which the wall expansion ratio α may be varied from α 0 to α 1 through some given functions, and concluded that the solutions quickly reached the steady state. More recently, Majdalani and Xuan [14] improved the results in [12] and obtained a complete asymptotic solution for the problem of channel flow with moving walls. In their work, a viscous boundary layer correction was provided to overcome the singular pressure distribution and its normal gradients near the midsection plane of the expanding porous channel.
Later, a wavelet-homotopy method was developed by Chen and Xu [15] to give solutions to this problem. For a porous tube with an expanding or contracting sidewall, analytical solutions for both large and small Reynolds number with small-to-moderate α were obtained by Saad and Majdalani [16] recently.
As for the stability of the solutions, Durlofsky and Brady [17] investigated the spatial stability of the solutions for two-dimensional porous wall channel and accelerating-wall channel flows under linear symmetric perturbations. For the same porous wall problem, Ferro and Gnavi [18] extended the results of Durlofsky and Brady to symmetric and asymmetric solutions, and analyzed the 3 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0051846
spatial stability of small perturbations of arbitrary shape. The temporal stability of these flows was examined by Zaturska et al. [19]. They proved that most of these flows were temporally unstable to two-dimensional antisymmetric perturbations. Later, Taylor et al. [20] generalized the work of Zaturska to three-dimensional flows. Watson et al. [21] investigated the temporal stability of asymmetric flows arising from a channel with porous and accelerating walls. For porous channel flows with expanding or contracting walls, the temporal stability analysis was presented in [1]. It is noted that all the perturbations used in the above temporal stability analysis are constrained to the form of the similarity transformation. We would also like to mention a few other recent work on linear stability analysis of relevant channel or duct flow problems. For flows in solid rocket motors, the stability was investigated in [22,23,24]. The temporal stability analysis of pressure-driven flows in channels patterned with superhydrophobic surfaces containing periodic grooves and ribs aligned longitudinally to the flow direction was performed by Yu et al. [25], and the stability of a pressure driven flow in a duct heated from below and subjected to a vertical magnetic field was studied by Qi et al. [26].
From both physical and mathematical points of view, the perturbation of a flow solution (no matter whether it is a similarity solution or not) is not necessarily in the similarity form. So to study the stability properly and accurately we have to consider the perturbation in a general form (including similarity and non-similarity forms). This is the purpose of this paper, that is, to investigate the temporal stability of similarity solutions (for flows in a channel with expanding or contracting walls) under perturbations of general form. The basic equations of the problem and the multiple solutions are described in Section 2.
The linear stability analysis of these solutions by numerical means is carried out in Section 3. The linear stability theory is based on linear approximation of the nonlinear equations, which does not cover the nonlinear temporal development of an initial perturbation. So in Section 4 a non-linear analysis is conducted by directly solving the nonlinear Navier-Stokes equations with small-amplitude initial perturbations of general form. Section 5 is devoted to the conclusions.

4
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0051846 2. Mathematical formulation of the flow problem The plotted streamlines correspond to a symmetric steady flow pattern.
Consider the two-dimensional, laminar and incompressible flow in a rectangular channel with two permeable and moving walls. As shown in Fig. 1, which depicts the cross section of the simulated domain. The channel height is 2d and the channel length is semi-infinite. Both sidewalls have the same permeability and expand or contract uniformly at a time-dependent rateḋ, where˙means the derivation oft. Additionally, withx representing the streamwise direction andȳ the normal direction, the corresponding streamwise and normal velocity components are defined asū andv, respectively. The over-bar is used to denote dimensional variables. Under these assumptions, let the velocity vector v = (ū,v), the general continuity and motion equations are given as where ∇ is the gradient operator and △ is the Laplace operator,p, ρ,t and ν are the dimensional pressure, density, time and kinematic viscosity, respectively.
The boundary conditions arē

5
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0051846
where v w is the injection velocity at the wall, which is assumed to be independent of position. A = v w /ḋ is a constant which is a measure of the wall permeability.
The condition (5) can be achieved by making the flow symmetrical with respect to the planex = 0, wherev is left free.
Next, we introduce the following scalings: then the following dimensionless equations are obtained.
The original boundary conditions become Here, v = (u, v) and is the wall expansion ratio. α > 0 implies the expansion and α < 0 the contraction. R is the cross-flow Reynolds number defined by R = dv w /ν. We can infer that R > 0 is for the injection and R < 0 for the suction. In the current study, we just consider the case for which R is time invariant. It follows that α is constant and can be specified by its initial valueḋ 0 d 0 /ν, where d 0 andḋ 0 are the initial channel half-height and expansion rate, respectively. Integrating (12), the channel height of the present solution will vary in time according to d = √ d 2 0 + 2ναt. The dimensionless flow problem admits an exact similarity solution of the

6
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0051846
For flows symmetric with respect to the midsection plane (ȳ = 0), the velocity satisfies the boundary conditions (3) and which are the same as (3) and (4) in [1], respectively. Further, the dimensionless boundary conditions become (9) and By using (13) into (7) and (8) , we obtain a differential equation for the similarity function f , the boundary conditions are In particular, there are symmetric steady solutions with where and Here, a prime denotes differentiation with respect to y. Particularly, equation (19) is Berman's classic equation in [5] when α = 0.
Remark 1. Some researchers also consider the asymmetric solutions satisfying the boundary conditions In this paper we shall mainly consider stability of symmetric steady solutions which satisfy the boundary conditions (20).

7
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0051846
The numerical solutions of (19) and (20)   I is in R 1 3 < R < ∞ and type IV is in R 1 3 < R < −0.796, types II and III exist in a common semi-infinite domain −∞ < R < R 2 3 . Where R 1 3 = −4.25 and R 2 3 = −9.545 are the common points for types I, IV and types II, III, respectively. These results are presented in Table 1. Note that ∞ and −∞ in this paper stand for relatively large or negatively large values, respectively.

Remark 2.
We can also present the bifurcation diagram as a function of α.
For example, Fig. 4 shows bifurcation diagrams for a couple of R values. For R = −14 (see Fig. 4(a)), there is just one symmetric solution (type I solution) for −5 ≤ α < −3.225, and there are three symmetric solutions (types I, II and III solutions) for −3.225 < α ≤ 2. The common point of types II and III is α = −3.225. For R = −30 (see Fig. 4(b)), there are three types (types I, II and III) solutions for −5 ≤ α ≤ 2. In this paper we consider the stability of steady solutions for a range of values of α, so we will not explore this type of bifurcation diagrams further.  This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.  This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0051846 α number of solutions found type designation existence ranges This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.  As described in [1], the characteristics of the four types of flows represented by the solutions are as follows:

PLEASE CITE THIS ARTICLE AS
(1) Type I covers the flows whose axial velocity profiles have a maximum at the center of the channel.
(2) Type II includes the flows whose axial velocity profiles have an inflection point and a maximum between the center of the channel and the wall and whose centerline velocity is positive for negative R far away from 0.
(3) Type III contains axial velocity profiles with the same form as type II but with reverse flow at the center of the channel.
(4) Type IV includes the flows which have reverse flow near the wall of the channel, and the wall velocity gradient (F ′′ (1)) for these flows increases rapidly with the increase of R.
The axial velocity profiles F ′ (y) for type I solutions with some injection and suction cross-flow Reynolds numbers over a range of wall expansion ratios are described in Figs. 5 and 6. In each case of R, increasing wall expansion ratio increases the axial velocity near the center of the channel and decreases near the wall. This behaviour is reversed for the case of contracting channel. In addition, the profiles for each case of α have a maximum at the center of the channel.

12
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0051846
For R > 0 (injection), they monotonically decrease to 0 at the wall, and the velocity at the centerline (F ′ (0)) is approximately equal to 1.57 as R → ∞. The axial velocity profiles for type II solutions are described in Fig. 7. For each case of α, it can be observed in Fig. 7(a) that the profiles have a minimum at the centerline and then pass through a maximum before going to zero at the wall. Further, for negative R far away from 0, the velocity for α ≥ 0 is close to 1 everywhere except in a boundary layer (see Fig. 7(b)), which is similar to that described for type I solutions with α ≤ 0. The axial velocity profiles for type III solutions are depicted in Fig. 8. For each case of α, these profiles have the same shape as those of type II, except that there is a region of reverse flow near the center of the channel at any R where these solutions exist. For types I, II and III solutions, the profiles (shown in Fig. 9(a)) have similar characteristics with those for types I, II and III solutions at α = 2, respectively.
For type IV solutions, the profiles (shown in Fig. 9(b)) are characterized by a rapid increase in the centerline velocity and the wall velocity gradient (F ′′ (1)) as R increases, and the development of reverse flow near the wall of the channel.

13
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.   14 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0051846

15
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0051846  The temporal stability analysis of above steady flows (denoted as U) under perturbations of the similarity form (13) is given in [1]. Although the similarity solutions are considered, from both physical and mathematical points of view, the perturbations are not necessarily of the similarity form. So it is not complete to examine the stability of the flows only for the perturbations of the similarity form. In this paper we investigate the temporal stability with perturbations of general form (including similarity and non-similarity forms). We shall adopt the numerical means later, and for numerical study, we can only deal with finite domain. We truncate the infinite domain to an artificial boundary at x = x r , and develop and impose a proper boundary condition at x = x r in order for the resulted steady solutions to be consistent with the similarity solutions and facilitate comparison with previous analysis in [1]. The conditions read:

16
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0051846
It is not difficult to verify that all steady state similarity solutions satisfy the proposed condition (25) at the artificial boundary.

Temporal stability analysis
Here we examine the linear temporal stability of above steady flows under perturbations of general form (including similarity and non-similarity forms), in order to determine whether such perturbations could destabilize a flow which is stable under perturbations of the similarity form (13). We write the perturbed velocity and pressure fields where P is the unperturbed pressure, v 1 and p 1 are infinitesimal perturbations for the steady flow U and P , respectively. Substituting (26) into the dimensionless equations (7), (8) and boundary conditions (22)- (25), and linearizing (8) for v 1 , we obtain the following linearized perturbation equations and boundary conditions The perturbations (v 1 ) are of general form and include those of the similarity form (13) considered in [1].
Based on the method of separation of variables, the perturbations v 1 and 17 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0051846 p 1 can be expressed in the following forms In (33),û(x, y),v(x, y) andp(x, y) are the amplitudes of the corresponding perturbations, s is the complex eigenvalue. The real part of s (Re(s)) represents the growth or decay rate of the perturbation. When α = 0, Re(s) represents the growth rate for R > 0, while for R < 0, the sign of t becomes negative and Re(s) represents the decay rate. When α < 0 (for contraction), t 1 is positive and Re(s) is the decay rate. When α > 0 (for expansion), we note that t is finite, and when t → (R/2α), the channel height d has already reached infinity. Hence, t 2 is also positive and Re(s) is the decay rate. That is, for α = 0 and R > 0, eigenvalues with positive real parts (Re(s)) indicate growing perturbations, so the instability is implied if there is an eigenvalue such that Re(s) > 0; while for the case of α ̸ = 0 and the case of α = 0 and R < 0, eigenvalues with negative real parts indicate growing perturbations, so the instability is implied if there is an eigenvalue such that Re(s) < 0. Especially when α > 0 (for expansion), the instability occurs at t → (R/2α). The imaginary part of s (Im(s)) represents the dimensionless frequency of the corresponding perturbation. If s is real, the perturbations either grow or decay monotonically.

18
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0051846
Substituting (33) into (27) and (28), we have the following eigenvalue problems: where G = −1 for α < 0 and G = +1 for α > 0, all associated with the boundary To overcome the difficulty of lacking a boundary condition for the pressure, the discretization of the eigenvalue problem (37) (for α > 0) associated with (38) is done on the staggered grid ( Fig. 20) introduced by Harlow and Welch [27] .
The corresponding finite difference scheme is given in Appendix A. Similar finite difference schemes have been constructed for eigenvalue problems associated with other α.
The eigenvalue pencil of the real unsymmetric eigenvalue problem (A3) satisfying (A4) contains real values and complex conjugate pairs. To detect the flow instability, we need to seek the eigenvalues with the maximal or minimal 19 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0051846
real part that corresponds to the least stable eigenvalues. For the case of α ̸ = 0 and the case of α = 0 and R < 0, the least stable eigenvalues are those with the minimal real part (i.e., the minimum decay rate). For α = 0 and R > 0, the opposite is true, and the least stable eigenvalues are those with the maximal real part (i.e., the maximum growth rate).
We also need to determine a proper artificial boundary (the truncated channel length) x = x r for the eigenvalue computation. It should not be too large in order to save the computational time. In the meanwhile it should not affect the stability study. Fig. 10 shows the minimal real part of the eigenvalues (marked by q) versus R for α = 1/2 with x r = 5, 10 and 20 and with a 10 × 800 mesh.
The figure suggests that the stability of the three types (types I, II and III) solutions behaves the same for the three choices of artificial channel lengths. We thus choose a smaller x r = 5 in all the following computations so as to reduce the overall computational cost. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0051846
The real part of the least stable eigenvalues for (36) and (37) is plotted versus R in Figs. 11-13. We mark the minimal real part of the eigenvalues as q in the case of α ̸ = 0. When α = 0, q represents the maximal real part for R > 0 and the minimum real part for R < 0.

21
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.   This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0051846
part of each of the eigenvalues is positive for α ̸ = 0, and negative for α = 0, namely no amplification of perturbations occurs and the injection flows are always stable. In addition, we note that except for a small range of R > 0, the least stable eigenvalues are not the same as that in the case of perturbations in similarity form (13), while the stability range is the same. As an example, for type I solutions with α = 2 and 10 ≤ R ≤ 50, the comparison of eigenvalues with the minimum decay rate is shown in Table 2. For each case of R, the decay rate of the present least stable perturbation is smaller than that of the least stable perturbation of [1], implying that the present perturbation decays more slowly. When R < 0 (that is, when there is suction), the stability of type I solutions varies for different α.  Table 3, which are consistent with those in [1]. For α < 0 (in the case of wall contraction), we find that the critical cross-flow Reynolds number decreases and the stability region increases as the contraction ratio (|α|) increases; for α > 0 (in the case of wall expansion), the critical R decreases and the stability region increases as the expansion ratio (α) increases. One possible explanation for this behaviour is as follows: in the case of wall contraction (α < 0), for larger contraction ratio (|α|), the channel half- This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0051846 decay more rapidly. We further note that, when α is larger, i.e., α ≥ 1, and Consequently, types II and III solutions with −5 ≤ α ≤ 2 are always unstable.
Further, we note that q → −1 as R → −∞ for type II solutions with α = 0.
When α = 3, the minimum decay rate (q) for a range of R is shown in Fig. 13.
We note that the value of q for type I solutions is positive for R 1 3 < R < ∞, the least stable perturbations therefore decay with time and type I solutions (with α = 3) are stable in this region. As opposed to type I solutions, the solutions of types II and III for −∞ < R < R 2 3 are unstable, since the perturbations are expected to grow in time as indicated by the negative value of q. Similarly, for type IV solutions, the value of q is negative for R 1 3 < R < −0.796 and hence type IV solutions are also unstable. 24 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. Remark 3. The vanishing of q(I) suggests the existence of bifurcation. For example, in the case of α = 1/2, the zero real eigenvalue at R = R ′ 1/2 = −5.905 corresponds to a "pitch-fork" bifurcation. Two types of asymmetric steady solutions of (19) subject to (21) appear at R ′ 1/2 and form two branches of the "pitch-fork" bifurcation. We name the types of these solutions I 1/2 and I ′ It is noted that we consider here not only perturbations of the similarity form (13), but also perturbations of the non-similarity (general) form. Although the linear stability (or instability) of the symmetric flows (with −5 ≤ α ≤ 3) obtained here is the same as that under the perturbations of the similarity form shown in [1], the least stable eigenvalues or the most unstable eigenvalues are not all the same. For the case of α ̸ = 0 and the case of α = 0 and R < 0, the minimal real part of the eigenvalues (i.e., the minimum decay rate) of the perturbations for some flows is smaller than that of [1]; for the case of α = 0 and R > 0, the maximal real part of the eigenvalues (i.e., the maximum growth rate) of the perturbations for some flows is larger than that of [1]. As a result, the least stable perturbations are expected to decay more slowly or grow faster.

PLEASE CITE THIS ARTICLE AS
A few more examples are given below.
For type I solutions with α = −0.5 and −1 ≤ R ≤ −0.2, the comparison of the least stable eigenvalues (i.e., the eigenvalues with the minimum decay rate) with those of [1] is shown in Table 4. In addition, the real parts of the streamwise velocity eigenfunctionû(x, −0.00125) and the normal velocity eigenfunctionv(x, 0) corresponding to the present results of the least stable eigenvalues (in Table 4) are illustrated in Fig. 15. The eigenvectorsv are normalised by using the corresponding 2-norm, so that the 2-norm of the eigenvectors is 1. It can be seen from Table 4 that the least stable eigenvalues for R = −1 and R = −0.8 are in 26 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0051846
good agreement with those of [1]. For each case of R = −1 and R = −0.8, we note that the perturbation (corresponding to the least stable eigenvalue) is of the similarity form (13), andû(x, y) of streamwise velocity perturbation andv(x, y) of normal velocity perturbation are real functions. Therefore,û(x, −0.00125) plotted in Fig. 15(a) is proportional to x, andv(x, 0) plotted in Fig. 15(b) is independent of x. Nevertheless, in each case of −0.6 ≤ R ≤ −0.2, the minimum decay rate is smaller than that of [1] (see Table 4), implying that the corresponding least stable perturbation decays more slowly. Further, we note that the perturbation is not of the similarity form, since Re(û(x, −0.00125)) plotted in Fig. 15(a) is not linear with respect to x, and Re(v(x, 0)) plotted in Fig.   15(b) changes with x.
For type III solutions with α = 2 and −55 ≤ R ≤ −35, Table 5 shows the comparison of the most unstable eigenvalues (i.e., the eigenvalues with the minimum decay rate) with those of [1]. In addition, the real parts of the streamwise velocity eigenfunctionû(x, −0.00125) and the normal velocity eigenfunction v(x, 0) for the present results of the most unstable eigenvalues (in Table 5) are illustrated in Fig. 16. The eigenvectorsv are also normalized by using the corresponding 2-norm, and the 2-norm of the eigenvectors is 1. The most unstable eigenvalue for R = −35 is in good agreement with that of [1] (see Table 5).
Moreover, the perturbation corresponding to this eigenvalue is of the similarity form. This is reflected in the results of the real parts of the eigenfunction componentsû(x, −0.00125) (plotted in Fig. 16(a)) andv(x, 0) (plotted in Fig.   16(b)) for R = −35. However, for each case of −55 ≤ R ≤ −40, the real part of the most unstable eigenvalue is smaller than that of [1] (see Table 5), indicating that the corresponding most unstable perturbation grows faster. Further, the perturbation is not of the similarity form as indicated by the results of the real parts ofû(x, −0.00125) (plotted in Fig. 16(a)) andv(x, 0) (plotted in Fig.   16(b)) for −55 ≤ R ≤ −40.

27
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.   Table 4. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.  Table 5.
When α ̸ = 0 and R < 0, a new variable t * = R 2α ln(1 − 2αt R ) is introduced so as to have a usual time interval 0 ≤ t * < ∞. Substituting t * into (7) and (8)   29 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
We solve (41) This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0051846
appear to be stable. Thus, we can obtain that a critical R exists between −6.4 and −6.3 for α = −1/2. Similarly, when α = 1/2, it can be seen from Fig. 17(e) that the axial velocity at R = −5.9 does not change significantly with time.
While the axial velocity at R = −6 shown in Fig. 17(f)

18(d) and 19(d))
, it can be seen that after a period of time, the flows of types II and III turn into the symmetric steady flow of type I.

32
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

Conclusion
In this numerical study, the multiple symmetric similarity solutions of a flow problem occurring in a uniformly porous channel with expanding (or contracting) walls are considered in a range of the wall expansion ratios α, say, [−5, 3].
We examine the linear temporal stability of these solutions under perturbations of general form (including similarity and non-similarity forms). Through a finite difference method on a staggered grid we solve two-dimensional eigenvalue problems associated with the linear stability analysis and the stability of these solutions is then obtained. That is, type I solutions in each case of −5 ≤ α ≤ 1/2 are only stable for a range of R (cross-flow Reynolds number), and type I solutions with 1 ≤ α ≤ 3 are stable for all R where they exist. Further, it is found that for α < 0 (in the case of wall contraction), the stable region of R increases as the contraction ratio (|α|) increases; for 0 < α ≤ 1 (in the case of wall expansion), the stable region increases as the expansion ratio (α) increases. So the expansion ratio α has a great influence on the stability of the flows of type I, and it seems that the presence of an inflection point of axial velocity (or the flow acceleration changes from decrease to increase) near the wall may stabilize the flow. In addition, other types of flows whose axial velocity profiles have an inflection point near the center of the channel are always unstable, suggesting 33 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0051846
that these flows may transition to turbulence prior to physically attaining these shapes. In other words, these flows may not be physically observable.
Although the stability (or instability) of these steady flows obtained here under perturbations of general form is the same as that under the perturbations of the similarity form shown in [1], the minimum decay rate or maximum growth rate of the perturbations are not all the same. For the case of α ̸ = 0 and the case of α = 0 and R < 0, the minimal real part of the eigenvalues (i.e., the minimum decay rate) of the perturbations for some flows is smaller than that of [1]; for the case of α = 0 and R > 0, the maximal real part of the eigenvalues (i.e., the maximum growth rate) of the perturbations for some flows is larger than that of [1]. As a result, the least stable perturbations are expected to decay more slowly or grow faster.
On the other hand, non-linear analysis has been carried out by directly solving the original nonlinear time dependent problem with an initial perturbation of general form. It is found that the stability results agree well with those obtained from the linear stability analysis.

Acknowledgements
This work is partially supported by the National Natural Science Foundation This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0051846 [27]. The finite difference scheme is as following: , h and k are grid sizes. The set of grid points in the xy plane is given by ( To fix an arbitrary constant associated with the solution of the pressure, without losing generality, we letp 1,1 = 0 and ignore the first discretized conti-

35
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
Here, n = (M −1)×N +M ×(N −1) and m = M ×N −1. Following the approach in [28,29], do QR decomposition of W T : where Q is n × n orthogonal, Q 1 is n × m, Q 2 is n × (n − m), R is n × m and R 1 (which is composed of the first m rows of R) is m × m nonsingular and upper triangular. Eliminatingp using Wv = 0, we thus essentially obtain the . Therefore, the original eigenvalue problem has precisely n − m eigenvalues, which can be obtained by solving the eigenvalues of the matrix

36
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. The approximation of (41) and (42) is This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0051846 are interpreted as and the initial condition (39) is interpreted as where h and k are the dimensions of the grids, and τ is the time increment. )k) at the time t * = (n +1)τ , where n = 0, 1, 2, · · · . We impose p n+1 1,1 = 0 as before to fix the arbitrary constant associated with the pressure solution. In the meantime, the discretized continuity equations for i = 1, j = 1 are ignored. Then the unknown values at time n+1 are uniquely determined and can be solved step by step.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0051846