The poroelastic theory considers the system formed by an elastic material with pores that can be penetrated by a fluid and/or dissolved solutes. A typical biological example is the glucose solute penetrating a tissue layer. The basic theory of poroelasticity was developed by Biot [1] and its state-of-the-art can be found in the well-known books [1,2,3] (see also recent studies, e.g., [4,5,6,7,8]). A poroelastic material is considered as the superposition of two continuous media: the matrix (skeleton), occupying the fractional volume ${\theta }_M$, and the system of pores saturated by a fluid, occupying the fractional volume ${\theta }_F({\theta }_F+{\theta }_M=1$). The deformation of the system under the fluid pressure is described by a deformation vector, and the dynamics of the deformation under the forces needs in general to be described by second order tensors. The relationship between stress and strain is usually assumed to be linear. The flux of fluid depends on the hydrostatic pressure and solute gradient (called osmotic pressure). The diffusive and convective transport mechanisms should be taking into account as well. The general three-dimensional theory describing transport in poroelastic materials (PEM) is very complex because relevant mathematical models involve 3D nonlinear partial differential equations (PDEs). As a result, in order to obtain analytical results, one-dimensional versions are usually discussed [9,10,11]. Here we study a one-dimensional model for fluid and solute transport in poroelastic materials(PEM) that was derived in [10] and generalized [11]. The governing equations of the model read as where ${\sigma }_1=\sigma RT$ and ${\rho }_F={\rho }_F^0$ are positive constant and the lower subscripts t and x denote differentiation with respect to these variables. The physical/biological meanings of the notations used above are presented in Table 1.

The nonlinear system of PDEs (1)–(5) was integrated in the stationary case. As a result, all steady-state solutions were identified and examples of their application for the glucose fluid transport in a biological tissue was provided [10]. In the nonstationary case, system (1)–(5) is not integrable, therefore the classical Lie method [12,13,14,15,16] was adopted for search exact solutions. Nowadays this method is widely used for construction of exact solutions of nonlinear PDEs arising in real-world applications and the most remarkable works are cited in the above cited books. However, it can be easily noted that there are not many studies devoted to multicomponent systems of PDEs because essential technical difficulties occur if one intends to find exact solutions for such systems by applying Lie symmetries. Taking into account the above observation, we refer the reader to the recent works [17,18,19,20,21,22], devoted to applications of the classical Lie method to the nonlinear three-component systems of PDEs.

Although several nontrivial Lie symmetries were identified and successfully applied for finding exact solutions in [10] (see also generalisations in [11]), some limiting cases have been not analysed therein. In particular, the special cases $D=0$ and/or $S=0$ were not examined in [10,11]. It is proved here that the special cases listed above lead to a rich Lie symmetry involving symmetry operators that do not occur for system (1)–(5) with $DS\ne 0$. It should be stressed that the governing equations with $D=0$ and/or $S=0$ can still be used as a real-world model in some cases. For example, it is known that the diffusion term in equation (5) is negligible (i.e. $D=0$) if molecules of large size (because of high atomic weight) are dissolved. A typical example is albumin, for which the diffusivity is circa 50 times smaller than the glucose diffusivity [23].

This paper is organised as follows. In Section 2, all possible extensions of the Lie algebra of invariance of the nonlinear system of PDEs (1)–(5) are derived. It is proved that there are four inequivalent cases when the system admits additional Lie symmetries depending on the values of the parameters D and S. Section 3 is devoted to the constructions of exact solutions of system (1)–(5) with $S=0$. Because the system admits an additional Lie symmetry, new reductions to systems of ODEs and, as a result, new exact solutions of the nonlinear model in question are derived. In Section 4, an exact solution is analysed in order to show its applicability for modeling solute transport in PEM. The analysis is supported by 3D plots of the solution. In Section 5, new exact solutions of the nonlinear system (1)–(5) with arbitrary D and S are constructed. The solution obtained were not identified in the previous study [10]. In Section 6, we present conclusions highlighting the main results obtained.

Here we start from the governing equations (1)–(5) of the model for fluid and solute transport in poroelastic materials that was derived in [10]. By the application of the substitution where $p^{*}$ is a so-called effective pressure, the above system takes the form In (7), parameters $k,{\lambda }^{*},D$, and S are constants satisfying the restrictions

According to the standard terminology, Lie algebra (9) is called the principal algebra of system (7) (see, e.g., Chapter 1 in [16]). Because the latter involves several parameters, a Lie symmetry classification problem arises that was not solved in [10].

We remind the reader that Lie symmetry classification (LSC) problems (group classification problems) are the most difficult those in Lie symmetry analysis. One may claim that Sophus Lie carried out some theoretical foundations for solving LSC problems, in particular he has done classification of Lie symmetries for a class of linear two-dimensional PDEs. During the last few decades, a vast number of papers were published devoted to theoretical foundations and applications of algorithms for classification of Lie symmetries of specific PDEs (systems of PDEs). Unfortunately, a misleading terminology is used in many of them because instead of a complete solving of the LSC problem for a given PDE, only examples of extensions of the principal algebra are presented. The current state-of-the-art and the relevant list of most important publications can be found in [16] (see Chapter 2 therein).

Proof of this theorem is based on the infinitesimal criteria of invariance, which was formulated by S. Lie in his pioneering works [12,13]. In the case of a system of PDEs of arbitrary order, this criteria can be found, e.g., in [14] (see Section 1.2.5). Note, system (7) consists of five PDEs of the second order. So, a linear first-order operator of the form (here the coefficients ${\xi }^0,{\xi }^1$ and ${\eta }^i,i=1,\cdots ,5,$ to-be-determined functions of independent and dependent variables) is a Lie symmetry (operator of Lie’s invariance, point symmetry) of system (7) provided the following equalities are simultaneously satisfied for each solution $\left(u(t,x),\rho (t,x),p^{*}(t,x),{\theta }_F(t,x),c(t,x)\right)$ of the PDE system (7). Here $X^{\left(2\right)}$ is the second-order prolongation of the operator X, which is again the first-order operator with coefficients defined by the well-known formulae via the first- and second-order derivatives of unknown coefficients ${\xi }^0,{\xi }^1$ and ${\eta }^i$ (see, e.g., [14], Section 1.2.1).

After relevant calculations, formulae (10) are reducible to a linear system of PDEs, called the system of determining equations (DEs), for finding the functions ${\xi }^0,{\xi }^1$ and ${\eta }^i$. Some of the equations belonging to the system of DEs, do not contain parameters S and D, therefore solving them, we obtain the most general form of the Lie symmetry operator where ${\alpha }_i(i=0,\cdots ,3)$ are arbitrary constants and $g\left(t\right)$ ia an arbitrary smooth function, while the constant ${\alpha }_4$ and the functions ${\eta }^4$ and ${\eta }^5$ must satisfy the restriction $S{\alpha }_4=0$ and the remaining DEs

If one assumes that all parameters arising in (7) are arbitrary constants then the result presented in Theorem 1 in a straightforward way is derived. However, here we need to determine all possible values of parameters (under restrictions (8)) when a wider Lie symmetry occurs. Looking at the system of DEs (12)–(15), one notes (see equations (12)) that two different cases, $S\ne 0$ and $S=0$, should be separately examined.

In the case$S\ne 0$, we immediately obtain ${\alpha }_4=0,{\eta }^4=0,{\eta }_{\rho }^5={\eta }_{{\theta }_F}^5=0.$ So, the system of determining equations (12)–(16) is reduced to the form If $S\ne 1$ then one readily finds ${\eta }^5={\alpha }_{5}c$, hence the principal algebra (9) is only obtained. However, assuming $S=1$, Cases 1 and 3 of Table 2 are obtained when $D\ne 0$ and $D=0$, respectively.

In the case$S=0$, an extension of the principal algebra (9) via the operator $x{\partial }_{p^{*}}+\frac{x^2}{2{\lambda }^{*}}{\partial }_u$ can be easily identified. If $D\ne 0$ then Case 2 in Table 2 is obtained.

Finally, special values $S=0$ and $D=0$ lead to the widest Lie symmetry. In fact, equations (12)–(14) simply disappear and the remaining equations are Solving equation (18), one finds Substituting the above function into equation (17), the function is obtained. Thus, the Lie symmetry operator (11) with the above functions ${\eta }^4$ and ${\eta }^5$ produce exactly the Lie symmetries listed in Case 4 of Table 2.

Although the cases listed in Table 2 are very special and are questionable from the applicability point of view, new Lie symmetries are very interesting. In particular, system (7) with $S=0$ admits the symmetry $x{\partial }_{p^{*}}+\frac{x^2}{2{\lambda }^{*}}{\partial }_u$ that can be not predicted using any physical laws. Here we study this case in detail and our goal is to construct a wide range of time-dependent solution.

As it was noted above, the most interesting cases from the symmetry point of view are the second and the fourth cases (see Table 2) because they involve the Lie symmetry $Y=x{\partial }_{p^{*}}+\frac{x^2}{2{\lambda }^{*}}{\partial }_u$ that cannot be predicted by any physical/biological consideration. The fourth case is rather artificial from applicability point of view. Indeed, the restrictions $D=S=0$ means that both diffusive and convective transports are neglected what is very strong assumption (however, one cannot claim that such a special case is absolutely unrealistic). Thus, we study the second case, $S=0$ in what follows.

The general technique for constructing exact solution using Lie symmetries is well-known from the theoretical point of view. There are two different approaches described, e.g., in Section 1.3 [16]. One of them is based on construction of optimal systems of inequivalent (nonconjugated) subalgebras of the given Lie algebra of symmetries. Although this technique is very popular, one successfully works only in the case of Lie algebras of low dimensionality (up to 4). In fact, all optimal systems for such algebras were described in a seminal work [24]. However, a pure algebraic problem occurs if one deals with a Lie algebra of the dimensionality 5 and higher. To the best of our knowledge, there is no generalisations of the results presented in [24] on higher-dimensional Lie algebras. Notably, the Lie algebra of symmetries of system (7) is infinity-dimensional (see the operator $g\left(t\right){\partial }_{p^{*}}$ in Theorem 1).

The second technique for constructing exact solution via Lie symmetries is based on the most general linear combination of basic operaintors of the Lie algebra in question. Of course, many technical difficulties arise in the case of higher-dimensional Lie algebras, however, these difficulties can be overcome by a careful analysis of the relevant invariant surface condition. Here this is demonstrated in the case of Lie algebra with the basic operators

Let us consider the most general linear combination of Lie symmetries listed above where ${\alpha }_1,{\alpha }_2$ and ${\beta }_i(i=0,1,2,3)$ are arbitrary constants, ${\alpha }_1^2+{\alpha }_2^2\ne 0$ (if ${\alpha }_1={\alpha }_2=0$ then there is no reduction to an ODE system).

Depending on the value of the constants ${\alpha }_1$ and ${\alpha }_2$, operator (19) produce three inequivalent ansätze for reduction of system (7) to ODE systems. Thus, one need to consider such cases :

$I.$${\alpha }_1{\alpha }_2\ne 0\Rightarrow {\alpha }_1=1,{\alpha }_2=\alpha \ne 0$;

$II.$${\alpha }_1\ne 0\Rightarrow {\alpha }_1=1,{\alpha }_2=0$;

$III.$${\alpha }_1=0,{\alpha }_2\ne 0\Rightarrow {\alpha }_2=1$.

From the very beginning, it can be noted that the function $g\left(t\right)$ does not play role in the first two cases and we may set $g\left(t\right)=0$ without losing a generality. In fact, ansätze obtained contain the added term $g^{*}\left(t\right)$ in expressions for the function $p^{*}$ and this term vanishes after substituting into system (7). So, one always can generalise an arbitrary solution by inserting into the pressure component the additional term $g^{*}\left(t\right)$ that simply follows from the Lie group generated by the operator $g\left(t\right){\partial }_{p^{*}}$: In the third case, the function $g\left(t\right)$ does play a role.

In the case I operator (19) leads to the ansatz where $f_i\left(\omega \right)(i=1,\cdots ,5)$ are new smooth functions.

Substituting ansatz (20) into PDE system (7), one arrives at the ODE system where new notations $F={f_1}^{\prime }$ and $\varphi ={f_3}^{\prime }$ are used.

We are looking for exact solutions of the ODE system (21) in the two essentially different subcases $\left(\mathbf{i}\right)$ $1+F=0$ and $\left(\mathbf{ii}\right)$$1+F\ne 0$.

Subcase$\left(\mathbf{i}\right)$. System (21) after the relevant calculations leads to the result where $C_i$ are arbitrary constants and constant $D>0$.

If ${\beta }_2=0$, then the general solution of equation (22) has the form $f_5=C_5+C_6\omega$. Thus, taking into account ansatz (20) and the above formulae, the exact solution of the PDE system (7) with $S=0$ (here ${\theta }_F$ is an arbitrary smooth function) is obtained.

If ${\beta }_2\ne 0$, then the general solution of equation (22) can be found in an explicit form only for some fixed function $f_4$. Let us solve equation (22) assuming $f_4\left(\omega \right)={\theta }_0{\omega }^n$ (here ${\theta }_0>0$ and n are arbitrary constants). So, one need to construct exact solutions of the equation Equation (26) has essentially different general solutions depending on the value of the constants ${\beta }_2$ and n. In the simplest case when $n=0$ the general solution of equation (26) has the form

The obtained functions $f_4={\theta }_0$ and $f_5$ from (27) lead to the components of the exact solution of the PDE system (7) with $S=0$ and $D>0$, while other components of the corresponding solution are presented in formulae (23)–(25).

Now let us return to equation (26). In the case $n\ne -2,$ this equation is reduced to the Bessel equation if ${\beta }_2<0$, and to the modified Bessel equation if ${\beta }_2>0$, by the substitution Thus, the general solution of equation (26) with $n\ne -2$ can be expressed in the following forms if ${\beta }_2<0$, and if ${\beta }_2>0$. Here J and Y are the Bessel functions, while I and K are the modified Bessel functions.

Note that solution (27) can be derived from (28)–(29) by setting $n=0$ in the last ones.

In the case $n=-2,$ the general solution of equation (26) has the form Note that in the case $n=-2$ and $4{\beta }_2{\theta }_0+D<0$ equation (26) has only complex solutions, not real ones.

The functions $f_5$ (from (28)–(30)) and $f_4={\theta }_0{\omega }^n$ produce the components and of the solutions of the PDE system (7) with $S=0$. The remaining components of the corresponding solutions are still given by formulae (23)–(25).

Subcase$\left(\mathit{ii}\right)$. In this case, the functions $f_1,f_2,f_3$ and $f_4$ can be expressed via the function F, namely : where the function F is a solution of the equation while the function $f_5$ must be founded from the equation

ODE (32) is a nonlinear Abel type equations and its full integration is questionable, moreover, even particular solutions are unknown. However, assuming that the function F is linear, one can find the exact solution provided the following restrictions are satisfied

Another possibility to integrate ODE (32) arises if one sets In this case ODE (32) takes the form The latter is integrable via the Lambert function: where The exact solution (36) can also be presented in the implicit form It can be easily seen that a special case occurs when $\frac{{\lambda }^{*}}{{\alpha }^2{\rho }_F^0}=1.$ In this case, one again arrives at the linear function Each of the exact solutions (34)–(37) can be used for finding the function $f_5$ by solving the linear ODE (33). Because the formulae obtained are very cumbersome (in particular, the Heun functions are obtained if one applies (37)), here we present only the details concerning solution (34). In this case, ODE (33) can be rewritten in the form Although this equation is very awkward, one may find the substitution transforming ODE (38) into the well-known Whittaker equation (see, e.g., [25]) where The general solution of equation (40) has the form where M and W are the Whittaker functions. In the case $\mu =\frac{1}{2}+\nu$ (that leads to the conditions ${\beta }_2=0$ or ${\beta }_2=C_0\alpha$), the function M takes the form $M_{\frac{1}{2}+\nu ,\nu }\left(\tau \right)=e^{-\frac{\tau }{2}}{\tau }^{\frac{1}{2}+\nu }.$ Thus, the general solution of the Whittaker equation in this special case has the form

Taking into account (39) and (42), we arrive at the solution of ODE (38) under restrictions ${\beta }_2=0$ and ${\beta }_2=C_0\alpha$, namely :

Using the obtained functions $f_5$ and the above formulae, one can construct exact solutions of the PDE system (7). In particular, setting $C_6=0$ in (43) for simplicity, an exact solution of the system in question can be written down in the form where

Similarly, setting $C_6=0$ in (44), another solution of the PDE system (7) with $S=0$ can be derived in the form while other components have the same form as in (45).

In the case $II$ operator (19) leads to the ansatz where $f_i\left(x\right)(i=1,\cdots ,5)$ are new smooth functions, while $P_2\left(x\right)=\frac{{\beta }_3}{2}{x}^2+{\beta }_{1}x+{\beta }_0.$

The special case is $P_2\left(x\right)=0,$ it means that the operator $x{\partial }_{p^{*}}+\frac{x^2}{2{\lambda }^{*}}{\partial }_u$ is not taken into account, hence one may apply ansatz (46) to system (7) with an arbitrary parameter S (see Section 5).

Here we substitute ansatz (20) with $P_2\ne 0$ into PDE system (7) with $S=0$ to obtain a reduced ODE system. As a result, one arrives at the system

The first four equations of the ODE system (47) can be easily integrated, producing the functions while the last equation on the function $f_5\left(x\right)$ takes the form

It is difficult to construct exact solutions of equation (49) without additional restrictions. Note that equation (49) in the case $C_4=0$ is the Heun equation and its general solution can be constructed using software such as Maple. Since this case leads to very cumbersome formulae, we omit those here.

Let us set ${\beta }_0={\beta }_1={\beta }_2=C_4=0,{\beta }_3\ne 0$, then equation (49) takes the form The latter is integrable in the terms of modified Bessel functions. Thus, we arrive at the exact solution of the PDE system (7) with $S=0$ :