Topological Optimization of a Generalized Model of Pollution in Porous Media without Density Assumptions

Mouhamadou Baïdy Dia^1*, Ibrahima Sagno¹, Lansana Toure¹, Mouctar Ndiaye¹, Alassane Sy², Ibrahima Faye²

¹Laboratoire de Mathématiques Appliquées (LABOMA), Faculté des Sciences de la Nature, Université Julius Nyerere de Kankan, Kankan, Guinée.
²Laboratoire d’Informatique, de Mathématiques et Applications (LIMA), UFR S.A.T.I.C, Université Alioune Diop de Bambey, Bambey, Sénégal.

DOI: 10.4236/ojop.2025.143006                

Abstract

This study focuses on the numerical topological gradient applied to a model of pollutant transfer in porous media, obtained by removing assumptions regarding critical parameters such as density and conductivity. A generalized model of pollutant transport dynamics in porous media was developed without any assumptions regarding the density and conductivity. Topological optimization is then applied to minimize the least-squares function associated with the system representing the model. This study demonstrates the effectiveness of this approach in both two-dimensional (2D) and three-dimensional (3D) configurations, providing insights into the optimal structural design of porous media to mitigate pollution transfer.

Keywords

Topological Gradient, Pollutant Transfer, Porous Media, Least-Squares Functional, Topological Optimization, Topological Gradient Descent of Direction

1. Introduction

Understanding pollutant transfer in porous media is a major challenge in environmental modeling, with applications in soil decontamination, groundwater management, and filtration processes. Classical models of pollutant transport often assume that key physical parameters, such as density and conductivity, remain constant [1] [2]. Although these assumptions simplify theoretical and numerical analyses, they fail to capture the complexity of heterogeneous porous media observed in real-world environments.

In our previous work [3] [4], we extended these classical models by incorporating spatial variability in density and conductivity, distinguishing between cases in which their divergence is zero or not. This approach provided a more realistic representation of pollutant transport but still imposed structural constraints on the physical properties of the medium.

This study further generalizes our previous approach by removing all restrictive assumptions regarding density and conductivity. Unlike previous studies that treated these parameters as constant or variable in some assumptions [1] [2], our model fully accounted for their spatial variability with no restriction, offering a more flexible and comprehensive framework. In this model, we applied topological optimization techniques, specifically the topological gradient method, to identify optimal configurations that minimize pollutant dispersion in porous structures.

First, we formulate a generalized mathematical model for pollutant transport, highlighting its flexibility and applicability to a wide range of scenarios. Next, we define the least-squares function as the target for minimization, ensuring a rigorous optimization framework. We then present the theory and implementation of the topological gradient method and demonstrate its effectiveness in exploring complex solution spaces. Finally, we provide numerical simulations in both 2D and 3D settings, illustrating the efficiency of our approach for optimizing pollutant transport structures.

This study highlights the potential of topological optimization for tackling complex environmental challenges by seamlessly integrating mathematical modeling, computational optimization, and numerical simulation.

The work proceeds as follows: In Section 2, we first present the model of pollutant transportation in porous media without any assumptions made on key parameters such as density and conductivity. Section 3 summarizes the topological optimization method and its application to the model of pollutant transportation in porous media. Section 4 describes the proposed numerical method and numerical simulations in both 2D and 3D. Section 5 provides some concluding remarks and possible extensions.

2. Modelization of Pollutant Transfer in Porous Media Without any Assumption on Density

The pollution model in porous media is commonly described by several laws governing the interaction of fluid motion through porous media. For more details, see [1] [5]-[10].

2.1. Equations Governing Pollutant Transfer

The transfer of pollutants in porous media can be described [1] [3] by the following system of equations:

$\{\begin{array}{l}\frac{\partial \left(\epsilon d\right)}{\partial t}+\nabla \cdot \left(dq\right)=0\text{ (2}\text{.1a)}\\ \frac{\partial \left(\epsilon d\mathcal{W}\right)}{\partial t}+\nabla \cdot \left(d\mathcal{W}q+J\right)=0\text{ (2}\text{.1b)}\\ J=-dD\nabla \mathcal{W}\left(\text{Fick's law}\right)\text{ (2}\text{.1c)}\\ d=d_0\exp \left({\beta }_T\left(T-T_0\right)+{\beta }_P\left(P-P_0\right)+\gamma \mathcal{W}\right),\left[1,6\right]\text{ (2}\text{.1d)}\\ q=-\frac{k}{\mu }\left(\nabla P+dg{e}_3\right)\left(\text{Darcy's law}\right)\text{ (2}\text{.1e)}\end{array}$

There were several important variables and parameters in the model. The flow density $d\left(\text{kg}/{\text{m}}^3\right)$ is a function of the temperature $\left(T\right)$ , pressure $\left(P\right)$ and the mass pollutant ratio transported $\mathcal{W}$ . The porosity of the porous medium, $\left(\epsilon \right)$ , is a dimensionless value from 0 to 1, which is the ratio of the void volume to the medium volume. The fluid vector indicates $q\left(\text{m}/\text{s}\right)$ and indicates the Darcys’ velocity, which is a function of the pressure of the fluid $P\left(Pa\right)$ . The pollutant mass transport is diffusive owing to the diffusion coefficient $D\left({\text{m}}^2/\text{s}\right)$ and is expressed in the diffusive flux $J\left(\text{kg}/\left({\text{m}}^2/\text{s}\right)\right)$ , that follows Fick’s law. Time $t\left(\text{s}\right)$ as the execution of the system continues. Some reference values for temperature $T_0$ and pressure $P_0$ and sensitivity coefficients $\left({\beta }_T\right),\left({\beta }_P\right)$ , and $\gamma$ for the temperature, pressure and $\left(\mathcal{W}\right)$ , affect the fluid density $\left(d\right)$ . The permeability of the porous medium $\left(\frac{k}{\mu }\right)$ represents the medium’s ability to allow fluid transmission, and ${ }^{m}u$ fluid dynamic viscosity (it’s value is in (Pa∙s)) means the degree of the fluid flow resistance. Gravitational effects are introduced by means of gravitational acceleration $g\left({\text{m}}^2/\text{s}\right)$ and unit vector $e_3=\left[0;0;1\right]$ that indicates how gravity acts in the model.

2.2. Generalized Model of Pollution in Porous Media

In this approach, we avoid making assumptions about density d. Instead, which is treated as a variable influenced by external factors such as pressure, temperature, and pollutant concentration. This generalization allowed the model to remain applicable to a wider range of porous media and scenarios.

The system of Equations (2.1) can be rewritten according to [1]-[3] as:

$\begin{array}{c}\frac{\partial \left(\epsilon d\right)}{\partial t}-\nabla \cdot \left(\frac{k}{\mu }\nabla Pd\right)-\nabla \cdot \left(\frac{k}{\mu }{d}^{2}g{e}_3\right)=0,\\ \frac{\partial \left(\epsilon d\mathcal{W}\right)}{\partial t}-\nabla \cdot \left(d\mathcal{W}\frac{k}{\mu }\nabla P\right)-\nabla \cdot \left(d^2\mathcal{W}\frac{k}{\mu }g{e}_3\right)-\nabla \cdot \left(dD\nabla W\right)=0\end{array}$ (2.2)

REMARK 1 The porosity ε of the medium can be modeled using Garner’s law:

$\epsilon =\{\begin{array}{ll}\frac{{\epsilon }_s-{\epsilon }_r}{1+{\left(\alpha h\right)}^{\beta }},\hfill & h\le 0\hfill \\ {\epsilon }_s,\hfill & h>0\hfill \end{array}$

where $h$ represents the pressure relative to the atmospheric pressure.

By developing the second equation in (2.2), we have

$\mathcal{W}\left[\frac{\partial \left(\epsilon d\right)}{\partial t}-\nabla \cdot \left(\frac{k}{\mu }\left(\nabla P\right)d_s\right)-\nabla \cdot \left(\frac{k}{\mu }{d}^{2}g{e}_3\right)\right]+\epsilon d\frac{\partial \mathcal{W}}{\partial t}-d\frac{k}{\mu }\nabla \mathcal{W}\cdot \nabla P-d^2\frac{k}{\mu }g\nabla \mathcal{W}\cdot e_3-\nabla \cdot \left(dD\nabla \mathcal{W}\right)=0$ (2.3)

According to the first equation of (2.2), we have:

$\epsilon d\frac{\partial \left(\mathcal{W}\right)}{\partial t}-d\frac{k}{\mu }\nabla \mathcal{W}\nabla P-d^2\frac{k}{\mu }g{\partial }_z\mathcal{W}-\nabla \cdot \left(dD\nabla \mathcal{W}\right)=0$

Due to (2.1d) and the temperature does not vary, one gets

$\begin{array}{c}\nabla P=-\frac{1}{\beta }\nabla W+\frac{1}{\gamma {\beta }_p}\nabla \left(\ln \left(\frac{d}{d_0}\right)\right)\end{array}$ (2.4)

$\epsilon d\frac{\partial \mathcal{W}}{\partial t}-d\frac{k}{\mu }\nabla \mathcal{W}\cdot \left[-\frac{1}{\beta }\nabla \mathcal{W}+\frac{1}{\gamma {\beta }_p}\nabla \left(\text{ln}\left(\frac{d}{d_0}\right)\right)\right]-d^2\frac{k}{\mu }g{\partial }_z\mathcal{W}-\nabla \cdot \left(d D\nabla \mathcal{W}\right)=0$

It follows that:

$\epsilon d\frac{\partial \mathcal{W}}{\partial t}+\frac{dk}{\mu \beta }{\left|\nabla \mathcal{W}\right|}^2-\frac{dk}{\mu \gamma {\beta }_p}\nabla \left(\text{ln}\left(\frac{d}{d_0}\right)\right)\cdot \nabla \mathcal{W}-d^2\frac{k}{\mu }g{\partial }_z\mathcal{W}-\nabla \cdot \left(d D\nabla \mathcal{W}\right)=0$

By rearranging, we have

$\frac{\partial \mathcal{W}}{\partial t}+\frac{k}{\epsilon \mu \beta }{\left|\nabla \mathcal{W}\right|}^2-\frac{k}{\mu \epsilon \gamma {\beta }_p}\nabla \left(\ln \left(\frac{d}{d_0}\right)\right)\nabla \mathcal{W}-\frac{kd}{\epsilon \mu }g{\partial }_z\mathcal{W}-\frac{\nabla \left(dD\right)}{d\epsilon }\nabla \mathcal{W}-\frac{D}{\epsilon }\text{Δ}\mathcal{W}=0$

we have

$\frac{\partial \mathcal{W}}{\partial t}+\frac{k}{\epsilon \mu \beta }{\left|\nabla \mathcal{W}\right|}^2-\left[\frac{k}{\mu \epsilon \gamma {\beta }_p}\nabla \left(\text{ln}\left(\frac{d}{d_0}\right)\right)+\frac{\nabla \left(dD\right)}{d\epsilon }\right]\nabla \mathcal{W}-\frac{kd}{\epsilon \mu }g{\partial }_z\mathcal{W}-\frac{D}{\epsilon }\text{Δ}\mathcal{W}=0$

Setting ${\alpha }_1=\frac{k}{\epsilon \mu \beta },\overline{{\alpha }_2}=\frac{k}{\mu \epsilon \gamma {\beta }_p}\nabla \left(\text{ln}\left(\frac{d}{d_0}\right)\right)+\frac{\nabla \left(dD\right)}{d\epsilon },{\alpha }_3=\frac{kd}{\epsilon \mu }g$, and ${\alpha }_4=\frac{D}{\epsilon }$ ,

one has

$\frac{\partial \mathcal{W}}{\partial t}+{\alpha }_1{\left|\nabla \mathcal{W}\right|}^2-\overline{{\alpha }_2}\cdot \nabla \mathcal{W}-{\alpha }_3{\partial }_z\mathcal{W}-{\alpha }_4\text{Δ}\mathcal{W}=0$

We aim to fix the boundary and initial conditions to guarantee the good posedness of the pollution model [11] [12]. Thus

$\begin{array}{c}\{\begin{array}{l}\frac{\partial \mathcal{W}}{\partial t}+{\alpha }_1{\left|\nabla \mathcal{W}\right|}^2-\overline{{\alpha }_2}\cdot \nabla \mathcal{W}-{\alpha }_3{\partial }_z\mathcal{W}-{\alpha }_4\text{Δ}\mathcal{W}=0,\text{ in }ℬ\times \left(0,T_1\right),\hfill \\ \mathcal{W}={\mathcal{W}}_1{\text{ on Γ}}_1\times \left(0,T_1\right),\hfill \\ \frac{\partial \mathcal{W}}{\partial n}=0\text{ on }\partial ℬ\setminus {\text{Γ}}_1\times \left(0,T_1\right),\hfill \\ \mathcal{W}\left(x,0\right)={\mathcal{W}}_i,x\in ℬ\times \left\{0\right\}.\hfill \end{array}\end{array}$ (2.5)

REMARK 2 The pollution model we obtain is a generalization of the models obtained in [1]-[3]. In fact, assuming the same hypotheses used in [3] of our generalized model, we obtain the same pollution model.

3. Topological Optimization

We applied the topological optimization tools to the model of pollution in porous media that we obtained. Among several methods used to compute the topological derivative [13]-[20], we adopted the generalized adjoint method developed in [13] [21].

3.1. A Generalized Adjoint Method

More details of the generalized adjoint method can be found in [22] [23].

Let $\mathscr{H}$ be a fixed Hilbert space and $ℒ\left(\mathscr{H}\right)\left(\text{resp}{ℒ}_2\left(\mathscr{H}\right)\right)$ denotes the spaces of linear (resp bilinear) forms on $\mathscr{H}$ . Let us assume the following hypotheses.

HYPOTHESIS 1: There exists a real function $f$ , a bilinear form ${\delta }_a\in {ℒ}_2\left(\mathscr{H}\right)$ and a linear form ${\delta }_l$ such that

$f\left(r\right)\to 0,r\to 0$ (3.1)

${\Vert a_r-a_0-f\left(r\right){\delta }_a\Vert }_{L_2\left(\mathscr{H}\right)}=o\left(f\left(r\right)\right)$ (3.2)

${\Vert l_r-l_0-f\left(r\right){\delta }_l\Vert }_{L\left(\mathscr{H}\right)}=o\left(f\left(r\right)\right)$ (3.3)

where $a_r,a_0$ resp ($l_r,l_0$ ) denotes the bilinear resp (linear) form.

HYPOTHESIS 2: The bilinear form $a_0$ is coercive. There exists a constant $\alpha >0$ such that

$\begin{array}{c}{a}_0\left(u,u\right)\ge \alpha {\Vert u\Vert }^2,\forall u\in \mathscr{H}\end{array}$ (3.4)

According to (3.2) the bilinear form $a_r$ depends continuously on $r$ , hence there exists $r_0$ and $\beta >0$ such that for $r\in \left[0,r_0\right]$ the following uniform coercive holds $a_r\left(u,u\right)\ge \beta {\Vert u\Vert }^2,\forall u\in \mathscr{H}$ . Moreover, according to Lax-Milligram’s theorem, for $r\in \left[0,r_0\right]$ , the problem finds $u_r\in \mathscr{H}$ , such that

$\begin{array}{c}{a}_r\left(u_r,v\right)=l_r\left(v\right),\forall v\in \mathscr{H}\end{array}$ (3.5)

has a unique solution and the following result holds.

LEMMA 1 If Hypotheses 1 and 2 are satisfied, there exists a unique solution $u_r$ to (3.5). Furthermore, the following estimate holds:

${\Vert u_r-u\Vert }_{{\mathscr{H}}_r}\begin{array}{c}=O\left(f\left(r\right)\right)\end{array}$ (3.6)

where $u$ is the solution of: $a_0\left(u,v\right)=l_0\left(v\right),\forall v\in \mathscr{H}$ .

Proof. Because of the ellipticity of $a_r$ we obtain

$\alpha {\Vert u_r-u_0\Vert }^2\le a_r\left(u_r-u_0,u_r-u_0\right)$

From

$\begin{array}{l}{a}_r\left(u_r-u_0,u_r-u_0\right)=\left|a_r\left(u_r,u_r-u_0\right)-a_r\left(u_0,u_r-u_0\right)\right|\\ =\left|l_r\left(u_r-u_0\right)-a_r\left(u_0,u_r-u_0\right)\right|\\ =\left|l_0\left(u_r-u_0\right)+\left(l_r-l_0\right)\left(u_r-u_0\right)-a_r\left(u_0,u_r-u_0\right)\right|\\ \le \Vert a_0\left(u_0,u_r-u_0\right)-a_r\left(u_0,u_r-u_0\right)\Vert +\Vert \left(l_r-l_0\right)\left(u_r-u_0\right)\Vert \end{array}$

Thus

$\left|\alpha \right|{\Vert u_r-u_0\Vert }^2\le \Vert \delta a\left(u_0,u_r-u_0\right)\Vert f\left(r\right)+\Vert \delta l\left(u_r-u_0\right)\Vert f\left(r\right)+\Vert u_0\Vert \Vert u_r-u_0\Vert o\left(f\left(r\right)\right)+\Vert u_r-u_0\Vert o\left(f\left(r\right)\right)$ .

Thus,

$\alpha {\Vert u_r-u_0\Vert }^2\le \left(\Vert \delta a\Vert \Vert u_0\Vert +\Vert \delta l\Vert \right)\Vert u_r-u_0\Vert f\left(r\right)+\left(1+\Vert u_0\Vert \right)\Vert u_r-u_0\Vert o\left(f\left(r\right)\right)$ .

which implies that

$\alpha \Vert u_r-u_0\Vert \le \left(\Vert \delta a\Vert \Vert u_0\Vert +\Vert \delta l\Vert \right)f\left(r\right)+\left(1+\Vert u_0\Vert \right)o\left(f\left(r\right)\right)$ ,

this confirms the proof. ■ We also assume the following hypothesis.

HYPOTHESIS 3: Let $j$ be a functional such that $j\left(r\right)=J\left(u_r\right)$ where $J$ is differentiable. For any $u\in \mathscr{H}$ , there exists a linear and continuous form $DJ\left(u\right)\in ℒ\left(\mathscr{H}\right)$ and ${\delta }_J$ such that:

$\begin{array}{c}J\left(u\right)-J\left(v\right)=DJ\left(u\right)\left(v-u\right)+f\left(r\right){\delta }_J\left(u\right)+o\left({\Vert u-v\Vert }_{\mathscr{H}}\right)\end{array}$ (3.7)

The Lagrangian $ℒ$ is given by

$\begin{array}{c}L\left(u,v\right)=a\left(u,v\right)-l\left(v\right)+J\left(u\right),\forall u,v\in \mathscr{H}\end{array}$ (3.8)

Its variation is given by

$\begin{array}{c}{ℒ}_r\left(u,v\right)=a_r\left(u,v\right)-l_r\left(v\right)+J\left(u\right),\forall u,v\in \mathscr{H}\end{array}$ (3.9)

We also have the following generic theorem.

THEOREM 1 Under hypotheses (1), (2) and (3), we have the following asymptotic expansion of $j$ :

$\begin{array}{c}j\left(r\right)-j\left(0\right)=f\left(r\right)\delta L\left(u_0,v_0\right)+o\left(f\left(r\right)\right)\end{array}$ (3.10)

where $u_0$ is the solution of (3.5) with $r=0,v_0$ is the solution to the adjoint problem, and $v_0$ is such that

$\begin{array}{c}{a}_0\left(ℬ,v_0\right)=-DJ\left(u_0\right)\int ,\forall f\in \mathscr{H}\end{array}$ (3.11)

and $\delta ℒ\left(u,v\right)={\delta }_a\left(u,v\right)-{\delta }_l\left(v\right)+{\delta }_J\left(u\right)$ , denotes the topological derivative.

Proof. For the proof of this theorem, see [23]. ■

3.2. Application of Topological Optimization to the Generalized Model of Pollution in Porous Media

To optimize the configuration of the porous media for pollution control, we employed topological optimization methods. This approach aims to minimize the function representing the pollutant mass ratio while satisfying physical constraints.

For a real number $r\ge 0$ , we define the perturbed domain ${\text{Ω}}_r=\text{Ω}\setminus \overline{{ℬ}_r}$, where ${ℬ}_r=x_0+rℬ$ and $ℬ\subset {ℝ}^n$ , with $n=2$ or 3.

We consider the sets $Q=\left(0,T\right)\times \text{Ω},{\displaystyle {\sum }_{{}_r}=}\left(0,T\right)\times {ℬ}_r$ , and $Q_r=\left(0,T\right)\times {\text{Ω}}_r$ .

The optimization process involves minimizing the objective functional $J\left({\mathcal{W}}_r\right)$ , defined as

where ${\mathcal{W}}_r$ is the solution to the perturbed problem (3.12):

$\begin{array}{c}\{\begin{array}{ll}\frac{\partial {\mathcal{W}}_r}{\partial t}+{\alpha }_1{\left|\nabla {\mathcal{W}}_r\right|}^2-\overline{{\alpha }_2}\cdot \nabla {\mathcal{W}}_r-{\alpha }_3\frac{\partial {\mathcal{W}}_r}{\partial z}-{\alpha }_4\text{Δ}{\mathcal{W}}_r=0,\hfill & {\text{in Ω}}_r\times \left(0,T\right),\hfill \\ {\mathcal{W}}_r={\mathcal{W}}_1,\hfill & {\text{on Γ}}_1\times \left(0,T\right),\hfill \\ \frac{\partial {\mathcal{W}}_r}{\partial n}=0,\hfill & \text{on}\left(\partial {\text{Ω}}_r\setminus {\text{Γ}}_1\right)\times \left(0,T\right),\hfill \\ B_{ℬ}{\mathcal{W}}_r=0,\hfill & \text{on }\partial ℬ,\hfill \\ {\mathcal{W}}_r\left(x,0\right)={\mathcal{W}}_0,\hfill & \text{for }x\in \text{Ω}\times \left\{0\right\}.\hfill \end{array}\end{array}$ (3.12)

Here, $B_{ℬ}$ is an operator defined as $\partial ℬ$ , mapping ${\mathcal{W}}_r$ to either $\mathcal{W}$ or ${\partial }_n\mathcal{W}$ , and ${\mathcal{W}}_d$ is the desired pollutant mass reference.

To simplify the problem, we introduce a change of variable:

$u=\psi \left(\mathcal{W}\right)=\lambda e^{\frac{{\alpha }_1}{{\alpha }_4}\mathcal{W}}$

This leads to the transformed system (3.13):

$\begin{array}{c}\{\begin{array}{ll}\frac{\partial u_r}{\partial t}-\overline{{\alpha }_2}\cdot \nabla u_r-{\alpha }_3\frac{\partial u_r}{\partial z}-{\alpha }_4\text{Δ}{u}_r=0,\hfill & {\text{in Ω}}_r\times \left(0,T\right),\hfill \\ u_r=\psi \left({\mathcal{W}}_1\right),\hfill & {\text{on Γ}}_1\times \left(0,T\right),\hfill \\ \frac{\partial u_r}{\partial n}=\psi \left({\mathcal{W}}_2\right),\hfill & \text{on }\left(\partial {\text{Ω}}_r\setminus {\text{Γ}}_1\right)\times \left(0,T\right),\hfill \\ B_{ℬ}{u}_r=0,\hfill & \text{on }\partial ℬ,\hfill \\ u_r\left(x,0\right)=\psi \left({\mathcal{W}}_0\right),\hfill & \text{for }x\in \text{Ω}\times \left\{0\right\}.\hfill \end{array}\end{array}$ (3.13)

The optimization problem then becomes solving:

where $u_r$ is the solution of (3.13) which weak formulation is to find $u_r\in ℱ\left(\text{Ω}\right)$ such that:

for all $v\in \mathcal{D}\left(\text{Ω}\right)$ and

$ℱ\left(\text{Ω}\right)=\left\{u_r\in H^1\left(\text{Ω}\right),u_r=u_1{\text{ on Γ}}_1,\frac{\partial u_r}{\partial n}=0{\text{ on Γ}}_2\text{ and }{B}_{{ℬ}_r}{u}_r=0,\text{ on }\partial {ℬ}_r\right\}.$

The bilinear and linear forms associated with (3.15) are

(3.14)

and

and

For $r=0$ , the objective functional simplifies to:

where $u_{\text{Ω}}$ is the solution to the following problem (3.15):

$\begin{array}{c}\{\begin{array}{ll}\frac{\partial u_{\text{Ω}}}{\partial t}-\overline{{\alpha }_2}\nabla u_{\text{Ω}}-{\alpha }_3\frac{\partial u_{\text{Ω}}}{\partial z}-{\alpha }_4\text{Δ}{u}_{\text{Ω}}=0,\hfill & \text{in Ω}\times \left(0,T\right)\hfill \\ u_{\text{Ω}}=\psi \left({\mathcal{W}}_1\right),\hfill & {\text{on Γ}}_1\times \left(0,T\right)\hfill \\ \frac{\partial u_{\text{Ω}}}{\partial n}=\psi \left({\mathcal{W}}_2\right),\hfill & \text{on }\left(\partial \text{Ω}\setminus {\text{Γ}}_1\right)\times \left(0,T\right)\hfill \\ u_{\text{Ω}}\left(x,0\right)=\psi \left({\mathcal{W}}_0\right),\hfill & \text{for }x\in \text{Ω}\times \left\{0\right\}\hfill \end{array}\end{array}$ (3.15)

The weak formulation of (3.15) is to find $u\in ℱ\left(\text{Ω}\right)$ such that:

for all $v\in \mathcal{D}\left(\text{Ω}\right)$ and

$ℱ\left(\text{Ω}\right)=\left\{u\in H^1\left(\text{Ω}\right),u=u_1{\text{ on Γ}}_1\text{ and }\frac{\partial u}{\partial n}=0{\text{ on Γ}}_2\right\}$ .

The bilinear and linear forms associated with (3.15) are

(3.16)

and

Under these conditions, we obtain the following results.

LEMMA 1 If $u_0$ and $u_r$ are the respective solutions to problems (3.16) and (3.14), then there exists $f\left(r\right)\ge 0$ , with , such that

$\begin{array}{c}{u}_r-u_0{}_{ℱ\left(\text{Ω}\right)}=O\left(f\left(r\right)\right)\end{array}$ (3.17)

Proof. See [3] ■

Thus, we have the following result:

THEOREM 2 (Main Result) Under the conditions of Lemma 1, the topological derivative of functional (3.2), where $ℬ$ denotes the unit ball centered at $x_0$ , is given by

(3.18)

for Neumann boundary conditions

, for Dirichlet boundary conditions.

State $u_{\Omega }$ solves the direct problem (3.13) and $v_{\Omega }\left(x_0,t\right)$ is the solution of the corresponding adjoint problem.

Proof. To prove this theorem we determine the variation of the bilinear form, objective function and linear form. Subsequently Theorem 2 is applied.

Variation of the bilinear form

The variation of the bilinear form $a_r-a_0$ associated with problem (3.14) and (3.16) is given by

However, from Lemma 1, the following estimate holds

.

Thus,

(3.19)

in the case of a Neumann condition, since $ℬ\in {ℝ}^n,n=2$ or 3, we take $f\left(r\right)=r^n$ .

It follows by applying the change of variables (3.2) that we have

Using the mean value theorem, we obtain the following result:

Since

,

it follows that

.

In the case of a Dirichlet condition at the boundary of the hole, a similar reasoning applies.

Variation of the functional form

We have

Let us denote:

(3.20)

and

(3.21)

We then have:

(3.22)

It follows that:

(3.23)

From Lemma 1, we know that $u_r-u_{\text{Ω}}{}_{ℱ\left(ℬ\right)}=O\left(f\left(r\right)\right)$ . Hence:

${\psi }^{-1}\left(u_r\right)-{\psi }^{-1}{\left(u_{\text{Ω}}\right)}_{ℱ\left(ℬ\right)}=O\left(f\left(r\right)\right),$

since $\psi$ is a bijective and continuous function.

Thus:

(3.24)

We therefore have:

(3.25)

To complete the proof, it suffices to show that:

(3.26)

• For the Neumann boundary condition around the of the hole, applying the variable change

$\begin{array}{c}C:B\subset {ℝ}^n\to {ℬ}_r\subset {ℝ}^n,\\ y\mapsto C\left(y\right)=x=x_0+ry,\end{array}$ (3.27)

we get:

(3.28)

Rearranging gives:

$$

By setting $f\left(r\right)=r^n$ and applying the mean value theorem, we obtain:

Since

the variation of the functional is

Therefore, when $ℬ$ represents the unit ball centered at $x_0$ , we have:

• For the Dirichlet boundary condition around the of the hole, using a Taylor expansion near $x_0$ , we obtain:

For the Dirichlet condition, it is well known that $f\left(r\right)=r$ in three dimensions. Thus, ${ℐ}_2=o\left(f\left(r\right)\right)$ , and the variation becomes:

$\delta J\left(u_{\text{Ω}}\left(x,t\right)\right)=0.$

According to Theorem 2, the expression for the topological derivative is:

(3.29)

for Neumann boundary conditions, or:

for Dirichlet boundary conditions.

4. Numerical Results

The simulation of the pollutant dynamics was performed with FreeFem++ [24], which allows the handling of different types of boundary conditions (Neumann and Dirichlet). The states and $v_{\text{Ω}}$ represent the solutions of the direct and adjoint problems, respectively. The time step $\text{Δ}t=0.1$ and the time $\tau =T_{\text{max}}=50$ were used.

The topological gradient descent algorithm seems to be a powerful tool for the analysis of several physical problems such as pollution [4] [25]. Therefore, we can write Algorithm (1) for the optimal design minimizing the least-squares functional.

Algorithm 1 (Topological Optimization Algorithm)

1) Initialization:

• Choose an initial domain ${\text{Ω}}_0$ ,

• Set a small threshold value $r_0>0$ and a radius $R_k$ ,

• Initialize the iteration index: $k=0$ .

2) Iterations: Repeat

• Solve the direct and adjoint problems to compute $u_{{\text{Ω}}_k}$ and $v_{{\text{Ω}}_k}$ in the domain ${\text{Ω}}_k$ .

• Compute the topological gradient $g_k\left(x\right)$ for all $x\in {\text{Ω}}_k$ .

• Identify the minimum of $g_k\left(x\right)$ , denoted $g_k\left(x_k^{\text{*}}\right)$ , and its corresponding location $x_k^{\text{*}}$ within ${\text{Ω}}_k$ .

• While $\left|\text{min}\left(g_k\left(x\right)\right)\right|\ge r_0$ , perform:

• Update domain ${\text{Ω}}_{k+1}={\text{Ω}}_k\setminus B\left(x_k^{\text{*}},R_k\right)$ , where $B\left(x_k^{\text{*}},R_k\right)$ is the ball centered at $x_k^{\text{*}}$ and radius $R_k$ .

• Solve the updated direct and adjoint problems in ${\text{Ω}}_{k+1}$ to compute $u_{{\text{Ω}}_{k+1}}$ and $v_{{\text{Ω}}_{k+1}}$ .

• Recompute the topological gradient: $g_{k+1}\left(x\right)$ for all $x\in {\text{Ω}}_{k+1}$ .

• Increment the iteration index: $k=k+1$ .

3) Finalization: Print the obtained optimal design.

In this section, we numerically apply topological gradient descent to the general pollution model. To discretize this partial differential equation with an implicit scheme (backward Euler), we use the relation:

(4.1)

where $\mathcal{A}\left(u,v\right)$ is the bilinear operator given by:

(4.2)

The weak formulation at the time step $n+1$ becomes:

(4.3)

The problem is solved iteratively over the time steps, where $u^n$ represents the solution at step n.

The topological gradient ${\mathcal{G}}_{\tau }\left(x_0\right)$ is given by:

(4.4)

in the Neumann condition and

for Dirichlet boundary conditions.

4.1. Numerical Results in Two Dimensions Case with a Dirichlet Boundary Condition on the Hole $\partial ℬ$

In two dimensions, we consider the computational domain $\text{Ω}=\left(0,H\right)\times \left(0,L\right)$ with $H=10$ and $L=10$ . The simulations were performed using the following parameters: the permeability is given by $\kappa =0.03y^2+0.3x^2+0.2$ , while the solid density is $d_s=y^2+30x^2+20$ . The diffusion coefficient is defined as $D=0.1(3y^2+30x^2+20)$ , and the drift velocity is defined $U_d=\left|\text{cos}\left(x\right)\text{sin}\left(x\right)\right|+t\left(x^2+y^2+1\right)$ . Additionally, we set $d_0=100$ , the initial condition as $u_i=100\left|\text{sin}\left(xy\right)\right|$ , and the boundary condition as $u_1=2\text{cos}\left(x\right)\text{sin}\left(x\right)\text{exp}\left(x^2+y^2\right)+1$ .

According to these data we obtain Figure 1 that represents the direct and adjoin states, and the topological derivative in the Domain.

Figure 1. 2D pollution problem (Dirichlet condition) at initial step. After 35 iterations, the topological derivative is almost null everywhere in the perforated domain. Then we achieve the following optimal domain in Figure 2.

Figure 2. 2D pollution problem (Dirichlet condition) after 35 iterations.

4.2. Numerical Results in Three Dimensions with Neumann Condition

In 3-dimensional cases, we consider the initial domain ${\text{Ω}}_0=\left(-5,5\right)\times \left(-5,5\right)\times$ $\left(0,6\right)$ . The simulations were performed under the following parameter values: the permeability is given by $\kappa =7x^2+3y^2+4z^2+1$ , while the density is $d_s=\frac{100}{7x^2+3y^2+4z^2+1}$ . The diffusion coefficient is given as $D=0.0003\left(7x^2+3y^2+4z^2+1\right)$ , and the desired state is $u_d=0.5\text{cos}\left(x\right)+0.5\text{sin}\left(y\right)+0.5\text{sin}\left(z\right)$ . In addition, we set the initial condition as $u_i=0.5\text{cos}\left(x\right)+0.5\text{sin}\left(y\right)+0.5\text{sin}\left(z\right)$ , and the homogen boundary condition as $u_1=x^2{(x-1)}^2+2+t\left(-2+12x-11x^2-2x^3+x^4\right)$ .

Figure 3. 3D pollution transportation at the first iteration.

Figure 4. 3D pollution transportation after 25 iterations.

Based on these data, Figure 3 shows the distributions of pollutant mass fraction and topological gradient in the perforated domain by inserting the spherical hole $ℬ$ centered at the minimum of ${\mathcal{G}}_{\tau }$ . We notice that creating a spherical hole, at the negative minimum of the topological derivative, decreases the objective function.

Thus, after 25 iterations, the topological gradient is almost null throughout the perforated domain (see Figure 4). The gradient descent algorithm converges, achieving the optimal domain.

5. Conclusions

This study presents a generalized framework for modeling and optimizing pollutant transfer in porous media without imposing restrictive assumptions on density or other physical parameters. By leveraging an optimization-based approach, we demonstrated that pollutant dispersion can be effectively controlled while maintaining flexibility in the governing equations. The key contributions of this study is the successful application of our methodology to both two-dimensional and three-dimensional domains, showcasing its robustness and adaptability to real-world scenarios.

Through three-dimensional simulations, we validated the effectiveness of the proposed approach in capturing the complex behavior of pollutant transfer in porous media. These simulations highlight how topological optimization techniques can identify optimal pollutant diffusion patterns and improve environmental management strategies. The numerical results confirm that our model can handle intricate spatial variations and provide valuable insights into the pollution control in heterogeneous porous structures.

Future work will involve further optimization of the numerical scheme, to improve the computational efficiency and address large-scale 3D problems. In addition, we sought to integrate real data with our models to enhance their utility in real-world applications, such as industrial-waste processing, groundwater contamination and soil bioremediation scenarios. We will also combine field measurements with established experimental data to quantitatively test these theoretical predictions and to develop the best control strategies.

Furthermore, future work will include cutting-edge numerical techniques, such as machine learning and data-driven optimization, which could improve the quality of the model to facilitate the control of pollutants in porous media and allow better strategies of environmental protection, detection, and exploitation of non-renewable energies.

Conflicts of Interest

The authors declare that they have no conflict of interest.

References

[1]	Faye, I., Sy, A. and Seck, D. (2008) On Topological Optimization and Pollution in Porous Media. In: Konaté, D., Ed., Mathematical Modeling, Simulation, Visualization and e-Learning, Springer, 209-237. https://doi.org/10.1007/978-3-540-74339-2_14
[2]	Ndiaye, L., Sy, A. and Seck, D. (2012) Pollution in Porous Media: Non Permanent Cases. Journal of Computations and Modelling, 2, 33-51.
[3]	Dia, M.B., Faye, I. and Sy, A. (2022) Topological Optimization Applied to a Variable Density and Conductivity Model of Pollution in Porous Media. Engineering Reports, 5, e12557. https://doi.org/10.1002/eng2.12557
[4]	Baidy Dia, M. (2023) Étude des cristaux photoniques, phononiques et des problèmes de transfert de polluants en milieux poreux par les méthodes d’optimisation topologique. UCAD.
[5]	Abdelfettah, F. (1995) Modélisation de la propagation de polluants dans les milieux poreux saturés et non saturés. Ph.D. Thesis, Institut National Polytechnique de Lorraine.
[6]	Gujisaite, V. (2008) Transport réactif en milieux poreux non saturés. Master’s Thesis, Institut National Polytechnique de Lorraine, Nancy université.
[7]	Schotting, R.J. (1998) Mathematical Aspects of Sault Transport in Porous Media. Ph.D. Thesis, Techniische Universiteit Delft.
[8]	Tlili, I., Shahmir, N., Ramzan, M., Kadry, S., Kim, J., Nam, Y., et al. (2020) A Novel Model to Analyze Darcy Forchheimer Nanofluid Flow in a Permeable Medium with Entropy Generation Analysis. Journal of Taibah University for Science, 14, 916-930. https://doi.org/10.1080/16583655.2020.1790171
[9]	Ramzan, M., Bilal, M., Chung, J.D. and Farooq, U. (2016) Mixed Convective Flow of Maxwell Nanofluid Past a Porous Vertical Stretched Surface—An Optimal Solution. Results in Physics, 6, 1072-1079. https://doi.org/10.1016/j.rinp.2016.11.036
[10]	Singh, M., Singh, C. and Dasaroju, G. (2016) Modelling for Flow through Unsaturated Porous Media with Constant and Variable Density Conditions Using Local Thermal Equilibrium. 4th International Conference on Advancements in Engineering and Technology (ICAET-2016), Sangrur, 18-19 March 2016, 24-30.
[11]	Brezis, H. (1987) Analyse Fonctionnelle: Théorie et application. Masson.
[12]	Evans, L.C. (2010) Partial Differential Equations. 2nd Edition, American Mathematical Society.
[13]	Masmoudi, M., Pommier, J. and Samet, B. (2005) The Topological Asymptotic Expansion for the Maxwell Equations and Some Applications. Inverse Problems, 21, 547-564. https://doi.org/10.1088/0266-5611/21/2/008
[14]	Schumacher, A. (1995) Topologie optimiering von banteilstrukturen unter verwendung von lopschposition-ierungkieterien. Ph.D. Thesis, Universitat-GesamthochsschuleSiegen.
[15]	Sy, A. (2018) Une méthode d’analyse asymptotique en optimisation topologique. Editions Universitaires Européennes.
[16]	Novotny, A.A. andSokolowski, J. (2013) Topological Derivatives in Shape Optimization. Springer-Verlag.
[17]	Samet, B., Amstutz, S. and Masmoudi, M. (2003) The Topological Asymptotic for the Helmholtz Equation. SIAM Journal on Control and Optimization, 42, 1523-1544. https://doi.org/10.1137/s0363012902406801
[18]	Kozlov, V., Maz’ya, V. and Movchan, A. (1999) Asymptotic Analysis of Fields in Multistructures. Oxford University Press.
[19]	Nazarov, S.A. and Sokołowski, J. (2003) Asymptotic Analysis of Shape Functionals. Journal de Mathématiques Pures et Appliquées, 82, 125-196. https://doi.org/10.1016/s0021-7824(03)00004-7
[20]	Céa, J., Garreau, S., Guillaume, P. and Masmoudi, M. (2000) The Shape and Topological Optimizations Connection. Computer Methods in Applied Mechanics and Engineering, 188, 713-726. https://doi.org/10.1016/s0045-7825(99)00357-6
[21]	Amstutz, S. (2003) Aspects théoriques et numériques en optimisation de forme topologique. Ph.D. Thesis, INSA Toulouse.
[22]	Murat, F. and Simon, J. (1976) Quelques résultats sur le contrôle par un domaine géeométrique. Publ. Labo. D’Anal. Num de l’Université de Paris VI.
[23]	Masmoudi, M. (2002) The Topological Asymptotic. In: Kawarada, H. and Periaux, J., Eds., Computational Methods for Control Applications, Internat. Ser. GAKUTO.
[24]	Hecht, F. (2012) New Development in FreeFem++. Journal of Numerical Mathematics, 20, 251-265. https://doi.org/10.1515/jnum-2012-0013
[25]	Dia, M.B., Faye, I. and Sy, A. (2018) Topological Optimization for Photonic and Phononic Crystals Problems. In: Euclid, P., Lo, G.S., Dia, G., Seydi, H. and Diakhaby, A., Eds., A Collection of Papers in Mathematics and Related Sciences, a festschrift in honour of the late Galaye Dia, SAPS Editions, 523-558. https://doi.org/10.16929/sbs/2018.100-05-03