Numerical simulation of fluid-saturated piezoelectric porous media¶

Mathematical model¶

We consider a porous piezoelectric medium which consists of a piezoelectric matrix, embedded metallic electrodes (conductors) and fluid-filled inclusions. These components are arranged in a periodic lattice so that the medium can be generated by copies of the reference unit cell, see Fig. 1. The mechanical behavior of such a structure can be described using the two-scale asymptotic homogenization method, see [RohanLukes2018].

The mechanical properties of the piezoelectric solid are given by the following constitutive equations

(1)¶ $\sigma^\veps_{ij}(\ub^\veps,\vphi^\veps) & = A_{ijkl}^\veps e_{kl}(\ub^\veps) - g_{kij}^\veps E_k(\vphi^\veps), \\ D_k^\veps(\ub^\veps,\vphi^\veps) &= g_{kij}^\veps e_{ij}(\ub^\veps) + d_{kl}^\veps E_l(\vphi^\veps),$

where $\sigmab^\veps$ is the Cauchy stress tensor, $\vec{D}^\veps$ is the electric displacement, $\eb(\ub^\veps) = {1\over2}\left(\nabla\ub^\veps + (\nabla\ub^\veps)^T\right)$ is the strain tensor, $\ub^\veps$ is the displacement field, $\vec{E}(\vphi^\veps) = \nabla \vphi^\veps$ is the electric field and $\vphi^\veps$ is the electric potential. On the right hand side of (1), we have the fourth-order elastic tensor $A^\veps_{ijkl}$ ( $A^\veps_{ijkl} = A^\veps_{klij} = A^\veps_{jilk}$ ), the third-order tensor $g^\veps_{kij}$ ( $g^\veps_{kij} = g^\veps_{kji}$ ), which couples mechanical and electric quantities, and the permeability tensor $d^\veps_{kl}$ . The superscript ${}^\veps$ denotes the quantities oscillating within the heterogeneous structure with the period equal to the size of the periodic unit.

The quasi-static problem of the piezoelectric medium is given by the following equilibrium equations

(2)¶ $-\nabla\cdot \sigmab^\veps(\ub^\veps,\vphi^\veps) & = \fb^\veps, \quad \mbox{ in } \Om_{m*}^\veps, \\ -\nabla\cdot \vec D^\veps(\ub^\veps,\vphi^\veps) & = q_E^\veps,\quad \mbox{ in } \Om_m^\veps,$

end by the mass conservation equation for the k th fluid inclusion occupying domain $\Om_c^{k,\veps}$

(3)¶ $\int_{\partial \Om_{k,\veps}} {\ub^\veps}\cdot \nb \dS + \gamma p^{k,\veps} |\Om_c^{k,\veps}|= 0, \quad \forall k \in \{1,\dots,\bar k^\veps\}.$

where $p^{k,\veps}$ is the inclusion pressure and $\gamma$ is the fluid compressibility.

Two-scale homogenization¶

Fig. 1 Macroscopic domain $\Omega$ and decomposition of microscopic domain $Y$ ¶

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, and to the global problem describing the behavior of the homogenized medium at the macroscopic level. The global problem involves the homogenized material coefficients which are evaluated using the solutions of the local problems. Due to linearity of the problem, the microscopic and macroscopic problems are decoupled.

As we assume given potentials $\bar\vphi^k$ in each of the electrode networks, the dielectric properties must be appropriately rescaled in order to preserve the finite electric field in the limit:

(4)¶ $\gb^\veps = \veps \bar\gb,\\ \db^\veps = \veps^2 \bar\db.$

The local microscopic responses of the piezoelectric structure 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}$ , $\eta^{ij}$ such that for all $\vb$ , $\psi$ and for any $i, j = 1, 2, 3$

(5)¶ $\int_{Y_{m*}} \left[\Ab \eeb{\omegab^{ij} + \Pib^{ij}}\right]: \eeb{\vb}\,\dV - \int_{Y_m} \left[\bar\gb^T\cdot\nabla \eta^{ij}\right]: \eeb{\vb} \, \dV &= 0, \\ \int_{Y_m} \left[\bar\gb:\eeb{\omegab^{ij} + \Pib^{ij}} + \bar\db \cdot\nabla \eta^{ij}\right]\cdot\nabla\psi \, \dV &= 0,$

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

Find $\omegab^P$ , $\eta^P$ such that for all $\vb$ , $\psi$ satisfying

(6)¶ $\int_{Y_{m*}} \left[\Ab \eeb{\omegab^P}\right]: \eeb{\vb}\,\dV - \int_{Y_m} \left[\bar\gb^T\cdot\nabla \eta^P\right]: \eeb{\vb} \,\dV &= -{1\over \vert Y\vert}\int_{\Gamma_c} \vb \cdot \nb\, \dS , \\ \int_{Y_m} \left[\bar\gb:\eeb{\omegab^P} + \bar\db \cdot\nabla \eta^P\right]\cdot\nabla\psi \,\dV &= 0.$

Find $\omegab^\rho$ , $\eta^\rho$ such that for all $\vb$ , $\psi$ satisfying

(7)¶ $\int_{Y_{m*}} \left[\Ab \eeb{\omegab^\rho}\right]: \eeb{\vb}\,\dV - \int_{Y_m} \left[\bar\gb^T\cdot\nabla \eta^\rho\right]: \eeb{\vb} \,\dV &= 0 , \\ \int_{Y_m} \left[\bar\gb:\eeb{\omegab^\rho} + \bar\db \cdot\nabla \eta^\rho\right]\cdot\nabla\psi \,\dV &= -{1\over \vert Y\vert}\int_{\Gamma_{m*}} \psi\, \dS.$

Find $\hat\omegab^k$ , $\hat\eta^k$ such that for all $\vb$ , $\psi$ and for any $k = 1, 2, \dots, k^c$ ( $k^c$ is the number of conductors)

(8)¶ $\int_{Y_{m*}} \left[\Ab \eeb{\hat\omegab^{k}}\right]: \eeb{\vb}\,\dV - \int_{Y_m} \left[\bar\gb^T\cdot\nabla \hat\eta^k\right]: \eeb{\vb} \,\dV &= 0, \\ \int_{Y_m} \left[\bar\gb:\eeb{\hat\omegab^k} + \bar\db \cdot\nabla \hat\eta^k\right]\cdot\nabla\psi \,\dV &= 0.$

The microscopic sub-problems are solved with the periodic boundary conditions and $\hat\eta^k = \delta_{ki}$ on $\Gamma^i_{m*}$ for $i = 1, 2, \dots, k^c$ , $\Gamma^i_{m*} = \overline{Y_m} \cap \overline{Y_c^i}$ is the interface between the matrix part $Y_m$ and $i$ -th conductor $Y_c^i$ . Functions $\eta^{ij}$ , $\eta^P$ and $\eta^\rho$ are equal to zero on $\Gamma_{m*}$ .

With the characteristic responses obtained by solving local sub-problems, the homogenized material coefficients $\Ab^H$ , $\Bb^H$ , $\Hb^{H,k}$ , $M^H$ and $Z^{H,k}$ can be evaluated using the following expressions:

(9)¶ $A_{ijkl}^H & = {1\over \vert Y\vert} \left[\int_{Y_{m*}} \left[\Ab \eeb{\omegab^{ij} + \Pib^{ij}}\right]: \eeb{\omegab^{kl} + \Pib^{kl}}\,\dV + \int_{Y_m} \bar\db \nabla\eta^{ij} \cdot\nabla\eta^{kl}\,\dV\right],\\ B_{ijkl}^H & = {1\over \vert Y\vert} \left[\int_{Y_{m*}} \left[\Ab \eeb{\omega^{P}}\right]: \eeb{\Pib^{kl}}\,\dV - \int_{Y_m} \bar\gb \nabla \Pib^{ij} \cdot\nabla\eta^{P}\,\dV\right] + \Phi \delta_{ij},\\ H^{H,k}_{ij} & = {1\over \vert Y\vert} \left[\int_{Y_{m*}} \left[\Ab \eeb{\hat\omegab^k}\right]: \eeb{\Pib^{ij}}\,\dV - \int_{Y_m} \left[\bar\gb:\eeb{\Pib^{ij}}\right]\cdot\nabla\hat\eta^k\,\dV\right],\\ M^H & = {1\over \vert Y\vert} \left[\int_{Y_{m*}} \left[\Ab \eeb{\omega^{P}}\right]: \eeb{\omega^{P}}\,\dV + \int_{Y_m} \bar\db \nabla\eta^{P} \cdot\nabla\eta^{P}\,\dV\right] + \Phi \delta_{ij},\\ Z^{H,k} &= -{1\over \vert Y\vert}\int_{\Gamma_c} \hat\omega^k \cdot \nb\, \dS$

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} [\Ab^H \eeb{\ub^0} - p^0 \Bb^H]: \eeb{\vb}\,\dV &= - \int_{\Omega} \eeb{\vb} : \sum\limits_k \Hb^{H,k} \bar\vphi^k \,\dV + \int_{\Omega} f \cdot \vb \dV,\\ \int_{\Omega} q^0(\Bb^H:\eeb{\ub^0} + p^0 M^H) \dV &= \int_{\Omega} q^0 \sum\limits_k Z^{H,k} \bar\vphi^k \,\dV.$

We assumed that the volume electric charge is equal to zero, otherwise we would need extra coefficients and right hand side terms.

Numerical simulation¶

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

./simple.py example_poropiezo-1/poropiezo_macro_dfc.py

This invoke the simply.py script which calculates the macroscopic problem (10) and calls the homogenization engine that solves the local subproblems (5)–(6), evaluates the homogenized coefficients (9) and finally performs the reconstruction of the solutions at the microscopic level. See [CimrmanLukesRohan2019] for more details related to the SfePy homogenization engine.

The macroscopic sample is fixed on its left side, so that no displacements are allowed, see Fig. 2 left. The defromation is induced due to piezoelectric effect, as the responce to the prescribed electric potentials $\bar\vphi^1$ , $\bar\vphi^2$ , see (10).

Fig. 2 Left - boundary conditions applied at the macroscopic level; right - revocered part of the macroscopic sample¶

The deformed shape of the sample, pressure field $p$ and the magnitude of macroscopic strain $\eb(\ub^0)$ are depicted in Fig. 3. The reconstructed strain and electric fields for a given $\veps^0$ in the part of the macroscopic domain are shown in Fig. 4.

Fig. 3 Deformed macroscopic sample and the resulting fields: left - pressure $p$ ; right - magnitude of macroscopic strain $\eb$ ¶

Fig. 4 Magnitudes of reconstructed fields: left - strain field $\eb^{mic}$ ; right - electric field $\vec{E}^{mic}$ ¶

References¶

RohanLukes2018: Rohan E., Lukeš V. Homogenization of the fluid-saturated piezoelectric porous media. International Journal of Solids and Structures, 147:110-120, 2018, DOI:10.1016/j.ijsolstr.2018.05.017
CimrmanLukesRohan2019: Cimrman R., Lukeš 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