Numerical simulation of viscous flow in deformable double porous media

Mathematical model

We consider a double porous medium which consists of an elastic solid matrix \Om_s perforated by a system of channels filled with an incompressible fluid \Om_f with interface \Gamma_{fs}=\overline{\Om_f\cap\Om_s}. These components are arranged in a periodic lattice at both the micro- and mesoscopic level. Thus, the porous matrix at the mesoscopic level can be generated medium can be generated by copies of the microscopic reference unit cell Y and, subsequently, macroscopic body can be generated as alattice of the mesoscopic reference unit cell Z, see Fig. 1. Two small scale parameters \veps and \delta chracterize micro- and meso-porosities. At the mesoscopic scale, the periodic structure is formed by fluid filled channels occupying domain Om_c^delta and by domain \Om_m^\delta=\Om\setminus\Om_c^\delta which is constituted by a microporous material. In particular, domain \Om_p^{\veps,\delta}\in\Om_m^\delta represents micro pores saturated by fluid, whereas \Om_s^{\veps,\delta}=\Om_m^\delta\setminus\Om_p^{\veps,\delta} is the skeleton, see Fig. 1..

To summarize the decompositions,

&\Om=\Om^\delta_c\cup\Om^\delta_m\cup\Gamma^\delta,\\
\mbox{solid }&\Om_s\in\Om_m^\delta,\\
\mbox{microporosity }&\Om_p^{\veps,\delta}=\Om_m^\delta\setminus\Om_s^{\veps,\delta},\\
\mbox{fluid }&\Om_f=\Om_c^\delta\cup\Om_p^{\veps,\delta}\cup\Gamma_{cp}^{\veps,\delta}.\\

The superscripts {}^\veps and {}^\delta denote the quantities oscillating within the heterogeneous structure with the period equal to the size of the micro- and mesoscopic periodic unit. However, we drop the superscript in following text to simplify the notation.

_images/fig-2level-homog-eps-del.JPG

Fig. 1 Macroscopic domain \Omega and decomposition of microscopic domain Y and mesoscopic domain Z.

The mechanical behavior of such a structure can be described using the two-level asymptotic homogenization method, (for more detailed explanaition we refer to [RohanTurjanicovaLukes2019]).

The mechanical properties of the deformable matrix are given by elasticity tensor \Db=(D_{ijkl}) which satisfies the usual symmetries.

The fluid is characterized by viscosity \eta which is given by a piece-wise constant function,

(1)\eta&=\veps^2\eta_p\, \mbox{ in }\Omega_p,\\
\eta&=\eta_c\, \mbox{ in }\Omega_c.

The scaling of the viscosity in micropores \Omega_p is the standart consequence of the non-slip boundary condition on he pore wall.

The problem of the fluid flow in deformable media at microscopic level is given by the following equilibrium equations and boundary conditions governing displacement of the solid \ub and both the fluid pressure and velocity fields (p,\vb^f):

(2)-\nabla\cdot \Db\eeb{\ub} & = \fb^s, \quad \mbox{ in } \Om_{s}, \\
\nb\cdot \Db\eeb{\ub} & =\nb\cdot\sigmab^f, \quad \mbox{ on } \Gamma_{fs}, \\
\nb\cdot \Db\eeb{\ub} & =\nb\cdot\gb^s, \quad \mbox{ on } \partial_\sigma\Omega_{s}, \\
\ub & =\0b, \quad \mbox{ on } \partial_u\Omega_{s}, \\
-\nabla\cdot(2\eta\eeb{\vb^f}-p\Ib) & =\fb^f, \quad \mbox{ in } \Omega_{f}, \\
\nabla\cdot\vb^f & =0, \quad \mbox{ in } \Omega_{f}, \\
\vb^f & =\dot\ub, \quad \mbox{ on } \Gamma_{fs}, \\
\vb^f-\dot{\tilde{\ub}} & =:\wb, \quad \mbox{ on } \partial_{ext}\Om_{f},

where \sigmab^f=\veps^2\eta_p\eeb{\vb^f}-p\Ib is the fluid stress, \eeb{\ub}=(e_{ij}(\ub)) is the strain in the solid with components e_{ij}(\ub)=1/2(\partial_ju_i+\partial_iu_j), \fb^{s,f} denotes the volume forces in the solid or in the fluid, and \gb^{s,f} is the surfacetraction stresses acting on the solid part. The relative fluid velocity \wb = \vb^f - \dot{\tilde{\ub}} in the fluid-filled pores \Om_f is defined whit use of a smooth extention \tilde{\ub} of the dislacement field \ub from solid \Om_s to whole domain \Om.

Two-level homogenization

Due to the double-porous nature of the medium, we performe two levels of homogenization. The 1st-level of homogenization concerns the asymptotic analysis \veps \longrightarrow 0 related to the fluid-structure interaction in microporous structure situated in \Omega_m. We apply the standard homogenization techniques to the above problem. It results in the limit model for \veps \longrightarrow 0, where \veps is the scale parameter relating the microscopic and macroscopic length scales. The homogenization process leads to local microscopic problems, defined within a reference periodic cell Y , and to the mesoscopic problem describing the behavior of the homogenized matrix at the mesoscopic level. The mesoscopic problem involves the homogenized material coefficients which are evaluated using the solutions of the local problems. The 2nd-level of homogenization deals with upsacling from meso- to macroscopic scale. It results in the limit model for \delta \longrightarrow 0 and subsequently in the local mesoscopic problems on a reference periodic cell Z, and in the derivation of the homogenized problem at macroscopic level. Due to linearity of the problem, the microscopic, mesoscopic and macroscopic problems are decoupled.

The local microscopic responses are given by the following sub-problems which are solved within the periodic reference cell Y, see Fig. 1, that is decomposed similarly to the decomposition of domain \Omega:

  • Find \omegab^{ij}, \omega^{P} such that for all \vb for any i, j = 1, 2, 3

(3)\int_{Y_s} \Db \eeby{\omegab^{ij} + \Pib^{ij}}: \eeby{\vb}\,\dV   &= 0, \\
\int_{Y_{s}} \Db \eeby{\omegab^P}: \eeby{\vb}\,\dV &=
-{1\over \vert Y\vert}\int_{\Gamma_Y} \vb \cdot \nb^{[s]}\, \dS,

where \Pi^{ij}_k = y_j \delta_{ik}.

  • Find \psib^{k}, {\pi}^{k} such that for all \vb, q satisfying

(4)\int_{Y_f} \eta_p \nabla_y{\psib^{k}}: \nabla_y{\vb}\,\dV - \int_{Y_f} \pi^{k}\nabla_y\cdot\vb \,\dV &=
-\int_{Y_f} v_k\,\dV , \\
\int_{Y_f} q\nabla_y\cdot\psib^{k} \,\dV &= 0.

The microscopic sub-problems are solved with the periodic boundary conditions and \Gamma_Y=\overline{Y_s} \cap \overline{Y_f} is the interface between solid and fluid part of the cell Y.

With the characteristic responses obtained by solving local sub-problems, the homogenized material coefficients \Ab, \Bb, \Kb and M can be evaluated using the following expressions:

(5)A_{ijkl} & = {1\over \vert Y\vert} \left[ \int_{Y_{s}} \Db \eeby{\omegab^{kl} + \Pib^{kl}}: \eebz{\omegab^{ij} + \Pib^{ij}}\,\dV\right],\\
B_{ij} & = \phi_f\delta_{ij}-{1\over \vert Y\vert} \left[\int_{Y_s} \Db \eeby{\omegab^{P}}:\eeby{\Pib^{ij}}\,\dV\right],\\
K_{ij} & = {1\over \vert Y\vert} \int_{Y_f} \nabla_y\psi^i:\nabla_y\psi^i\,\dV,\\
M & = {1\over \vert Y\vert} \left[\int_{Y_s} \Db \eeby{\omegab^{P}}:\eeby{\omegab^{P}} \,\dV \right].

Homogenization - 2nd level

At the 2st-level of homogenization, the asymptotic analysis \delta \longrightarrow 0 is related to the interaction between the homogenized microporous matrix in \Omega_m and fluid in channels \Omega_c at mesoscopic level. By same upscaling procedure as in 1st level of homogenization, we obtain local mesoscopic problems, defined within a reference periodic cell Z, where we enter the homogenized coefficients obtained by 1st level homogenization. We also arrive to the global problem describing the behavior of the homogenized matrix at the macroscopic level. The homogenized material coefficients describing whichdescribe behavior at macroscopic level are evaluated using the solutions of the local mesoscopic problems.Due to linearity of the problem, the microscopic, mesoscopic and macroscopic problems are decoupled.

The local mesoscopic responses are given by the following sub-problems which are solved within the periodic reference cell Z, see Fig. 1, that is decomposed similarly to the decomposition of domain \Omega:

  • Find \omegab^{ij}, \omega^{P} such that for all \vb for any i, j = 1, 2, 3

(6)\int_{Z_m} \Ab \eebz{\omegab^{ij} + \Pib^{ij}}: \eebz{\vb}\,\dV   &= 0, \\
\int_{Z_{m}} \Ab \eebz{\omegab^P}: \eebz{\vb}\,\dV - \int_{Z_m} \Bb: \eebz{\vb} \,\dV &=
 -{1\over \vert Z\vert}\int_{\Gamma_z} \vb \cdot \nb^{[m]}\, \dS ,

where \Pi^{ij}_k = z_j \delta_{ik}.

  • Find \pi^k, \vphi^k such that for all q satisfying

(7)\int_{Z_{m}} \nabla_z q \Kb \nabla_z \pi^k \,\dV &= -\int_{Z_{m}} \nabla_z q \Kb \nabla_z z_k , \\
\int_{Z_{m}} \nabla_z q \Kb \nabla_z \vphi^k \,\dV &= -{1\over \vert Z\vert}\int_{\Gamma_z} q \nb^{[c]}_k\, \dS .

  • Find \psib^{ij}, \hat{\vphi}^{ij} such that for all \vb, q satisfying

(8)\int_{Z_c} 2\eta \eebz{\psib^{ij}}: \eebz{\vb}\,\dV - \int_{Z_c} \hat{\vphi}^{ij}\nabla_z\cdot\vb \,\dV &=
-\int_{Z_{c}} 2\eta \eebz{\Pib^{ij}}: \eebz{\vb}\,\dV , \\
\int_{Z_c} \nabla_z\cdot\psib^{ij} q \,\dV &= -\int_{Z_c} \nabla_z\cdot\Pib^{ij} q \,\dV.

The mesoscopic sub-problems are solved with the periodic boundary conditions and \Gamma_{z} = \overline{Z_m} \cap \overline{Z_c} is the interface between the matrix part Z_m and canal Z_c.

With the characteristic responses obtained by solving local sub-problems, the homogenized material coefficients \Acalb, \Bcalb, \Hcalb, \Pcalb, \Scalb, and \Mcal can be evaluated using the following expressions:

(9)\Acal_{ijkl} & = {1\over \vert Z\vert} \left[ \int_{Z_{m}} \left[\Ab \eebz{\omegab^{kl}
+ \Pib^{kl}}\right]: \eebz{\omegab^{ij} + \Pib^{ij}}\,\dV\right],\\
\Bcal_{ij} & = \phi_c\delta_ij+{1\over \vert Z\vert} \left[\int_{Z_m} \Bb: \eebz{\Pib^{ij}+\omegab^{ij}} \,\dV
- \int_{\Gamma_Z} \nb^{m}\cdot \omegab^{ij} \,\dV\right],\\
\Hcal_{ij} & = {1\over \vert Z\vert} \int_{Z_m} \nabla_z\vphi^j\cdot\Kb\nabla_z\vphi^i\,\dV,\\
\Kcal_{ij} & = {1\over \vert Z\vert} \int_{Z_m} \nabla_z(z_j+\pi^j)\cdot\Kb\nabla_z(z_i+\pi^i)\,\dV,\\
\Pcal_{ij} &= \phi_c\delta_{ij}-{1\over \vert Z\vert}\int_{\Gamma_Z} \pi^i n_j^{m}\, \dS
=\phi_c\delta_{ij}+{1\over \vert Z\vert} \int_{Z_m} \nabla_z\vphi^j\cdot\Kb\nabla_z\vphi^i\,\dV,\\
\Scal_{{ijkl}} &={1\over \vert Z\vert}\left[2\eta\int_{Z_c}\eebz{\psib^{kl}+\Pib^{kl}}:\eebz{\psib^{ij}
+\Pib^{ij}}\,\dV-\int_{Z_c} \hat{\vphi}^{kl}\nabla_z\cdot\Pib^{ij}\, \dV\right],\\
\Mcal &={1\over \vert Z\vert}\left[\int_{Z_m}M\,\dV+ \int_{Z_m} \Bb: \eebz{\omegab^{P}} \,\dV
- \int_{\Gamma_Z} \nb^{c}\cdot \omegab^{P} \,\dV\right].

All coefficients are symmetric with respect to indices related to strain and strain rate tensors, i.e. a_{ij}=a_{ji}.

The global macroscopic problem is defined in terms of the homogenized coefficients as: Find the macroscopic displacements \ub^0 and p^0 such that for all \vb^0 and q^0

(10)\int_{\Omega} [\Acalb \eebx{\ub^0} - p^0 \Bcalb^T]: \eebx{\vb^0}\,\dV - \int_{\Omega}\vb^0\cdot\Pcalb^T(\nabla_xp^0-\fb^f) &
- \int_\Omega \vb^0 \Hcalb^T\wb^0\,\dV +\int_{\partial_\sigma\Omega}\bar{\phi}_c\bar{p^0}\nb\cdot\vb^0\,\dV\\
&=  \int_{\Omega} \fb^{blk}\cdot\vb^0 \,\dV+\int_{\partial_\sigma\Omega}\bar{\phi}_m\bar{\phi}_s\gb\cdot\vb^0\,\dV,
+ \int_{\Omega} f \cdot \vb \dV,\\
\int_{\Omega} q^0(\Bcalb:\eebx{\dot\ub^0} + \dot p^0 \Mcal) \dV +\int_\Omega (\Kcalb(\nabla_xp^0-\fb^f)-\Pcalb\wb^0)\,\dV
&= -\int_{\partial_w\Omega} q^0 (\bar\phi_m\bar\phi_f\bar w_n^{mic}+\bar\phi_c\bar w_n^{mes}) \,\dV,\\
\int_{\Omega}\eebx{\thetab^0}:\Scalb\eebx{\wb^0+\dot\ub^0}\,\dV+\int_\Omega\thetab^0\cdot\Pcalb(\nabla_xp^0-\fb^f)\,\dV
+\int_\Omega\Hcalb\wb^0 &=\int_{\partial_p\Omega}\bar\phi_c\bar p^0\nb\cdot\thetab^0\,\dV.

The Dirichlet boundary conditions prescribed for \ub^0 and \wb^0 can be imposed,

(11)\ub^0&=\bar{\ub}^0\; \textrm{ on }\partial_u\Om,\\
\wb^0&=\bar{\wb}^0\; \textrm{ on }\partial_w\Om.

The pressure p^0 fulfils zero-means conditions \int_\Om p^0=0 in the whole domain \Om.

The complementary Neumann-type boundary condition are specfied as follows

(12)\nb\cdot(\Acalb \eebx{\ub^0} - p^0 \Bcalb^T)&=-\bar{\phi}_c\bar{p^0}\nb+\bar{\phi}_m\bar{\phi}_s\gb=\0b
\; \textrm{ on }\partial_\sigma\Om,\\
\nb\cdot(\Kcalb(\nabla_xp^0-\fb^f)-\Pcalb\wb^0)&=\bar\phi_m\bar\phi_f\bar w_n^{mic}+\bar\phi_c\bar w_n^{mes}
\; \textrm{ on }\partial_j\Om,\\
\nb\cdot\Scalb\eebx{\wb^0+\dot\ub^0}&=\bar\phi_c\bar{p}^0\nb=\0b\; \textrm{ on }\partial_p\Om.

For the purpose of this example, we simplify the problem (13). We consider only steady state of the problem and omitt all volume froces \fb^f and \fb^{blk} and surface tractions \gb. Also we will consider closed micropores on the whole boundary \partial\Om, i.e. w_n^{mic}=0. The problem (10) becomes: Find the macroscopic displacements \ub^0 and p^0 such that for all \vb^0 and q^0

(13)\int_{\Omega} [\Acalb \eebx{\ub^0} - p^0 \Bcalb^T]: \eebx{\vb^0}\,\dV - \int_{\Omega}\vb^0\cdot\Pcalb^T\nabla_xp^0
- \int_\Omega \vb^0 \Hcalb^T\wb^0\,\dV&= 0\\
\int_\Omega (\Kcalb\nabla_xp^0-\Pcalb\wb^0)\,\dV
&= -\int_{\partial_w\Omega} q^0 \bar\phi_c\bar w_n^{mes} \,\dV,\\
\int_{\Omega}\eebx{\thetab^0}:\Scalb\eebx{\wb^0}\,\dV+\int_\Omega\thetab^0\cdot\Pcalb\nabla_xp^0\,\dV+\int_\Omega\Hcalb\wb^0 &=0.

Numerical simulation

_images/comp_cells.png

Fig. 2 Left - geometric representation of microscopic domain; Right - geometric representation of mesoscopic domain

To run the numerical simulation, download the archive, unpack it in the main SfePy directory and type:

./simple.py example_perfusion_BDB-1/perf_BDB_mac.py

This invoke the simply.py script which calculates the macroscopic problem (13) and calls the homogenization engine that solves the local subproblems for given parameter \veps, viscosity \eta and elastic tensor \Db. First, it solves subproblems (3) and (4) on microscopic cell Y (see Fig. 2 left) and evaluates the homogenized coefficients (5). Then, using solution from previous step, it solves subproblems (6)(8) on mesoscopic cell Z (see Fig. 2 right) and evaluates the homogenized coefficients (9). See [CimrmanLukesRohan2019] for more details related to the SfePy homogenization engine.

The macroscopic sample is fixed on both left and right side, so that no displacements are allowed, see Fig. 3. The defromation is induced due to the flow through porous matrix, as the responce to the prescribed velocity \bar\wb_n, see (13).

_images/bcond.png

Fig. 3 Boundary conditions applied at the macroscopic level

The resultin macroscopic pressure field p, displacement \ub and the velocity field \wb are depicted in Fig. 4, where deformation is visualised by deformed wireframe.

_images/results_macro.png

Fig. 4 Macroscopic sample and the resulting macroscopic fields: left - 3D view of pressure p, displacement \ub and velocity field \wb; right - pressure p, displacement \ub and velocity field \wb in xy-crosssection.

References

RohanTurjanicovaLukes2019

Rohan E., Turjanicová J., Lukeš V. The Biot–Darcy–Brinkman model of flow in deformable double porous media; homogenization and numerical modelling. Computers and Mathematics with applications, 78(9):3044-3066, 2019, DOI:10.1016/j.camwa.2019.04.004

CimrmanLukesRohan2019

Cimrman R., Lukes V., Rohan E. Multiscale finite element calculations in Python using SfePy. Advances in Computational Mathematics, 45(4):1897-1921, 2019, DOI:10.1007/s10444-019-09666-0