Multiscale numerical modelling of perfusion in deformable double porous media described by the Biot-Darcy-Brinkman model¶

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. Two small scale parameters $\veps$ and $\delta$ chracterize micro- and meso-porosities.

At the mesoscopic scale, the periodic structure is formed by two systems of fluid filled channels occupying domains $\Om_\alpha^\delta, \alpha=1,2$ , whereby $\Om_c^\delta\cup_\alpha$ 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 mesoscopic channels $\Om_\alpha^\delta, \alpha=1,2$ are mutually separated by the micro-porous matrix $\Om_m^\delta$ , so that also two disconnected interfaces are defined $\Gamma_\alpha^\delta=\partial\Om_\alpha^\delta\cap\Om_m^\delta$ .

The microporous material is generated as a periodic lattice by repeating representative volume element (RVE) occupying domain $Y^\veps=\veps Y$ . It splits in the solid $Y_s$ and fluid part $Y_f=Y\setminus Y_s$ with interface $\Gamma_Y$ , as seen on Fig. 1.

Similarly, at the mesoscopic level the porous structure is generated by RVE $Z^\delta=\delta Z$ , with decomposition into the part filled by microporous matrix $Z_m$ , and fluid channels $Z_c=Z\setminus Z_m$ . The channels further split into two sub-parts $Z_\alpha$ , each with interface $\Gamma^\delta_\alpha=\overline{Z}_m\cap\overline{Z}_\alpha$ , see Fig. 1.

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.

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 [RohanTurjanicovaLukes2020]).

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_\alpha$ at mesoscopic level. By similar 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. However, when compared to the model of a mesoscopic structure containing only one system of channels, the upscaled mesoscopic structures involving two mesoscopic channels yields two macroscopic velocity fields $\wb^{0,\alpha}, \alpha=1,2$ which describe the two parallel flows.

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,\alpha}, \alpha=1,2$ 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,\alpha} \,\dV &= {1\over \vert Z\vert}\int_{\Gamma^Z_\alpha} q \nb^{[c]}_k\, \dS .$

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

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

The mesoscopic sub-problems are solved with the periodic boundary conditions and $\Gamma^Z_{\alpha} = \overline{Z_m} \cap \overline{Z_\alpha}$ is the interface between the matrix part $Z_m$ and system of channels $Z_\alpha$ .

With the characteristic responses obtained by solving local sub-problems, the homogenized material coefficients $\Acalb$ , $\Bcalb$ , $\Hcalb^{\alpha\beta}$ , $\Pcalb^\alpha$ , $\Scalb^\alpha$ , 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^{\alpha\beta}_{ij} & = {1\over \vert Z\vert} \int_{Z_m} \nabla_z\vphi^{j,\beta}\cdot\Kb\nabla_z\vphi^{i,\alpha}\,\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^\alpha_{ij} &= \phi_\alpha\delta_{ij}-{1\over \vert Z\vert}\int_{\Gamma_Z} \pi^i n_j^{\alpha}\, \dS =\phi_\alpha\delta_{ij}+{1\over \vert Z\vert} \int_{Z_m} \nabla_z\pi^i\cdot\Kb\nabla_z\vphi^{j,\alpha}\,\dV,\\ \Scal^\alpha_{{ijkl}} &={1\over \vert Z\vert}\left[2\eta\int_{Z_\alpha}\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$ , pressure $p^0$ and velocity fields $\wb^{0,\alpha}, \alpha=1,2$ such that for all $\vb$ , $q$ and $\psib$

(10)¶ $\int_{\Omega} \left[\Acalb \eebx{\ub^0} - p^0 \Bcalb^T\right]: \eebx{\vb}\,\dV - \sum\limits_\alpha \int_{\Omega}\vb\cdot \left[(\Pcalb^\alpha)^T(\nabla_xp^0-\fb^f) + \sum\limits_{\beta=1,2} \Hcalb^{\alpha\beta}\wb^{0,\beta}\right]&\,\dV +\int_{\partial_\sigma\Omega}\bar{\phi}_c\bar{p^0}\nb\cdot\vb^0\,\dV\\ &= \int_{\Omega} \fb^{blk}\cdot\vb \,\dV+\int_{\partial_\sigma\Omega}\bar{\phi}_m\bar{\phi}_s\gb\cdot\vb\,\dV, + \int_{\Omega} \fb^f \cdot \vb \dV,\\ \int_{\Omega} q\left[\Bcalb:\eebx{\dot\ub^0} + \dot p^0 \Mcal\right] \dV +\int_\Omega \nabla q\cdot \left[\Kcalb(\nabla_xp^0-\fb^f)-\sum\limits_{\beta=1,2}\Pcalb^\beta\wb^{0,\beta}\right]\,\dV &= -\int_{\partial_w\Omega} q (\bar\phi_m\bar\phi_f\bar w_n^{mic}+\sum\limits_{\beta=1,2}\bar\phi_\beta\bar w_n^{mes,\beta}) \,\dV,\\ \int_{\Omega}\eebx{\psib}:\Scalb^\alpha\eebx{\wb^{0,\alpha}+\dot\ub^0}\,\dV +\int_\Omega\psib\cdot\left[(\Pcalb^\alpha)^T(\nabla_xp^0-\fb^f) + \sum\limits_{\beta=1,2} \Hcalb^{\alpha\beta}\wb^{0,\beta} \right]\,\dV&=0.$

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\left[\Acalb \eebx{\ub^0} - p^0 \Bcalb^T\right]&=-\bar{\phi}_c\bar{p^0}\nb+\bar{\phi}_m\bar{\phi}_s\gb=\0b \; \textrm{ on }\partial_\sigma\Om,\\ \nb\cdot\left[\Kcalb(\nabla_x p^0-\fb^f)-\sum\limits_{\beta=1,2}\Pcalb^\beta \wb^{0,\beta}\right]&= \bar\phi_m\bar\phi_f\bar w_n^{mic}+\sum\limits_{\beta=1,2}\bar\phi_\beta\bar w_n^{mes,\beta} \; \textrm{ on }\partial_j\Om,\\ \nb\cdot\Scalb^\alpha\eebx{\wb^{0,\alpha}+\dot\ub^0}&=\bar\phi_\alpha P^0_\alpha\nb=\0b\; \textrm{ on }\partial_p\Om.$

For the purpose of this example, we simplify the problem (13). We 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$ , velocity $\wb^{0,\beta}$ and pressure fields $p^0$ such that for all $\vb^0$ , $q^0$ and $\psib^0$

(13)¶ $\int_{\Omega} \left[\Acalb \eebx{\ub^0} - p^0 \Bcalb^T\right]: \eebx{\vb}\,\dV - \sum\limits_\alpha \int_{\Omega}\vb\cdot \left[(\Pcalb^\alpha)^T(\nabla_xp^0-\fb^f) + \sum\limits_{\beta=1,2} \Hcalb^{\alpha\beta}\wb^{0,\beta}\right]\,\dV &= 0,\\ \int_{\Omega} q\left[\Bcalb:\eebx{\dot\ub^0} + \dot p^0 \Mcal\right] \dV +\int_\Omega \nabla q\cdot \left[\Kcalb(\nabla_xp^0-\fb^f)-\sum\limits_{\beta=1,2}\Pcalb^\beta\wb^{0,\beta}\right]\,\dV &= -\int_{\partial_w\Omega} q \sum\limits_{\beta=1,2}\bar\phi_\beta\bar w_n^{mes,\beta} \,\dV,\\ \int_{\Omega}\eebx{\psib}:\Scalb^\alpha\eebx{\wb^{0,\alpha}+\dot\ub^0}\,\dV +\int_\Omega\psib\cdot\left[(\Pcalb^\alpha)^T(\nabla_xp^0-\fb^f) + \sum\limits_{\beta=1,2} \Hcalb^{\alpha\beta}\wb^{0,\beta} \right]\,\dV&=0.$

The initial conditions of fields $\ub^0$ and $p^0$ are necessary for computation of time dependent problem and are computed using the steady state form of the boundary problem defined above.

Numerical simulation¶

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_BD2B-1/perf_BD2B_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.

Then, the script computes both the steady state and time evolution of the problem (13). For steady state, the macroscopic sample is fixed on both top and bottom 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 velocities $\bar\wb^\alpha_n, \alpha=1,2$ , see (13).

For simulation of the time evolution of the macroscopic problem we take the steady state as initial value of displacement and preassure fields. In this case, the gradually increasing displacement is prescribed on the top side of the macroscopic sample.

Fig. 3 Boundary conditions applied at the macroscopic level

The resulting macroscopic pressure field $p$ , displacement $\ub$ and the velocity fields $\wb^\alpha, \alpha=1,2$ are depicted in Fig. 4, Fig. 5 and Fig. 4. These figures depict distribution of macroscopic quantities on the deformingin at macroscopic specimen at different computational times ( $t=0.0\,s$ , $t=0.45\,s$ and $t=T=1.0\,s$ ). The nondeformed shape of macroscopic specimen is visualised by its outline.

Fig. 4 Macroscopic sample and the resulting macroscopic fields: left - displacement $\ub$ field at different computational times; right - pressure $p$ field at different computational times.

Fig. 5 Macroscopic sample and the resulting macroscopic fields: left - displacement $\ub$ field in $xy$ -crosssection at different computational times; right - pressure $p$ field in $xy$ -crosssection at different computational times.

Fig. 6 Macroscopic sample and the resulting macroscopic fields: left - 3D view of velocity field $\wb_1$ in $xy$ -crosssection at different computational times; right - velocity field $\wb_2$ in $xy$ -crosssection at different computational times.

References¶

[RohanTurjanicovaLukes2020]

Rohan E., Turjanicová J., Lukeš V. Multiscale modelling and simulations of tissue perfusion using the Biot-Darcy-Brinkman model. Submitted to Computers & Structures, 2020

[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