Topology optimization of unsaturated flows in multi-material porous media: application to a simple diaper model

TL;DR

This paper develops a topology optimization method for maximizing liquid absorption in multi-material porous media, demonstrated on a simplified diaper model with significant improvements over the reference design.

Contribution

It introduces a density-based multi-material topology optimization approach for unsaturated flow in porous media, applicable to complex, nonlinear transient problems.

Findings

Achieved up to 45% increase in absorbed liquid compared to reference design.

Developed sensitivity analysis enabling first-order optimization for nonlinear problems.

Validated the approach on two- and three-material configurations.

Abstract

We present a mathematical approach to optimize the material distribution for fluid transport in unsaturated porous media. Our benchmark problem is a simplified diaper model as an exemplary liquid absorber. Our model has up to three materials with vastly different properties, which in the reference configuration are arranged in parallel layers. Here, swelling is neglected and the geometry of a swollen diaper is assumed and treated as a porous medium of high porosity. The imbibition process is then simulated by solving an unsaturated flow problem based on Richards' equation. Our aim is to maximize the amount of absorbed liquid by redistributing the materials. To this end, a density based multi-material topology optimization (based on the well known SIMP model) is performed. A sensitivity analysis for the nonlinear transient problem is provided, which enables the application of first order optimization solvers. We perform two- and three-material optimization and discuss several variants of the problem setting. We present designs with up to 45% more absorbed liquid compared to the reference configuration.

Tables

Table 1. Table 1: Parameters defining the properties of the diaper materials A 𝐴 A , B 𝐵 B and C 𝐶 C (see ( 2 )–( 5 ) ).

	$A$	$B$	$C$	Description
$\varphi$	0.9	0.9	0.95	porosity
$A_h$ in 1/cm	10	5	3	coefficient in (2)
n	4	2.5	1.7	coefficient in (2)
m	1.0	1.5	1.7	coefficient in (2)
K0 in cm/s	2.0	0.5	0.15	maximum conductivity

Equations

$$\varphi \frac{\partial \theta }{\partial t}=\nabla \cdot \left[K\nabla (h+z)\right],$$

$$\frac{\theta -{\theta }_r}{1-{\theta }_r}=\frac{1}{[1+(A_{h}h{)}^n{]}^m}.$$

$$K=\frac{k\rho g}{\mu },$$

$$K(\theta )=K_0{\theta }^4.$$

$$C(h)=\frac{\text{d}\theta }{\text{d}h},$$

$$\varphi C(h)\frac{\partial h}{\partial t}=\nabla \cdot \left[K\nabla (h+z)\right],$$

$$\varphi \frac{{\theta }^{i+1}-{\theta }^i}{\Delta t}=\nabla \cdot \left[K^{i+1}\left(\nabla h^{i+1}+z\right)\right]$$

$${\theta }^{i+1,j+1}={\theta }^{i+1,j}+C^{i+1,j}\left[h^{i+1,j+1}-h^{i+1,j}\right]+O(h.o.t),$$

$$\varphi C^{i+1,j}\frac{h^{i+1,j+1}-h^{i+1,j}}{\text{Δt}}=\nabla \cdot \left[K^{i+1,j}\nabla (h^{i+1,j+1}+z)\right]-\varphi \frac{{\theta }^{i+1,j}-{\theta }^i}{\Delta t}.$$

$${\mathbf{w}}_j=(w_j^A,w_j^B,w_j^C)\in \{0,1{\}}^3,$$

$$||{\mathbf{w}}_j||=w_j^A+w_j^B+w_j^C=1,$$

$${\theta }_j({\mathbf{w}}_j,h)$$

$$C_j({\mathbf{w}}_j,h)$$

$$K_j({\mathbf{w}}_j,\theta )$$

$${\varphi }_j({\mathbf{w}}_j)$$

$${\mathbf{\rho }}_e\in [0,1]\times [0,1]\times [0,1]=[0,1{]}^3\subset {\mathbb{R}}^3$$

$${\rho }_e^A+{\rho }_e^B+{\rho }_e^C=1.$$

$${\rho }_e\in [0,1]$$

$$\varphi (\rho )C(\rho ,h)\frac{\partial h}{\partial t}-\nabla \cdot K(\rho ,h)\nabla h=0,$$

$$\mathbf{M}(\mathbf{\rho },\mathbf{h})\frac{\partial \mathbf{h}}{\partial t}-\mathbf{K}(\mathbf{\rho },\mathbf{h})\mathbf{h}=0.$$

$$\displaystyle\min_{\bm{\rho}}J({\bm{\rho}},{\bm{h}})\Big{|}_{t=T}$$

$$\text{s.t. }\mathbf{M}(\mathbf{\rho },\mathbf{h})\frac{\partial \mathbf{h}}{\partial t}-\mathbf{K}(\mathbf{\rho },\mathbf{h})\mathbf{h}$$

$$g(\mathbf{\rho })$$

$$\widetilde{\bm{M}}\big{|}_{t}\,\frac{\partial{{\bm{h}}(t)}}{\partial t}-\widetilde{\bm{K}}\big{|}_{t}\,{\bm{h}}(t)={\bm{0}}\quad\forall t\in[0,T],$$

$$\Phi=\Phi({\bm{\rho}},{\bm{h}})\Big{|}_{t=T}=J({\bm{\rho}},{\bm{h}})\Big{|}_{t=T}+\int_{0}^{T}{\bm{\lambda}}^{\top}\left(\widetilde{\bm{M}}\,\frac{\partial{{\bm{h}}}}{\partial t}-\widetilde{\bm{K}}\,{\bm{h}}\right)\,\mathrm{d}\,t,$$

$$\frac{\text{d}\Phi }{\text{d}{\rho }_e}$$

$$\tilde{\mathbf{M}}\frac{\partial \mathbf{\lambda }}{\partial t}+\tilde{\mathbf{K}}\mathbf{\lambda }=0$$

$$\widetilde{\bm{M}}\,{\bm{\lambda}}\bigg{|}_{t=T}=-\left(\frac{\partial{J}}{\partial{h}}\right)^{\top}\bigg{|}_{t=T}.$$

$$\frac{\text{d}\,{\Phi}}{\text{d}\,{\rho_{e}}}=\frac{\partial{J}}{\partial\rho_{e}}\Big{|}_{t=T}+\int_{0}^{T}\left({\bm{\lambda}}^{\top}\frac{\partial{\widetilde{\bm{M}}}}{\partial\rho_{e}}\,{{\bm{h}}}-{\bm{\lambda}}^{\top}\frac{\partial{\widetilde{\bm{K}}}}{\partial\rho_{e}}\,{\bm{h}}\right)\,\mathrm{d}\,t,$$

$$-\widetilde{\bm{M}}\frac{\partial{\overline{{\bm{\lambda}}}}}{\partial t}+\widetilde{\bm{K}}\,\overline{{\bm{\lambda}}}={\bm{0}}\;\Bigg{|}_{T-\tau},\quad\overline{{\bm{\lambda}}}(\tau)={\bm{\lambda}}(T-t),$$

$$\widetilde{\bm{M}}\,\overline{{\bm{\lambda}}}\bigg{|}_{\tau=0}=-\left(\frac{\partial{J}}{\partial{h}}\right)^{\top}\bigg{|}_{t=T}.$$

Full text

††volume: 00

\journalname

Communications in Nonlinear Science and Numerical Simulation \runauth \jidcnsns

\dochead

1 Introduction

Superabsorbent polymers (SAPs) are granular materials that revolutionized the market of disposable hygiene products such as diapers since their commercialization in the late seventies. The major feature of SAPs is that they can absorb and retain over thirty times their own weight of liquid while swelling. Because of this property, SAPs are widely used in sanitary products as well as in the agricultural industry. One particularly important product are modern diapers, which are composed of high-tech multi-material SAPs designed to maximize liquid absorption, while remaining dry at the contact surface with the skin. A diaper is typically composed of three parallel layers of materials with distinct properties [1]. The first layer is designed to evacuate the fluid of the surface as fast as possible, whereas the second one distributes the fluid in the longitudinal direction and the third is the actual SAP that stores it.

Building upon theoretical principles derived for deformable porous media in soil mechanics [2], Diersch et al. [1] derived a set of equations describing the flow of liquid and the swelling in absorbent hygiene products. The first equation is a generalized Richards’ equation extended to absorbing and swelling materials. Richards’ equation [3] describes the motion of liquid in an unsaturated media and finds widespread use in soil science to estimate water infiltration depth. It reduces a complex two-phase fluid problem to a partial differential equation for a single variable (the capillary pressure) and is the analogous of Darcy’s equation, but for unsaturated media. Note that in Darcy’s equation pressure is the only variable, whereas in Richards’ equation there is also the moisture content (saturation). For a given material there is a pressure-moisture curve, which must be measured experimentally and is typically fitted with highly nonlinear empirical laws [4]. However, a grain-scale modelling technique was recently developed to predict capillary pressure-moisture curves for SAP as a function of their granular structure [5]. Similarly, the permeability of the material depends also nonlinearly on the saturation. The second and third equation derived by Diersch et al. [1] are the conservation of mass in the solid phase and a relationship for the solid strain, respectively, supplemented additionally by several complex constitutive equations. While their equations model the physics of fluid absorption in SAP faithfully, simulating them accurately is a challenge. In particular, the large swelling deformations requires moving meshes with automated refinement and several adaptive techniques must be combined to obtain robust and accurate results [6].

An interesting question that can be numerically addressed is whether the performance of a diaper can be optimized by changing the spatial distribution of the three materials. In fact, it is common for density based topology optimization to introduce a porous model as differentiable design variable (pseudo density). Typically problems are solved where the final design consist of a topology of fully permeable and impermeable material. For example, Guest and Prévost [7] coupled a fluid problem and topology optimization with a Darcy-Stokes approach and found out the optimal structure of porous media based on minimizing dissipated power; Shou and Fan [8] investigated the fast and controlled capillary flows in homogeneous porous structure and optimized the porous structure for imbibition of capillary flow with viscous resistance; Hübner et al [9] employed the Biot model for shape optimization of microstructure saturated viscous fluids and reported a two-scale structure of porous material based on effective poroelastic coefficients. However, to the best of our knowledge, multi-layer optimization of transient unsaturated flows in porous media has not been done before and is the main focus of our work. For this purpose, we consider a highly simplified model of a diaper. The main simplification is that the swelling of the SAPs, as liquid is absorbed, is not accounted for. Instead we assume that the diaper is in its swollen state before the absorption of liquid starts and the materials are treated as porous media. The displacement of air as the liquid fills the diaper is modeled with Richards’ equation. This simplification is necessary because the complexity of the equations with swelling, and more importantly the numerical techniques necessary to solve them accurately, would render the optimization problem intractable.

In this paper, we describe a numerical method to compute the flow of liquid in an idealized diaper. The performance of this diaper is subsequently improved by optimizing the spatial distribution of the materials. The target is to maximize the volume of absorbed liquid. The assignment of material is at the level of the mesh discretization, permitting nearly arbitrarily complex designs. This approach, called topology optimization, belongs together with shape optimization to the field of structural optimization. The underlying optimization model is the so-called solid isotropic material with penalization (SIMP), dating back to Bendsøe [10]. Comparisons of its performance with alternative structural optimization approaches can be found in the literature [11, 12].

This paper is structured as follows. The model and the numerical method for the solution of the diaper problem are presented in Sec. 2 and Sec. 3, respectively. In Sec. 4 we introduce the topology optimization approach, including the formulation of the optimization problem. Our optimization results are presented in Sec. 5 and Sec. 6, with variants in Sec. 6.1. The conclusions are given in Sec. 7.

2 Simplified numerical model of a diaper

We model the diaper as a Cartesian box $\Omega$=$[0,L_{x}]\times[0,L_{y}]\times[-L_{z},0]$ and gravity acts in the negative $z$ direction. The length, width and depth of the box are $L_{x}$=25 cm, $L_{y}$=10 cm and $L_{z}$=2 cm, respectively, and the total diaper volume $v=500\,\text{cm}^{3}$. All materials are considered as non-swelling porous media, so our model represents a diaper swollen to its maximum volume even if no water is present. With this simplification a dry diaper can be seen as a porous medium of very high porosity, which models SAPs having the capacity to absorb a high volume of liquid. As the liquid is discharged, the pores get progressively filled.

A schematic of the computational model of the diaper is shown in Fig. 1. Liquid enters the box from the top at a 5 cm $\times$ 5 cm rectangle during a defined discharge time followed by an equilibration phase. The goal of the optimization problem is to propose a diaper design that can absorb more liquid than the reference three-layer design. The new design should use the same volume of material $A$, $B$ and $C$ but achieve a higher efficiency by optimizing their spatial distribution. The precise mathematical formulation of the optimization problem is given in Sec. 4.

2.1 Governing equations and material parameters

Richards’ equation [3] describes the motion of liquids in unsaturated porous media. In our formulation, the hydraulic pressure head $h=p/(\rho\,g)$ (in cm) is used as the main variable, where $p$ is the capillary pressure, $\rho$=1g/cm3 the liquid density (water at 20∘C) and $g$=981cm/s2 the gravitational acceleration. A full liquid saturation corresponds to zero capillary pressure, $h=0$, whereas completely dry material has $h\rightarrow-\infty$. Richards’ equation for the hydraulic pressure head reads

[TABLE]

where $\phi$ is the porosity of the material and $K$ the hydraulic conductivity. For a solid material (without pores) $\phi$=0, whereas $\phi=1$ corresponds to the limit in which there is no solid matrix. The dimensionless moisture content $\theta(h)$ can in principle take values between [math] and $1$, corresponding to dry and completely filled pores respectively. Here the van Genuchten parametrization is used to model the behavior of the SAPs

[TABLE]

In practice a small residual moisture content $\theta_{r}$=0.0025 is used in order to avoid that the pressure head diverges ($h\rightarrow-\infty$) in dry material. The values of the material parameters $A_{h}$, $n$ and $m$ have been chosen to qualitatively model the absorption behavior of diapers, while keeping the problem tractable numerically (see Tab. 1 for their specification). The values used for material C are similar to [6], who compared their results to laboratory experiments with diapers. The left panel of Fig. 2 shows the moisture content as a function of the pressure head for the three materials. Note how the uppermost material $A$ remains essentially dry up to $h$=-1 cm, whereas bottommost material $C$ starts to get significantly wet already at $h$=-10 cm. The buffer material $B$ features intermediate properties.

The materials are assumed isotropic. The hydraulic conductivity is a scalar

[TABLE]

where $k$ is the permeability of the material and $\mu$=0.01g/(cm s) the dynamic viscosity of water at 20∘C. As materials are progressively wetted, their hydraulic conductivity increases, and this is modeled by the expression

[TABLE]

The values of the maximum conductivity $K_{0}$ for each material can be found in Tab. 1. It is worth noting that the maximum conductivity of material $A$ is significantly larger than that of $C$, so that liquid is rapidly evacuated towards the bottom of the diaper. The mid panel of Fig. 2 shows the conductivity as a function of pressure head for the three materials.

Another key material property is the moisture capacity

[TABLE]

which quantifies how much liquid can be stored at a given pressure head. The three materials feature peak moisture capacity at different pressures. Material $A$ absorbs liquid efficiently only when it is quite wet, i.e. $h$=-0.1 cm, whereas material $C$ has a quite flat profile, allowing liquid to be efficiently absorbed for a wide range of pressure heads. Material $B$ features an intermediate behavior.

We note that equation (1) is known as the mixed form of Richards’ equation. By using the relationship between moisture content and pressured head $\theta(h)$, one obtains the head form of Richards’ equation

[TABLE]

which has the advantage that $h$ is the only unknown. Hence Richards’ equation is essentially a transient Darcy equation for the pressure head $h$, where the effective porosity $\phi C(h)$ and the hydraulic conductivity $K(h)$ depend nonlinearly on $h$. Because a direct discretization of the head form of the equation results in large mass losses [13], Richards’ equation was here solved in its mixed form (1) as detailed in Sec. 3.

2.2 Initial and boundary conditions

At time $t$=0 the diaper is dry and the liquid discharge phase begins. When liquid is deposited on top of a porous medium, such as the diaper considered here, it forms a puddle, which is then subsequently absorbed. The absorption speed and hence the flux at the surface are unknown a priory. In fact, the flux depends nonlinearly on the permeability and capacity, which depend on the moisture content itself. In a nutshell, the flux is a function of the material properties and of the liquid distribution inside the diaper. We model the puddle as a 5 cm $\times$ 5 cm rectangle centered at $(12.5,5,0)$, where fully wet material ($\theta=1$) is imposed as boundary condition. At the rest of the domain’s boundary the flux is set to zero to enforce that no liquid exits the diaper. During a subsequent equilibration phase, zero flux boundary conditions are imposed at all boundaries. In this problem, the wetting dynamics and the content of water which can be absorbed is mainly determined by the large moisture gradients near the inlet, and the large gradients of permeability at the interface between materials. This requires very fine meshing in the axial directions and small time-step sizes during the liquid discharges.

2.3 Simulation of the three-dimensional reference design

We considered a discharge resulting in a puddle of 24 s duration, which is followed by an equilibration phase of 120 s. The left panel of Fig. 3 shows the evolution of the total liquid content in the diaper with time. At $t=0$, there is a rapid increase of the liquid content, because of the strong pressure gradients between the wet inlet and dry diaper. As the liquid progressively fills the diaper, gradients are reduced and the absorption rate decreases. During the equilibration phase liquid is distributed throughout the diaper, while the total liquid content remains constant. This indicates that our simulation conserves mass. We note that because of the dominance of capillary forces, gravity was neglected in our simulations. This corresponds to eliminating the term $z$ from the right-hand-side of Richards’ equation (1) and is justified in the next section.

2.4 Influence of gravity

In order to study the influence of gravity on the wetting dynamics, a two-dimensional three-layer simulation of the reference structure was performed. We considered the 5 cm $\times$ 2 cm symmetric slice cut across the width and depth shown in Fig. 1 and the boundary conditions as given in Fig. 4. The liquid load was applied on the left half on top of the domain, the right half of the top boundary, as well as the right and bottom sides, are set to outflow to let the liquid pass through, the left boundary is set to be symmetric in order to reduce the computing time. The simulation started with the liquid coming from the loading region and ran until liquid started to pass through the right boundary. The whole simulation time was 7.5s.

Fig. 5 shows that the moisture content of the reference three-layer design at $t=$5s is indistinguishable in the cases with and without gravity. Fig. 6 shows the development of the deepest location of penetration with respect to time. It can be observed that although the deepest location is slightly larger with the effects of gravity, the difference is still fairly small. Furthermore, the simulation with gravity is substantially more expensive because it requires much higher grid resolution, resulting in too expensive simulations for three-dimensional optimization. In what follows, gravity is neglected.

3 Numerical method

Because of the nonlinearities in Richards’ equation, care needs to be taken to enforce mass conservation in the numerical simulations. Celia et al. [13] developed a modified Picard iteration of the mixed form of the equations (1), which conserves mass in free drainage problems. Although several improvements and extensions of this scheme have been proposed in the literature [14, 15, 16, 17, 18], we here employed the popular method of Celia et al. [13] because it gives an excellent compromise between accuracy and efficiency in simple geometries, and simple implementation. The reader is referred to List et al. [19] for a recent review on methods to linearize the Richards’ equation. In what follows, we briefly summarize the modified Picard scheme used in this paper. Time was advanced with the implicit Euler scheme,

[TABLE]

which required the solution of a nonlinear problem every time step $i+1$. This was solved iteratively with the Picard method [13], whereby the moisture content at time step $i+1$ is expanded in a Taylor series with respect to $h$ as

[TABLE]

and the superscript $j+1$ denotes the current Picard iteration. Neglecting high order terms and inserting this expression into (1), the Picard iteration reads

[TABLE]

Hence at each Picard iteration a (linear) Helmholtz equation for $h^{i+1,j+1}$ must be solved. The spatial derivatives were discretized with the second-order central finite-volume method. Our code was implemented in OpenFOAM [20] and was successfully benchmarked against [13]. Note how eq. (9) is identical to the Picard equation for the head form of the Richards’ equation, but with the additional term in the right-hand-side containing the moisture content.

An important difficulty in our problem is that the diaper contains three materials with very different properties. Numerical studies of two-layer [21, 22] and multi-layer [23] domains have been reported in the soil-science literature. Note however that the material properties of diapers varies significantly and are numerically challenging to treat. Our analysis domain has a discretization of $50\times 20\times 35$ cells, where the cell height in $z$-direction grows from very thin cells at the inlet plane towards rather thick ones at the other end.

A further difficulty in our problem is that the time-step required for the Picard iterations to converge varies across several orders of magnitude. At the beginning of the liquid discharge, small time-step sizes ($\Delta\,t\approx 10^{-4}$s) were required, whereas during the equilibration periods $\Delta t\approx 1$s. Hence the time-step size was dynamically adjusted based on the number of Picard iterations, whereby 20 iterations were considered as optimal. Then depending on whether 20 Picard iterations were exceeded or not in the previous time step, $\Delta t$ was divided or multiplied by 1.1. If the number of Picard iterations was below 6, the next time step was increased by 1.5. In addition to this criterion to determine the time-step dynamically, we set a hard limit of 50 Picard iterations. If this was exceeded, $\Delta t$ was reduced by $75\%$ and the time-step was re-computed. If necessary, this was repeated until the number of iterations dropped below 50. Finally, we note that $\Delta t$ was adjusted to resolve the switching between inflow and zero flux boundary conditions at $t=24$s. This was done to avoid wasting time by reaching the 50 iteration limit repeatedly. The time-step size and number of Picard iterations required in the simulation of the reference three-layer design is shown in the right panel of Fig. 3.

3.1 Treatment of the multi-material domain

Our solver was implemented to handle three materials in one function. For this purpose, we defined a three-element vector determining the material of each computational volume $j$,

[TABLE]

where

[TABLE]

and only one component of ${\bm{w}}_{j}$ was set to one and the others to zero at each computational volume, corresponding to an acceptable (discrete) material choice. The generalized material properties can be written as follows

[TABLE]

4 Topology Optimization

We applied topology optimization to distribute the materials within the design domain $\Omega$. In the $x$- and $y$-directions, the design discretization of $\Omega$ was identical to the finite-volume discretization used for the numerical simulation. In the $z$-direction, each design cell could contain several finite volumes. In the following, we refer to $N$ as the total number of design cells (with index $e$), whereas $D$ is the total number of computational volumes (with index $j$).

4.1 Parametrization

To perform gradient based optimization, we relaxed the selection variable ${\bm{w}}_{j}$ (10) to a continuous weighting variable

[TABLE]

with $1\leq e\leq N$ and the constraint

[TABLE]

In topology optimization, the design variables ${\bm{\rho}}_{e}$ are commonly called pseudo density. We use bold symbols and lowercase letters to indicate vectors. Hence the full vector of design variables of size $3\,N$ is given in bold as is the set of variables assigned to a single element for the multimaterial case. Three scalar variables per element in (17) are written as elementwise vector in (16) which spans the three-dimensional continuous space beween 0 and 1. In standard topology optimization, the distribution of a single material is optimized, whereas we here aim to solve a multi-material optimization problem similar to the class of problems treated by Erik Lund and co-workers using the discrete material optimization approach, see e.g. [24] and earlier work. To enforce condition (17), $N$ constraints are required when we model the optimization problem. This approach can be easily extended to more than three materials.

When we restrict ourselves to only two materials, a scalar design variable $\rho_{e}$ is sufficient with

[TABLE]

and the weights of two arbitrary materials, e.g. $A$ and $C$, are $\rho_{e}^{A}=\rho_{e}$ and $\rho_{e}^{C}=1-\rho_{e}$ with $\rho_{e}^{B}=0$. This approach is known as bi-material optimization (see [25]) and the norm condition is implicitly fulfilled.

The idea of the optimization is to construct material properties (13)–(15) using ${\bm{\rho}}_{j}$ instead of ${\bm{w}}_{j}$. We note that only discrete material selections ${\bm{\rho}}_{j}\in\{0,1\}^{3}$ allow a meaningful physical interpretation. With the final designs sufficiently close to discrete material distributions the effect of mixed material on the different material properties (13)–(15) plays only a minor role and is thus not needed to be analyzed it in detail.

4.2 Sensitivity Analysis

The purpose of the sensitivity analysis is to provide the first-order derivatives of the state dependent functions with respect to the design variables, entering the optimization problem as objective function or constraints. Having the sensitivities available, a suitable first-order optimizer toolbox can be used to solve the optimization problem. In general, this is vastly more efficient than utilizing derivative-free optimizers [26] and is here of particular importance, as we solve a high-dimensional nonlinear problem (with $3\,N$ design variables for the three-material case).

4.2.1 Semi-Discrete System

We consider the head form of Richards’ equation (1) and denote explicitly the dependencies on the state variable pressure head $h$ and the material parametrization by the design vector ${\bm{\rho}}$. If gravity is neglected, this reads as

[TABLE]

and after spatial discretization with the finite-volume method, we obtain the semi-discrete system of equations

[TABLE]

The system is still continuous in time. The mass matrix ${\bm{M}}$ and stiffness matrix ${\bm{K}}$ are discretizations of their continuous counterparts given in (19), see [13].

We can now formulate a semi-discrete algebraic generic optimization problem with objective function $J$ depending on the design ${\bm{\rho}}$ and state ${\bm{h}}$, as well as a constraint function $g({\bm{\rho}})$ independent of the state. The number of constraint functions and their actual type (equality or inequality) is not of importance for the following sensitivity analysis

[TABLE]

Note that we are interested in the function value for $J$ at final time $T$ only.

4.2.2 Derivation of the Gradients

The sensitivity analysis presented in the following closely follows Dahl, Jensen and Sigmund [27], however our system (20) is nonlinear and has no second time-derivative. To find the adjoint formulation, we consider the residual of the state equation point-wise

[TABLE]

where \widetilde{\bm{M}}\big{|}_{t}={\bm{M}}({\bm{\rho}},{\bar{\bm{h}}})\big{|}_{t}, \widetilde{\bm{K}}\big{|}_{t}={\bm{K}}({\bm{\rho}},{\bar{\bm{h}}})\big{|}_{t} and ${\bar{\bm{h}}}:[0,T]\to\mathbb{R}^{D}$ is the converged state solution of (20), where $D$ is the number of computational volumes resulting from the space discretization. The tilde denotes explicit dependence on ${\bm{\rho}}$.

Introducing a Lagrange multiplier ${\bm{\lambda}}:[0,T]\to\mathbb{R}^{D}$ we obtain

[TABLE]

where the integral term in the equation is zero because of the residual (22) and the objective function value is not changed. Differentiation yields

[TABLE]

The critical term to determine is the derivative of the state ${\bm{h}}$ with respect to the design variables $\rho_{e}$. In the following, we show that this vanishes with the help of the Lagrange multiplier. Applying integration by parts yields

[TABLE]

By collecting the terms with $\frac{\partial{{\bm{h}}}}{\partial\rho_{e}}$ we obtain

[TABLE]

The unknown terms $\frac{\partial{{\bm{h}}}}{\partial\rho_{e}}$ vanish if ${\bm{\lambda}}$ solves the terminal value adjoint equation

[TABLE]

and

[TABLE]

Here we make use of the fact that $\widetilde{\bm{M}}$ and $\widetilde{\bm{K}}$ are symmetric matrices. As there is no dependency of the initial state solution on the initial material distribution ([27]), also the term for $t=0$ vanishes.

The terminal value adjoint problem is transformed into an initial value problem by variable substitution. Hence the gradient is calculated as

[TABLE]

where ${\bm{\lambda}}$ solves the initial value adjoint problem

[TABLE]

with initial condition

[TABLE]

Note the change of sign for the time derivative in (24). In the adjoint problem the mass and stiffness matrices do not depend on $\overline{{\bm{\lambda}}}$, hence no Picard iterations are necessary for the numerical solution of the adjoint problem. On the other hand, we need to store for each time step the converged state ${\bm{h}}$ and store or reconstruct $\widetilde{\bm{M}}$ and $\widetilde{\bm{K}}$ as well.

Because of (25), we need to solve a separate transient adjoint problem for any state dependent function. The negative derivative of the function with respect to the state is the initial condition.

4.3 Cost and Constraint Functions

The total volume of absorbed liquid is

[TABLE]

where $\nu_{e(j)}$ is the volume of the computational cells $\Omega_{j}$ and $e:\mathbb{R}^{D}\to\mathbb{R}^{N}$ maps design cells to computational cells. Again, we skip explicit dependence on time in the notation. For the objective function, we are interested in the function value only at final time $t=T$. We recall that the reason why we distinguish between design and computational cells is that we must use a very fine non-uniform mesh discretization in $z$-direction for numerical reasons. We thus combine several computational cells to one design cell.

The amount of a given material, also known as resource constraint, is

[TABLE]

where $v_{e}$ is the volume of design cell $e$. For the bi-material case, $v({\bm{\rho}})$ simply measures the amount of material $A$. The amount of the second material is implicitly given by $1-v({\bm{\rho}})$. For the three-material case we require functions for two materials, here we choose the materials $A$ and $B$ such that the volume constraints $v^{A}({\bm{\rho}})$ and $v^{C}({\bm{\rho}})$ add to $\rho_{e}^{A}$ and $\rho_{e}^{C}$, respectively.

In many applications in topology optimization, it is necessary to regularize the design space [28]. Typically the change of the design from cell to cell is limited by filtering the design [29] or with the Heaviside projection method [30]. Here we chose to use explicit regularization by slope constraints [31]. The principal idea of slope constraints is to restrict the partial derivative of the design

[TABLE]

which are implemented with forward finite differences for every design element

[TABLE]

Choosing $c_{s}=1$ disables the slope constraints, $c_{s}=1/3$ requires at least two intermediate design cells for a full change from 0 to 1, or vice versa. Slope constraints come with a high number of linear constraints. Here we used the commercial SNOPT [32] package as optimizer to handle the high number of linear constraints. In our three-material case, it is sufficient to apply slope constraints on two materials. We want to mention that alternatively other forms of regularization, like density filtering, might be used.

Regularization still leaves room for isolated spots of intermediate material. One possible reason for this to happen might be that such material is irrelevant or negligible, e.g. in areas of the design domain too far away to be reached by moisture during the simulation. Moreover, surplus material stemming from the resource constraints might be distributed without benefit.

As a remedy against mixed materials, we used constraints of the type

[TABLE]

which are known as grayness constraints. Binary material (zero or one) has grayness zero, complete intermediate material (0.5) results in grayness one. Grayness constraints are known for the tendency to lock designs in local optima (which might be poor) as a smooth transition from void to solid or solid to void requires a temporary increase of the grayness which might violate the constraint.

A second problem is related to the concurrent application of slope constraints, which enforce intermediate densities at material interfaces. A too strict grayness constraint acts indirectly as control of the size of the interface, which might have a much more significant impact than mere removal of spots. For this reason, we applied grayness constraints only in a second stage of the optimization process starting from an already optimized design. In order to obtain a bound for the grayness constraint, we evaluated the grayness resulting from active slope constraints and used this value with a slight supplement.

4.4 Problem Formulation

Having introduced the objective and constraints, we can now formulate the optimization problem, the general work flow is shown in Fig. 7. For the bi-material case we optimize for the materials $A$ and $C$ with the design variable $\rho_{e}$ and $1-\rho_{e}$ expressing the fraction of material $A$ and $C$, respectively. Material $B$ has contribution zero. The bi-material optimization problem reads

[TABLE]

Here we effectively have a single variable $\rho_{e}$ for each of our $N$ design cells, see Sec. 4.1. Moreover, volume, grayness and slope constraints act on $\rho_{e}$ only.

In the full three-material optimization model we have $3\,N$ design variables $\rho_{1}^{A}$, $\rho_{1}^{B}$, $\rho_{1}^{C}$, $\rho_{2}^{A}$, $\ldots$, $\rho_{N}^{C}$, but it is sufficient to double constraints (35)–(40) for a second variable only, because the third one is implicitly controlled. With the additional local condition (17), we obtain the full three-material problem formulation

[TABLE]

It is noted that the third material variable could be expressed as $\rho_{e}^{B}=1-\rho_{e}^{A}-\rho_{e}^{C}$, but this would only reduce the number of variables and not the number of constraints. In topology optimization, the number of constraints is more critical in terms of numerical performance.

5 Results of the bi-material topology optimization

For simplicity, we began with a bi-material optimization problem, whose formulation and presentation of the results are much simpler than for the three-material problem. In this section, the spatial distributions of materials $A$ and $C$ were optimized, whereas no buffer material ($B$) was used in the designs. The corresponding two-layer reference design consisted here of a layer of material $A$ (25% volume fraction) directly on top of a layer of material $C$ (75% volume fraction). As for the reference three-layer design described in Sec. 2.3, the liquid discharge lasted for 24 s and this was followed by an equilibration phase of 120 s. With this setting, the reference two-layer design absorbed $64.5$ cm3, which is much more than the $41$ cm3 of the reference three-layer design, where $A$ has a volume fraction of only 5% and $B$ 10%. Clearly, increasing the volume of the most permeable material ($A$) results in a much larger volume of absorbed liquid.

5.1 Unregularized Problem

The first problem was posed without grayness and slope constraints, hence the problem consisted of equations (33), (34), (35) and (40). The optimized design is shown in Fig. 8a, where the discharge is from below to aid visualization of the structure. In this design, material $A$ is removed from the areas far away from the discharge inlet and is placed below it, forming a thick layer of complex form. Because of the high permeability of material $A$ this design allows evacuating the fluid faster from the discharge inlet, thereby increasing the flux of liquid into the diaper. This leads to a substantially better performance. Compared to the initial design, the total absorbed volume (objective function) is increased by $45\%$ ($64.5\to 93.4$cm3). Note that with a volume fraction of $25\,\%$ for material $A$, not all of it is necessary and as a result part of it remains on the top (bottom in the visualization) or is moved to other regions which remain dry during the whole simulation time. The obtained result shows some mixed material, especially close to the discharge area. For this reason, we performed a material-fraction preserving rounding, which led to the material distribution in Fig. 8b. Interestingly, this rounding has almost no impact on the objective value ($93.4\to 93.2$ cm3).

5.2 Regularized Problem

As a second problem, we added the slope constraints (37)–(39) to the bi-material problem. The slope parameter was set to $c_{s}=0.3$, hence approximately two intermediate elements were necessary to change from pure material $A$ to $C$ and vice versa. The resulting design is shown in Fig. 9a and is similar to that of Fig. 8a. However, because of the imposed slower transition between materials, there is more mixed material than enforced by the slope constraints and there is significantly less material $A$ left over on the top layer (bottom in the visualization) far away from the discharge inlet. The restriction of the design space by slope constraints leads to an objective function value of $91.1$ cm3, which is still an improvement of $41\%$ compared to the layered reference design.

Next, we calculated the grayness induced by the slope constraints for the design shown in Fig. 9a, as described in Sec. 4.3, and added the value as bound for the additional grayness constraint (36). We solved this second problem starting from the solution obtained from the slope constraint reference problem. The resulting binary design is shown in Fig. 9b. The additional grayness function has almost no impact on the obtained value of the objective function. For the binary design the volume of absorbed liquid even increases slightly to $92.5$ cm3. This is an indication that intermediate material, here enforced by the slope constraints, is obviously not necessarily advantageous for such bi-material problems.

5.3 Effect of the liquid discharge and equilibration times

Next, we studied the effect of changing the duration of the discharge and equilibration phases. First, we reduced the discharge time to 3 s and the equilibration time to 27 s. For this problem, regularization was not necessary. We furthermore skipped the resource constraint such that the problem reduces to (33), (34) and (40). The distribution of material $A$ of the obtained design is shown in Fig. 10a. This design has practically no mixed material and is therefore shown in gray to aid visualization. The objective function only increases from $23.9$ cm3 in the two-layer initial design with material $A$ filling $25\%$ of the domain, to an objective value of $26.7$ cm3 with material $A$ filling $30\%$ of the domain. As in the previous cases, the design shows a bubble of material $A$ below the discharge inlet, but appears simpler in topology. The initial layer of material $A$ at the top remains nearly untouched as it is remains dry within the short simulation time.

Finally, we considered a long discharge of 48 s followed by an equilibration phase of 240 s. We again show the design obtained without resource constraint and without slope constraints in Fig. 10b. There is a centered bubble of material $A$, with a small dent filled with material $C$. The bubble is surrounded by a ring of material $C$, which is itself surrounded by a thin ring of material $A$. Beyond this described structure, the design is in its initial state in regions which remain dry. The objective function increases from 118 to 176 cm3. Solving the same problem with material constraint and with or without slope constraints, yields a similar design. The presented one visualizes the design principle best.

6 Three-Material Optimization

We here focus on the optimization of the reference three-layer design described in Sec. 2.3, where material $A$ is restricted to $5\%$ (43), material $C$ to $85\%$ (44) and the same liquid discharge along with equilibration times as for the standard bi-material problem is used. Here the spatial distribution of the buffer material occupies 10% of the volume and its distribution is allowed to change as well. The problem without regularization and grayness constraints is given by the equations (41) – (44) and (53)–(56). The obtained design features almost no mixed materials and is shown in Fig. 11a. The objective function increases here from 41 to 86 cm3, which is achieved with structures reminiscent of those found in the bi-material optimization, but with important differences. Firstly, because material A (shown as green) is only 5% of the total volume, it cannot be employed to form a large bubble below the discharge inlet. Instead, the bubble of material $A$ is much smaller and is surrounded by a layer along the border of the inlet, which is thinner in the horizontal direction, while penetrating deeper. This allows to distribute liquid efficiently along the horizontal direction. Secondly, the central region of the bubble, which was formed of material $A$ in the bi-material designs, is now made of material $B$ (shown as orange). Although $B$ is less permeable than $A$, it is still much more permeable than $C$.

Adding regularization by slope constraints with $c_{s}=0.3$, we have with the exception of grayness constraints the full problem formulation (41) – (44) and (47)–(56). The enforced material mix is clearly disadvantageous, the design becomes much less distinct as depicted in Fig. 11b and the obtained objective function value reduces to 74 cm3.

6.1 Variants

We carried out numerical simulations with extended problem formulations, which we summarize in what follows. A real diaper must not only absorb as much liquid as possible during the discharge, but should also keep the top layer, which is in contact with the skin, as dry as possible. The function for the moisture content (26) can easily be defined for a subarea $\Omega_{q}$ of the complete computational domain $\Omega$, together with an own adjoint problem to be solved. We defined a top layer of 0.3 mm and limited the amount of fluid within the top layer by a constraint function in the sense of a wetness constraint. The result is that the top layer is filled with material $A$ and partly with $B$ (far away from the discharge inlet) because they have a small moisture capacity. Except for this difference, the optimized design remain similar.

A further point is that several liquid discharges may occur before the diaper is discarded. We performed optimization runs with several discharge and equilibration phases and again no significant impact on the design was observed. Finally, we considered several discharge locations associated with a discharge probability. This requires the separate solution of the state and adjoint problems for each location. For the objective function, the discharges were combined as a weighted sum (with the probability above) and a separate wetness constraint was associated to each discharge case. This complicates the problem formulation further. Material $A$ and $B$ need to be distributed to cover all load locations. The resulting design applies the general design principle according to their weighting factor for the individual load cases and is less clear. As a result, the gains in performance are much more limited, only of about 6%, in comparison to the reference three-layer design.

7 Discussion and Conclusion

A model for gradient based topology optimization of different materials within an unsaturated flow problem was presented. Our benchmark problem was a highly simplified computational model of liquid propagation in a diaper, where the spatio-temporal evolution of the moisture content is obtained from the solution of Richards’ equation. A versatile finite-volume solver was implemented in OpenFOAM for this purpose and was further extended to solve the adjoint equation resulting from the topology optimization. It was shown that gravitational force plays a negligible role in the imbibition process, which is mainly dominated by capillary force.

Our topology optimization yielded designs showing a substantially enhanced capability to absorb liquid with respect to a simple layered design. The optimized designs are similar in all investigated bi-material and three-material problems. Their main common feature is the clustering of highly permeable material close to the discharge inlet in order to allow a fast evacuation of the liquid from the surface. This results in a higher flux of liquid through the inlet and hence a larger total volume is absorbed at the end of the simulation.

All the bi-material optimized designs differ only in subtle details and present a similar performance. This is because of the large volume fraction of the most permeable material $A$ considered (25%). In fact, in most cases there is so much volume of $A$ that parts of it are left on the top of the diaper, in the regions which never hold any moisture. In all cases, a large bubble of material $A$ is formed and reaches from the discharge inlet at the top, to the bottom of the diaper. Because of the large permeability of $A$, these designs enable a rapid distribution of the liquid in the horizontal direction, where the liquid is then stored in the superabsorber (material $C$).

The three-material problem considered here is closer to a real diaper. In this case, only 5% of the total volume is composed of material $A$, 10% of the buffer material $B$ and the remaining 85% of the superabsorber $C$. The large bubble of high permeability below the inlet remains similar in shape and size, but because of the scarcity of material $A$, the bubble is composed not only of material $A$, but also of the buffer material $B$. As expected, $A$ forms a layer immediately below the inlet, which is slightly thicker than the one of the layered reference design. However, in a narrow region at the perimeter of the inlet, this layer of material $A$ protrudes down to the bottom of the diaper. This enhances the distribution of the liquid far away from the inlet in the horizontal direction. Material $B$ covers $A$ completely and reaches down to the bottom of the diaper. Therefore the role of material $B$ in the optimized design is similar to its role in the reference three-layer design, namely to act as a transition region between $A$ and $C$.

Very small volumes of mixed material appeared in the optimized designs. This can be explained by its poor performance. In fact, regions of mixed material occurred only in the transition between pure materials and merely because the volume fraction of each material must be preserved in the optimization. Hence, we expect that with finer discretization, the volume of mixed materials would be even smaller. At the end of the optimization procedure, mixed material was eliminated with a volume-preserving rounding, which resulted in a negligible loss of performance. Regularization of the problem with slope constraints led to worse results.

Finally, we note that some of the fine features of the optimized designs are not obvious. A detailed understanding of the underlying mechanisms and their relationship to the spatial and temporal settings and material properties, is an interesting question in itself, which is beyond the scope of this publication. In addition, the model may be adopted for arbitrary applications of unsaturated flow problems in heterogeneous porous materials. The objective here was to maximize the amount of absorbed fluid, but with the presented model, more complex problems with upper and lower bounds for the amount of absorbed fluids in different regions, or more complex discharge cases, can be easily formulated.

Acknowledgements

The authors gratefully acknowledge the support of the Cluster of Excellence ’Engineering of Advanced Materials’ at the University of Erlangen-Nuremberg, which is funded by the German Research Foundation (DFG) within the framework of its ’Excellence Initiative’.