Mesh generation for periodic 3D microstructure models and computation of effective properties

Abstract

Understanding and optimizing effective properties of porous functional materials, such as permeability or conductivity, is one of the main goals of materials science research with numerous applications. For this purpose, understanding the underlying 3D microstructure is crucial since it is well known that the materials' morphology has a significant impact on their effective properties. Because tomographic imaging is expensive in time and costs, stochastic microstructure modeling is a valuable tool for virtual materials testing, where a large number of realistic 3D microstructures can be generated and used as geometry input for spatially-resolved numerical simulations. Since the vast majority of numerical simulations is based on solving differential equations, it is essential to have fast and robust methods for generating high-quality volume meshes for the geometrically complex microstructure domains. The present paper introduces a novel method for generating volume-meshes with periodic boundary conditions based on an analytical representation of the 3D microstructure using spherical harmonics. Due to its generality, the present method is applicable to many scientific areas. In particular, we present some numerical examples with applications to battery research by making use of an already existing stochastic 3D microstructure model that has been calibrated to eight differently compacted cathodes.

Introduction

Porous media can be found in many natural as well as artificial physical, biological and chemical systems. From the composition of soils [1], [2], through which liquids seep into the ground water, to the mechanical stiffness of cements [3], [4], from battery electrodes [5], [6], [7], in which lithium ions are stored, to sponge-based filtration materials [8]: the porous microstructure of the respective system has a crucial impact on the overall behavior [9]. For example, the morphology of electrodes in lithium-ion batteries significantly influences the electrochemical properties [10], [11], [12], [13], which is the main reason why tailored structuring and manufacturing of anodes and cathodes is one promising approach to improve the performance of the cell [14], [15], [16]. Thus, it is a major issue in many research areas to design the microstructure in such a way that the overall performance, e.g. permeability, electrical conductivity, mechanical stiffness, energy density and further quantities, is optimized.

From a mathematical point of view, the impact of the 3D morphology of porous media on their macroscopic behavior, e.g., the flow rate of water through soil or the flux of lithium ions through a battery electrode, can be studied with homogenization techniques. A prominent and mathematically sound tool is periodic homogenization theory [17], which assumes that the porous medium, given as a certain domain Ω, is a periodic repetition of some representative volume element ω, see Fig. 1. This method allows to derive a set of partial differential equations (PDEs) for which the porous microstructure is not spatially resolved anymore. This significantly reduces the numerical complexity of the problem. The method is based on an asymptotic expansion of the balance equation in terms of ε, which is the ratio between a macro-scale length L and the cell-scale length ℓ, i.e., $\epsilon =\frac{\ell }{L}$. In the asymptotic limit, where $\epsilon \to 0$, a set of homogenized balance equations is then obtained, together with some porous media parameters.

Consider the decomposition $\mathrm{\Omega }={\mathrm{\Omega }}_{\mathtt{E}}\cup {\mathrm{\Omega }}_{\mathtt{S}}$, where the set ${\mathrm{\Omega }}_{\mathtt{E}}$ is simply connected and corresponds exemplarily to an electrolyte phase, and ${\mathrm{\Omega }}_{\mathtt{S}}$ is multiply connected, denoting exemplarily a solid phase. The interface between ${\mathrm{\Omega }}_{\mathtt{E}}$ and ${\mathrm{\Omega }}_{\mathtt{S}}$ is denoted by ${\mathrm{\Sigma }}_{\mathtt{E},\mathtt{S}}$. As already mentioned above, Ω is a periodic repetition of the unit cell $\omega ={\omega }_{\mathtt{E}}\cup {\omega }_S$, and the common interface ${\sigma }_{\mathtt{E},\mathtt{S}}={\omega }_{\mathtt{E}}\cap {\omega }_S$.

For a scalar balance equation with surface reactions, we have the general model

$$\mathbf{PDE1}:\{\begin{array}{cc}\frac{\partial u}{\partial t}={\text{div}}_{\mathbf{x}}\left({\mathbf{j}}_u\right)\hfill & \text{for all }\mathbf{x}\in {\mathrm{\Omega }}_{\mathtt{E}},\hfill \\ {\mathbf{j}}_u\cdot \mathbf{n}=\epsilon \cdot \underset{s}{r}{}_u\hfill & \text{on }{\mathrm{\Sigma }}_{\mathtt{E},\mathtt{S}}.\hfill \end{array}$$

This problem is formally solved via the introduction of a multi-scale expansion $u(\mathbf{x},t)=u^0(\mathbf{x},\mathbf{y},t)+\epsilon \cdot u^1(\mathbf{x},\mathbf{y},t)+\mathcal{O}({\epsilon }^2)$, where $\mathbf{y}=\frac{\mathbf{x}}{\epsilon }$, which yields a sequence of PDEs to determine the unknown functions $u^j$ in the orders ${\epsilon }^j$ of the scaling parameter ε. Briefly summarized, with periodic homogenization [17], [18], [19] one obtains the following statements:

Order ${\epsilon }^0$  yields essentially $u^0(\mathbf{x},\mathbf{y},t)=u^0(\mathbf{x},t)$, i.e., the leading order function $u^0$ is independent of the micro-scale y.

Order ${\epsilon }^1$  yields $u^1(\mathbf{x},\mathbf{y},t)={({\chi }_{\mathtt{E}}^1,{\chi }_{\mathtt{E}}^2,{\chi }_{\mathtt{E}}^3)}^T\cdot {\mathrm{\nabla }}_{\mathbf{x}}{u}^0$ as well as the condition

$$\mathbf{CP1}:\{\begin{array}{cc}{\text{div}}_{\mathbf{y}}{\mathrm{\nabla }}_{\mathbf{y}}{\chi }_{\mathtt{E}}^k\hfill & =0\text{for all }\mathbf{y}\in {\omega }_{\mathtt{E}},\hfill \\ \mathrm{\nabla }{\chi }_{\mathtt{E}}^k\cdot \mathbf{n}\hfill & =n_k\text{on }{\sigma }_{\mathtt{E},\mathtt{S}},\hfill \\ {\chi }_{\mathtt{E}}^k\hfill & \text{ periodic},\hfill \end{array}$$

for the (geometrical) corrector function ${\overrightarrow{\chi }}_{\mathtt{E}}$. The typeface ${\overrightarrow{\chi }}_{\mathtt{E}}={\overrightarrow{\chi }}_{\mathtt{E}}(\mathbf{y})$ emphasizes that ${\overrightarrow{\chi }}_{\mathtt{E}}$ is a (numerical) solution of the cell problem CP 1, which depends thus only on the micro-scale y. Since $u^0$ and thus also ${\mathrm{\nabla }}_{\mathbf{x}}{u}^0$ depend only on the macro-scale x, i.e. ${\mathrm{\nabla }}_{\mathbf{x}}{u}^0:=\overrightarrow{h}(\mathbf{x})$, the result $u^1={\overrightarrow{\chi }}_{\mathtt{E}}(\mathbf{y})\cdot \overrightarrow{h}(\mathbf{x})$ for the first successive term $u^1$ is a separation between microscopic geometrical effects and macroscopic gradients of the leading order term $u^0(\mathbf{x})$. This is a central feature of (periodic) homogenization theory and we show in Section 4 how ${\overrightarrow{\chi }}_{\mathtt{E}}$ is related to the tortuosity of a microstructure.

Order ${\epsilon }^2$  yields the PDE

$${\psi }_{\mathtt{E}}\frac{\partial u^0}{\partial t}={\text{div}}_{\mathbf{x}}({\psi }_{\mathtt{E}}{\mathit{\pi }}_{\mathtt{E}}\cdot {\mathbf{j}}_u^0)+a_{\mathtt{E},\mathtt{S}}\underset{s}{r}{}_u\text{for all }\mathbf{x}\in \mathrm{\Omega }$$

for the leading order term $u^0$ and the leading order flux ${\mathbf{j}}_u^0$, where the porous media parameters are given by

• the porosity (or phase fraction) of ${\mathrm{\Omega }}_{\mathtt{E}}$,

  $${\psi }_{\mathtt{E}}=\frac{1}{\text{vol}(\omega )}\underset{{\omega }_{\mathtt{E}}}{\int }1dV,$$

• the interfacial area of ${\mathrm{\Sigma }}_{\mathtt{E},\mathtt{S}}$,

  $$a_{\mathtt{E},\mathtt{S}}=\frac{1}{\text{vol}(\omega )}\underset{{\sigma }_{\mathtt{E},\mathtt{S}}}{\int }1dA,$$

• and the (flux) corrector,

  $${\mathit{\pi }}_{\mathtt{E}}=(\mathbf{1}-\frac{1}{\text{vol}({\omega }_{\mathtt{E}})}\underset{{\omega }_{\mathtt{E}}}{\int }\mathbf{\nabla }\left(\begin{array}{c}\hfill {\chi }_{\mathtt{E}}^1\hfill \\ \hfill {\chi }_{\mathtt{E}}^2\hfill \\ \hfill {\chi }_{\mathtt{E}}^3\hfill \end{array}\right)dV).$$

After the homogenization procedure the index ⁰ of the leading order term is typically dropped and considered as the macroscale variable. If ${\mathbf{j}}_u$ is a diffusion or heat flux, e.g., ${\mathbf{j}}_u=D_u\cdot \mathrm{\nabla }u$, the corrector ${\mathit{\pi }}_{\mathtt{E}}$ yields the effective diffusion coefficient (or conductivity) $D_{\mathtt{E}}^{\text{eff}}={\mathit{\pi }}_{\mathtt{E}}\cdot D_u$. The corrector ${\mathit{\pi }}_{\mathtt{E}}$ is thus also related to the tortuosity of the porous medium.

For the Stokes problem in a similar manner we obtain

$$\mathbf{PDE2}:\{\begin{array}{cc}\mathrm{\nabla }p-{\epsilon }^2\mu \text{div}\mathrm{\nabla }\mathbf{v}=\mathbf{f}\hfill & \text{for all }\mathbf{x}\in {\mathrm{\Omega }}_{\mathtt{E}},\hfill \\ \text{div}\mathbf{v}=0\hfill & \text{for all }\mathbf{x}\in {\mathrm{\Omega }}_{\mathtt{E}},\hfill \\ \mathbf{v}=0\hfill & \text{on }{\mathrm{\Sigma }}_{\mathtt{E},\mathtt{S}},\hfill \end{array}$$

where periodic homogenization leads to the Darcy flow

$$\mathbf{v}=\frac{1}{\mu }{\mathit{\kappa }}_{\mathtt{E}}(\mathbf{f}-\mathrm{\nabla }p)\text{for all }\mathbf{x}\in {\mathrm{\Omega }}_{\mathtt{E}}^{\text{Hom}},\text{div}\mathbf{v}=0\text{for all }\mathbf{x}\in {\mathrm{\Omega }}_{\mathtt{E}},\mathbf{v}\cdot \mathbf{n}{|}_{\partial {\mathrm{\Omega }}_{\mathtt{E}}^{\text{Hom}}}=0.$$

The corrector ${\mathit{\kappa }}_{\mathtt{E}}$ is frequently called a permeability tensor, where

$${({\mathit{\kappa }}_{\mathtt{E}})}_{j,k}=\frac{1}{\text{vol}({\omega }_{\mathtt{E}})}\underset{{\omega }_{\mathtt{E}}}{\int }\mathrm{\nabla }{\mathbf{w}}_j\cdot \mathrm{\nabla }{\mathbf{w}}_{k}dV,$$

and determined from the cell problem [20]

$$\mathbf{CP2}:\{\begin{array}{c}{\mathrm{\nabla }}_{\mathbf{y}}{q}_k-{\text{div}}_{\mathbf{y}}{\mathrm{\nabla }}_{\mathbf{y}}{\mathbf{w}}_k={\mathbf{e}}_k\text{for all }\mathbf{y}\in {\omega }_{\mathtt{E}},\hfill \\ \text{div}{\mathbf{w}}_k=0,\hfill \\ {\mathbf{w}}_k=0\text{on }{\sigma }_{\mathtt{E},\mathtt{S}},\hfill \\ q_k,{\mathbf{w}}_k\text{periodic, }k=1,2,3\text{.}\hfill \end{array}$$

For every PDE problem, e.g., PDE1 or PDE2 described above, periodic homogenization leads to a different cell problem, i.e., CP1 or CP2, which has to be solved in order to determine the effective porous media parameters. However, all of these cell problems do have in common that some stationary PDE system has to be solved on the periodic representative volume element ω. Since this is analytically possible only for a very tiny amount of geometries, the cell problems have in general to be solved numerically. And, in order to so, adequate discretizations of ω are required.

Various approaches for the discretization of the representative volume element ω exist and we briefly review exclusively those which ensure the periodicity of ω. The most simple approach is a voxel based discretization of ω, with equal edge length of the voxel in all three dimensions. This format of discretization is widely used in 3D imaging. However for the purpose of numerical calculations voxel based meshes are inappropriate, (i) because the computational degrees of freedom scale with $\mathcal{O}{(N)}^3$ for N equally sized voxels [3], or (ii) because local refinements (or coarsening), leading to $\mathcal{O}{(N)}^{3-\alpha },\alpha >0$, lead to hanging nodes,¹ which are numerically very problematic. To get rid of these problems a mesh for numerical calculations in 3D is typically built by tetrahedra, with which a 3D geometry can be discretized and locally refined without producing hanging nodes. This yields far more efficient numerical calculations and the corresponding meshes are said to be of high quality. Numerical simulations of 3D microstructures allow then for virtual materials testing [21], [22].

The generation of volumetric meshes for realistic microstructures is an interdisciplinary topic and various approaches are found throughout the literature [23], [24], [25]. A central aspect is the assumption regarding the geometrical shape of the particles or inclusions in the microstructure. For example a rather general approach of modeling inclusions in matrix materials assumes ellipsoidal shapes [26], which is explicitly exploited in the periodic mesh generation. However, tomographic imaging methods have shown that realistic 3D microstructures of various functional materials are significantly more complex. A flexible tool to model and simulate the morphology of such particle systems is hence desirable and we propose a complete pipeline for this issue.

In this paper, we propose a robust mesh generation for periodic representative volume elements of realistic microstructures, see Fig. 2, with geometrically more flexible star-shaped particles.

The method is based on a description of the microstructure in terms of spherical harmonics, a subsequent surface mesh generation of $\partial {\omega }_{\mathtt{E}}$ and $\partial {\omega }_{\mathtt{S}}$, and finally a volume mesh generation based on TetGen [27]. The proposed method can be applied to a broad spectrum of scenarios arising in different fields of research since numerous scientific problems involve solving a system of differential equations on periodic porous media. Another advantage of the presented approach is that periodic boundary conditions can be easily applied in x-, y- and z-direction as well as to an arbitrary subset of directions. This can be used for example in battery research, where the size of electrodes is typically several orders of magnitudes larger in in-plane direction compared to the thickness of the electrode such that it is reasonable to consider periodic boundary conditions in two directions.

The rest of this paper is organized as follows. In Section 2, we describe the generation of periodic 3D microstructures based on spherical harmonics and a stochastic microstructure model. Then, in Section 3, the generation of a quality volume mesh on the basis of the representation of the particle system via spherical harmonics is explained. In Section 4, some numerical examples are presented. Finally, in Section 5, the paper is concluded by a summary of the main results and an outlook to possible future research is given.