key: cord-0306092-z8zg3cr3
authors: Nakshatrala, K. B.
title: How to pose material design problems for flow through porous media applications?: Sensitivity of dissipation rate to medium's permeability holds the key
date: 2021-10-20
journal: nan
DOI: 10.1063/5.0076317
sha: aca40a5635803f5a77eda411201b49f409e1ef74
doc_id: 306092
cord_uid: z8zg3cr3

Recent studies have advocated using the total dissipation rate under topology optimization to realize material designs involving the flow of fluids through porous media. However, these studies decided how to pose the design problem, such as maximizing the total dissipation rate for some situations while minimizing for others, by solving one-dimensional problems and justifying their choices using numerical experiments. The rigor is lacking -- a bottleneck for further scientific advancements to computational material design. This paper provides the missing theoretical justification. We identify four classes of boundary value problems using the adjoint state method and analytically calculate the sensitivity of the total dissipation rate to the permeability field. For two of those classes in which the flow of fluids is pressure-driven, the sensitivity is positive -- the total dissipation rate increases if the medium's permeability increases. While for the other two classes, in which the flow is velocity-driven, the trend is the opposite. These sensitivities provide rigorous answers to the central question: how to pose a material design problem for flow through porous media applications. The impact of our work is multi-fold. First, this study further elevates the role of the dissipation rate in posing well-posed material design problems using topology optimization. Second, besides the theoretical significance, the results benefit computational scientists and practitioners to realize optimal designs. Third, given their simplicity yet far-reaching impact, both the approach and results possess immense pedagogical value.

The flow of fluids through porous media is ubiquitous in both natural and synthetic worlds. Examples include nutrient transport in tissues [Dey and Sekhar, 2021; Friedman, 2008] , the flow of hydrocarbons in a petroleum reservoir [Ewing, 1983] , water filtration using membranes [Chen et al., 2021; Herterich et al., 2015] , and workings of facial masks to protect against transmission of COVID-19 and air pollution [Krishan et al., 2021; Verma et al., 2020] , to name a few. Given the field's wide-ranging importance, researchers have been developing mathematical models [Boer, 2012; Coussy, 2004; Rajagopal, 2007] , numerical formulations [Brezzi and Fortin, 1991; Chang et al., 2017; Masud and Hughes, 2002] , computer packages [COMSOL Multiphysics, 2018] , and a posteriori accuracy measures [Shabouei and Nakshatrala, 2016] to study such problems. Material design is always a primary focus given the constant quest to realize new functionalities and optimize associated devices' performance.

However, "optimal" design will have different connotations depending on the application and the functionality one wants to realize. For instance, in water filtration using membranes, the surface of the pores is coated with (chemical) functional groups so that microbes, which have surface charges, are removed via electrostatic attraction or repulsion. The question in such applications is how to manipulate pore surfaces to extract pathogens optimally. In applications involving thermal regulation (e.g., heat exchangers), researchers alter the multi-phase characteristics of the fluids to maximize cooling efficiency. Yet, in some other applications, one designs the very geometrical nature of the pore structure. In this paper, the material design optimizes the dissipation rate via the spatial distribution of porous materials with different permeabilities.

In the past, the limitations of fabrication techniques did not allow researchers to take full advantage of novel designs, which often have complex microstructures. With recent advances in additive manufacturing, the fabrication of complex layouts is now possible [Liu and Ma, 2016; Sigmund, 2009] . So, now that manufacturing is not the bottleneck, researchers look for intelligent designs: realized via designers' ingenuity or mathematical-driven solution procedures.

Topology optimization is a leading mathematical framework for designing material systems [Bendsoe and Sigmund, 2013] . Although the early studies focused on structural optimization, topology optimization is increasingly popular in several other fields involving the flow of fluids through porous media; these emerging fields include the design of microfluidic devices [Alexandersen and Andreasen, 2020; Kreissl et al., 2011] and medical devices such as dialyzer (filter) used for hemodialysis [Alonso and Silva, 2021; Cancilla et al., 2022] . However, many works apply topology optimization in these relatively new application areas by mimicking the material design problem from solid mechanics. For example, some of these works replace the compliance with a similar-looking quadratic functional or replace the principle of minimum potential energy from structural mechanics with the principle of minimum power [Guest and Prévost, 2006; Phatak and Nakshatrala, 2021b] . However, such a translation from solid mechanics to flow through porous media may not be physically meaningful and sometimes not even consistent with the state equations. For example, the minimum power principle does not apply to nonlinear porous models but is limited to linear models such as Darcy and Darcy-Brinkman. More importantly, the said principle is not a physical law but rather a restatement of the governing equations for linear models using the calculus of variations [Phatak and Nakshatrala, 2021b] .

In pursuit of an objective function to drive the material design for flow through porous media applications, Phatak and Nakshatrala [2021b] have promoted the total dissipation rate as an ideal candidate to pose material design problems using topology optimization. This fundamental physical quantity has a universal appeal-applicable even to nonlinear models-which is not the case with the alternatives such as the principle of minimum power. Because of its recent popularity, the total dissipation rate still poses several unanswered questions for defining well-founded material design problems. Three such questions are:

(a) whether should we maximize or minimize the total dissipation rate, (b) should the volumetric bound constraint be on the use of high-or low-permeability material, and (c) when does the optimization problem render nontrivial designs?

The cited paper addressed the first two design questions by considering one-dimensional problems under the Darcy model and supported by numerical results on more complicated boundary value problems. Alluding to the success under the Darcy model, Phatak and Nakshatrala [2021a] have adopted a similar way of posing the material design for the Darcy-Brinkman model. Although the one-dimensional problems offered sound guidance and numerical results provided the proof of concept, a lack of rigor hinders the progress on theoretical and practical fronts.

Thus, the central aim of this paper is to fill this lacuna by providing an in-depth analysis of how to pose the material design problem for flow through porous media applications. We consider Darcy and Darcy-Brinkman models-the two prominent mathematical models describing the flow of incompressible fluids through porous media. Our approach first estimates the design sensitivitythe sensitivity of the total dissipation rate (i.e., the objective function in a material design problem) to the permeability field (i.e., the design variable)-using the adjoint state method 1 .

The adjoint state method has an eventful history, with its first application traced to the optimal control of systems governed by partial differential equations-a seminal work by Lions [1971] . Since then, researchers have applied this method in the design of structural systems [Haug et al., 1986] , seismology [Tromp et al., 2005] , topology optimization [Jensen et al., 2014; Nakshatrala and Tortorelli, 2016; Nakshatrala et al., 2013] , reliability studies and parameter identification [Tortorelli and Michaleris, 1994] , check-pointing in computer codes [Wang et al., 2009] , machine learning [Zhuang et al., 2020] , PDE-constrained optimization [Bradley, 2013; Hinze et al., 2008] , and inverse problems in subsurface modeling [Sun, 2013] . The main advantage of the adjoint state method is in the calculation of the sensitivity of an objective function to a set of design parameters. The method circumvents the explicit calculation of the sensitivity of the solution fields by introducing Lagrange multipliers, calculated by solving an adjoint problem (also referred to as the adjoint state problem), which itself is a system of partial differential equations.

Our analysis-based on the adjoint state method-identifies four principal classes of boundary problems. Two of them are pressure-driven, and the other two are velocity-driven; the adjoint problem is analytically solvable for these classes of problems. As shown in this paper, the answers to the questions raised above-how the total dissipation rate alters with the permeability field and how to pose material design problems-depend on the class of problems. Notably, for the two pressure-driven classes of problems, the total dissipation rate increases with an increase in the permeability field; in comparison, the reverse trend holds for the two velocity-driven classes of problems. Thus, for example, trivial designs occur under different scenarios for pressure-driven and velocity-driven problems. These findings have theoretical, numerical, practical, and pedagogical significance. The research will place the total dissipation rate for material design problems on a sound theoretical footing, achieved by rigorously answering the aforementioned design-related questions.

The plan for the rest of this paper is as follows. Section 2 presents Darcy equations and the dissipation rate under the Darcy model. Section 3 identifies four classes of boundary value problems, which are central to much of the discussion in this paper. Using the adjoint state method, Section 4 derives the sensitivity of the total dissipation to the permeability field and solves the adjoint problem analytically for the four classes of boundary problems. Section 5 shows that these findings translate even to the Darcy-Brinkman model. By synthesizing the results from the previous sections, Section 6 addresses the central and related design questions with mathematical justifications. Section 7 offers concluding remarks.

Consider the flow of an incompressible single-phase fluid in a porous domain Ω with boundary ∂Ω. The fluid's density and coefficient of viscosity are denoted as ρ and µ, respectively. k(x) denotes the permeability of the porous domain at a spatial point x ∈ Ω. The spatial divergence and gradient operators are denoted by div [·] and grad[·], respectively.

A general forward problem 2 under the Darcy model-the so-called Darcy equations-reads: Find velocity vector field v and pressure scalar field p satisfying

where b(x) is the body force vector field; p p (x) and v p n (x) are, respectively, the prescribed pressure and normal component of the velocity on the boundary; n(x) is the outward unit vector on the boundary; and Γ p and Γ v are the parts of boundary on which pressure and velocity boundary conditions are prescribed, respectively. For mathematical well posedness, the two parts of the boundary satisfy:

Solutions to Darcy equations are unique, except for the case when the velocity is prescribed on the entire boundary (i.e., Γ v = ∂Ω). For such a case, a solution exists only when the following compatibility condition is met:

The rationale behind the above compatibility condition is the divergence theorem. To wit,

For the case Γ v = ∂Ω with a compatible velocity boundary condition, the pressure field can be uniquely found only up to a constant; however, the velocity solution field is still unique. The uniqueness of solutions will be handy in ascertaining the solution for adjoint variables.

In a forward problem, one solves Darcy equations to get the solution (i.e., velocity and pressure) fields for a fixed set of inputs, including a prescribed permeability field. For such analyses (with assigned permeability fields), solution fields are mere spatial fields. However, in a material design problem, say topology optimization, the permeability field is the design variable-an unknown by itself. During design iterations, the permeability field changes, altering the solution fields in each iteration. Regarding this paper, we address changes to the total dissipation rate upon changes to the permeability field. To quantify such changes (i.e., to calculate sensitivities), one must account for the implicit dependence of the dissipation rate on the permeability field via the solution fields. Said differently, one cannot treat the solution fields as mere spatial fields, depending only on the spatial location; a change to the permeability field alters the solution fields.

Hence, we distinguish the solution fields from arbitrary spatial fields to emphasize their (implicit) dependence on the permeability field. From hereon, we denote the solutions fields-the Darcy velocity and the associated pressure-as v D x; k(x) and p D x; k(x) . Note the semicolon notation: the solution fields depend explicitly on the quantities on the left side of the semicolon, while the quantity on the right side of the semicolon is the design variable-the permeability field in this paper.

For given permeability k(x) and (arbitrary) vector v(x) fields, the rate of dissipation under the Darcy model takes the following mathematical form:

Given a dissipation rate density, the corresponding total dissipation rate, calculated over the entire domain, reads:

We introduce Φ D [k(x)]-referred to as the total dissipation rate functional 3 -to denote the total dissipation rate expended under the Darcy model by the Darcy velocity field with a permeability field k(x). Mathematically,

2.1. A notation for calculus of functionals. The (Fréchet) derivative of a functional Φ[k(x)] at k(x) is defined as [Spivak, 1971] :

for any perturbation field ∆k(x), and · is the standard Euclidean norm. For a given scalar field of the form p(x; k(x)) (e.g., the pressure field satisfying the forward problem), the derivatives D 1 p and D 2 p are respectively defined as follows:

For convenience, we denote these two derivatives as follows:

) for a vector field v(x; k(x)). As done in tensor algebra, the divergence can then be defined as follows:

where tr[·] is the trace operator for second-order tensors [Chadwick, 1999] .

3 A functional is a function of functions; see [Gelfand and Fomin, 2000 ].

For a given scalar field of the form ϕ x, k(x), v(x; k(x) , we define the derivatives D 2 ϕ and D 3 ϕ as follows:

The definition for D 1 ϕ is similar to the one provided before (cf. Eq. (9)), and grad[ϕ] ≡ D 1 ϕ.

Appealing to the chain rule, we define the total derivative of ϕ with respect to the permeability field as follows:

The first term on the right side of the above equation accounts for explicit dependence, while the second term accounts for implicit dependence.

For further analysis, we identify the following four classes of boundary value problems under the Darcy model: 4

(1) Class A: Prescription of a pressure loading on the entire boundary. Stated mathematically,

(2) Class B: Prescription of pressure loading on a part of the boundary while homogeneous velocity boundary conditions on the rest. That is,

(3) Class C: Prescription of a compatible velocity boundary condition on the entire boundary and requiring the prescribed body force to be a conservative vector field. Mathematically,

where ψ(x) is a scalar field. (4) Class D: Prescription of velocity boundary conditions on a part of the boundary with zero pressure loading on the rest, and requiring the body force is zero. Written mathematically,

The flow of fluids under the first two classes of boundary value problems is pressure-driven, while the flow is velocity-driven for the last two. Even a few other types of boundary value problems, which in their original form do not fall into any of these classes, can be recast into one of the above classes. For example, consider a boundary value problem with a non-zero but constant pressure applied on a part of the boundary; one can transform this problem-by shifting the datum of the pressure to the prescribed (constant) pressure-so that it belongs to Class D. As shown in the next section, design sensitivities can be calculated analytically for the four classes of boundary value problems.

We now estimate the sensitivity of the total dissipation rate to the permeability field. Mathematically, the task at hand is to find DΦ D [k(x) ]. For this estimate, we use the adjoint state method, which offers a systematic procedure to calculate sensitivities of an objective function to a set of design variables.

The objective function, in general, depends explicitly on a set of design parameters and implicitly on the solution fields, which satisfy the forward problem. In addition, the definition of the forward problem involves the design parameters. Thus, the sensitivity of the objective function depends on the sensitivity of the solution fields. The adjoint state method circumvents an explicit calculation of the said sensitivity of the solution fields by introducing Lagrange multipliers that are the solutions of an adjoint problem-a system of partial differential equations.

To crystallize the concept outlined above, we will relate the above description to our problem. As mentioned earlier, our task is to find the sensitivity of the total dissipation rate to the permeability field. The dissipation rate is a function of the permeability field and the velocity of the fluid. The fluid's velocity is a solution of, for example, Darcy equations-a system of partial differential equations. Darcy equations, in turn, are written in terms of permeability. Thus, the objective function is the total dissipation rate, the design parameter is the permeability, and the state equations are Darcy equations. In the rest of this section, we derive the associated adjoint problem and calculate the desired sensitivities.

4.1. Deriving the adjoint problem. The sensitivity of the total dissipation rate functional with respect to its argument-the permeability field-takes the following mathematical form: (21) Following the adjoint state method, we augment the above integral with additional terms, which are essentially zero. These new terms are in the form of integrals constructed from the derivative of the state equations with respect to the design variable (i.e., k(x)). After augmenting these terms, the sensitivity can be written as follows:

where Λ(x) and λ(x) are the newly introduced adjoint variables. The terms in the curly brackets are the residuals of the state equations. The factor 2 is introduced in the last four integrals for convenience, as it simplifies the adjoint problem and the expressions for the associated solution.

Noting that the coefficient of viscosity µ, the outward unit normal vector n(x), and the prescribed fields-ρb(x), v p n (x) and p p (x)-are all independent of the permeability field, we get the following:

By using the Green's identity on the third term of the third integral and on the fourth integral, and grouping the terms we get the following:

In the above equation, we do not know a priori the sensitivities of the solution field variables (i.e., v ′ D (x; k(x)) and p ′ D (x; k(x))). Until now the choice for the adjoint field variables is arbitrary; however, a judicious selection will circumvent the need to find the sensitivities of the solution field, and simplify the expression for design sensitivities. Specifically, the last four integral terms in the above equation can be eliminated if the adjoint field variables satisfy the following boundary value problem:

Under this specific choices for the adjoint variables, the sensitivity of the objective function with respect to the permeability field reduces to:

The above boundary value problem, in terms of the adjoint variables, is commonly referred as the adjoint problem. The above set of equations governing the adjoint variables is similar to Darcy equations with a pseudo body force µv D (x; k(x))/k(x) and homogeneous boundary conditions. In general, one has to use a numerical method to solve the adjoint problem, similar to the case with the forward problem. However, as shown below, for the four classes of boundary value problems defined earlier, we can find analytical solutions for the adjoint variables, and hence express the design sensitivities in terms of the solution fields.

Solutions to the adjoint problem and sensitivities for the four classes. Our task is to find the adjoint variables-Λ(x) and λ(x)-that satisfy equations (25a)-(25d).

is a plausible solution for the adjoint problem. Clearly, the above fields satisfy Eq. (25a). Λ(x) is divergence-free, satisfying Eq. (25b), as the Darcy velocity is divergence free in keeping with the continuity equation of the forward problem (1b). For Class A problems, Γ p = ∅, and hence Eq. (25d) is trivially satisfied. For Class B problems, v p n (x) = 0 on Γ v , implying that Accordingly, the sensitivity of the total dissipation rate to the permeability field, given by Eq. (26), for the Classes A and B problems can be simplified as follows:

Since µ > 0, we have

for any nonzero Darcy velocity field. 4.2.2. Class C boundary value problems. The solution for the adjoint variables will be:

where ψ(x) is the scalar potential for the body force. Recall that the pressure field can be determined up to an arbitrary constant when velocity boundary conditions are prescribed on the entire boundary-which is the case with Class C boundary value problems. The above fields trivially satisfy Eqs. (25b) and (25d). Since Γ p = ∅, Eq. (25c) is also met. Upon substituting the above fields in Eq. (25a), we get:

Noting that ρb(x) = −grad[ψ(x)], the above equation can be rewritten as:

which is equivalent to the balance of linear momentum equation (1a) under the forward problem. By definition, the Darcy velocity and the associated pressure satisfy Eq. (1a). Hence, fields given in Eq. (30) are the solution of the adjoint problem for Class C boundary value problems, appealing again to the uniqueness theorem for Darcy equations. Thus, for Class C boundary value problems, we have

The positivity of the coefficient of viscosity implies that

The above inequality becomes strictly negative for any nonzero Darcy velocity field. 4.2.3. Class D boundary value problems. The solution for the adjoint variables will be:

Clearly, Λ(x) = 0 satisfies Eqs. (25b) and (25d). Since for this class of problems, p D (x; k(x)) = p p (x) = 0 on Γ p , the above solution satisfies Eq. (25c). The above fields make the left side of Eq. (54a) to be −grad[p D (x; k(x))], which is equal to

using Eq. (1a) and noting that the body force is zero for this class of problems, satisfying Eq. (25a). Thus, the fields in Eq. (35) are a solution of the adjoint problem; by appealing to the uniqueness theorem, these fields are the only solution for Class D boundary value problems. Similar to the previous class, the sensitivity of the total dissipation rate to the permeability is:

Since µ > 0, for any nonzero Darcy velocity field we have Said differently, what we have established thus far is the total dissipation rate increases with an increase in permeability for the boundary value problems under Classes A and B; in contrast, the reverse trend holds for boundary value problems under Classes C and D.

The Darcy-Brinkman model-a modification to the Darcy model proposed by Brinkman [1949] accounts for internal dissipation within the fluid besides the drag at the liquid and solid interface, considered in the Darcy model. We now show that the results established for the Darcy model extend to the Darcy-Brinkman model.

We denote the symmetric part of gradient of a vector field w(x) as follows:

where the superscript 'T' denotes the transpose of a second-order tensor. We also define the following projection tensor:

where I is the second-order tensor. It is easy to check that P T = P , P P = P and P n(x) = 0

The second equation is the definition for a tensor to be a projection. The third equation implies that the normal vector n(x) is not in the range space of the projection tensor P ; said differently, n(x) is in the null space of P . Using this projection tensor, any vector field w(x) can be uniquely decomposed into:

Using the above notation, Darcy-Brinkman equations can be written as follows:

where p p (x) is the prescribed pressure, v p (x) is the prescribed tangential component of the velocity vector on Γ p , and v p (x) is the prescribed (full) velocity vector on Γ v . If velocity boundary conditions are prescribed on the entire boundary (i.e., Γ v = ∂Ω), similar to Darcy equations, solutions to Darcy-Brinkman equations exist only if the following compatibility condition is satisfied:

The solutions to Darcy-Brinkman equations are unique except when Γ v = ∂Ω, where, upon meeting the above compatibility condition (44), the pressure is unique only up to an arbitrary constant. See [Shabouei and Nakshatrala, 2016] for a proof of the unique theorem for Darcy-Brinkman equations. From hereon, the solutions fields-the Darcy-Brinkman velocity and the associated pressure-are denoted as v DB (x; k(x)) and p DB (x; k(x)), respectively.

The dissipation rate density for an arbitrary vector field v(x) under the Darcy-Brinkman model takes the following mathematical form:

The total dissipation rate functional-to denote the total dissipation rate expended by the Darcy-Brinkman velocity under the Darcy-Brinkman model-is defined as:

5.1. Four classes of boundary value problems under the Darcy-Brinkman model. Since the boundary conditions under the Darcy-Brinkman model differ from that of the Darcy model, we need to slightly modify the definitions of the classes of the boundary value problems. Still, the flow of fluids under Classes A and B will be pressure driven while Classes C and D consider velocity-driven flows.

(1) Class A: Prescription of a pressure loading on the entire boundary with the corresponding tangential velocity to be zero. Stated mathematically,

(2) Class B: Prescription of pressure loading on a part of the boundary while homogeneous velocity boundary conditions on the rest. That is,

(3) Class C: Prescription of a compatible velocity boundary condition on the entire boundary and requiring the prescribed body force to be a conservative vector field. Mathematically,

where ψ(x) is a scalar field. (4) Class D: Prescription of velocity boundary conditions on a portion of the boundary and homogeneous (pressure and tangential velocity) boundary conditions on the rest, and requiring zero body force. Written mathematically,

5.2. Sensitivities using the adjoint state method. We now calculate the sensitivity of the total dissipation rate to the permeability field, DΦ DB [k(x)], which takes the mathematical form:

Following the adjoint state method, we augment the above term with the derivative of the residuals of Darcy-Brinkman equations:

The rationale behind the incorporation of the factor 2 in the last four terms and the selection of the weights within the integrals, especially in the last two terms, is to simplify the resulting adjoint problem. We note that the coefficient of viscosity µ, the outward unit normal vector n(x), and the prescribed fields-ρb(x), v p (x), v p (x) and p p (x)-are all independent of the permeability field. Using Green's identity and the properties of projection tensor given in Eq. (41), performing algebraic simplifications, and grouping the terms, we get the following:

Following the same approach as in the case of the Darcy model, we choose the adjoint variables to make the last five integrals vanish. Accordingly, Λ(x) and λ(x) satisfy the following adjoint problem, corresponding to Darcy-Brinkman equations:

The above adjoint problem resembles Darcy-Brinkman equations (43a)-(43e) but with homogeneous velocity boundary conditions, a non-homogeneous body force of

and a pressure loading of

The negative sign in the above expression indicates that a pressure loading acts inwards, opposing the unit outward normal vector n(x).

For the adjoint variables satisfying the adjoint problem (54a)-(54e), the sensitivity of the total dissipation rate with respect to the permeability field is:

The above expression is similar to that of the one obtained under the Darcy model (cf. Eq. (26)).

To find the sensitivity, we need to know Λ(x), which is a solution of the adjoint problem. Note that v DB (x; k(x)) is known from the forward problem. 5.2.1. Classes A and B boundary value problems. Even under the Darcy-Brinkman model, the solution for the adjoint variables for these two classes of boundary value problem will be:

Following the same reasoning as done for Darcy equations and appealing to the uniqueness theorem for Darcy-Brinkman equations, the above pair is the only solution of the adjoint problem (54a)-(54e). Consequently,

Since µ > 0, we have

for any nonzero Darcy-Brinkman velocity field. 5.2.2. Class C boundary value problems. The solution for the adjoint variables will be:

For this class of boundary value problems, the sensitivity of the total dissipation rate to the permeability field can be compactly written as follows:

The coefficient of viscosity is positive, implying

for any nonzero Darcy-Brinkman velocity field.

The solution for the adjoint variables will be:

The sensitivity of the total dissipation rate to the permeability field is again:

For any nonzero Darcy-Brinkman velocity field, we have

as µ > 0.

5.3. Discussion. Although Darcy-Brinkman equations consider an additional mode of dissipationinternal friction in the fluid, the trend for the sensitivity of the total dissipation rate to the permeability field remains the same as that of Darcy equations. Moreover, expressions for the sensitivities are similar except for the apparent replacement of v D (x; k(x)) with v DB (x; k(x)). Thus, what we have established is, for the boundary value problems under Classes A and B, the total dissipation rate increases with an increase in the permeability. On the other hand, the opposite is true for boundary value problems under Classes C and D. Mathematically, the trend can be summarized as follows:

Classes C and D:

In topology optimization, the design variable is selection of the material. In porous media applications, the material property is the permeability. The extremization can be minimization or maximization. The questions is should we minimize or maximize, and should we place a volumetric bound on the usage of the high-permeability or low-permeability material. The actual choice depends on the nature of the boundary value problem (i.e., on the forward problem).

In a recently published paper, this question has been addressed by solving one-dimensional problems. However, no rigorous justification is provided, and we now provide the much needed theoretical justification. In general, four combinations of material design problems are possible: maximize/minimize the total dissipation rate with a volume bound constraint on the use of highor low -permeability material.

For Classes A and B, the total dissipation rate increases with the permeability field. For such problems, if the total dissipation rate is maximized with a volume bound constraint only on the low-permeability material-no bound on the use of the high-permeability material-then the design is trivial: the optimal design is the one with the high-permeability material everywhere. On the other hand, for the same classes of problems, if the volume bound constraint is placed on the use of the high-permeability material, the optimal design is nontrivial. If the total dissipation rate is minimized, then volume bound constraint on the high-permeability material produces a trivial design-the low-permeability material is placed everywhere in the domain, while a volume bound constraint is placed on the use of low-permeability material then the design is nontrivial. Table 1 summarizes all these said scenarios-valid for both Darcy and Darcy-Brinkman models. 

This paper addressed how to pose well-posed material design problems-using the rate of dissipation-for flow through porous media applications. The answer to this design question depends on specifics of the boundary value problem-the so-called forward/direct problem. Mathematical analysis using the adjoint state method identified four classes of boundary value problems: two of them involve pressure-driven flows, and the other two are velocity-driven. Using these four classes, this paper rigorously answered:

(1) whether to maximize or minimize the total dissipation rate, (2) should we place a bound constraint on high-or low-permeability material, and (3) conditions for nontrivial designs.

This newfound knowledge puts the material design problem using the dissipation rate under a sound footing.

While answering the design question, this study revealed a fundamental property of the flow of fluids through porous media: how the total dissipation rate varies with the permeability field for the four classes of boundary value problems. The said sensitivity is positive for the two classes with prescribed pressure loading driving the flow. However, the trend is opposite for the other two classes in which velocity boundary conditions drive the flow. Yet another use of the analytical sensitivities derived in this paper is that a computer implementation can use them without solving the adjoint problem numerically, decreasing the time to compute. These four classes of boundary value problems are fundamental, with far-reaching consequences beyond posing material design problems.

(1) Class A: Prescription of a traction on the entire boundary. Stated mathematically,

(2) Class B: Prescription of a traction on a part of the boundary while homogeneous velocity boundary conditions on the rest. That is,

(3) Class C: Prescription of a compatible velocity boundary condition on the entire boundary and requiring the prescribed body force to be a conservative vector field. Mathematically,

where ψ(x) is a scalar field. (4) Class D: Prescription of velocity boundary conditions on a part of the boundary with zero traction on the rest, and requiring the body force is zero. Written mathematically, Γ p = ∅, Γ v = 0, t p (x) = 0 on Γ p , and b(x) = 0 (A8)

A review of topology optimisation for fluid-based problems

Topology optimization for blood flow considering a hemolysis model. Structural and Multidisciplinary Optimization

Topology Optimization: Theory, Methods, and Applications

Theory of Porous Media: Highlights in Historical Development and Current State

PDE-constrained optimization and the adjoint method

Mixed and Hybrid Finite Element Methods

A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles. Flow, Turbulence and Combustion

A porous media CFD model for the simulation of hemodialysis in hollow fiber membrane modules

Continuum Mechanics: Concise Theory and Problems

Modification to Darcy-Forchheimer model due to pressuredependent viscosity: consequences and numerical solutions

Flow and fouling in elastic membrane filters with hierarchical branching pore morphology

Comsol User's Guide, Version 5.3. COMSOL AB

Mathematical modeling of electrokinetic transport through endothelial-cell glycocalyx

The Mathematics of Reservoir Simulation

Principles and Models of Biological Transport

Topology optimization of creeping fluid flows using a Darcy-Stokes finite element

Design Sensitivity Analysis of Structural Systems

Tailoring wall permeabilities for enhanced filtration

Optimization with PDE Constraints

On the consistency of adjoint sensitivity analysis for structural optimization of linear dynamic problems. Structural and Multidisciplinary Optimization

Topology optimization for unsteady flow

Efficacy of homemade face masks against human coughs: Insights on penetration, atomization, and aerosolization of cough droplets

Optimal Control of Systems Governed by Partial Differential Equations

A survey of manufacturing oriented topology optimization methods

A stabilized mixed finite element method for Darcy flow

Nonlinear structural design using multiscale topology optimization. Part II: Transient formulation

Nonlinear structural design using multiscale topology optimization. Part I: Static formulation

Effect of viscous shearing stresses on optimal material designs for flow of fluids through porous media

On optimal designs using topology optimization for flow through porous media applications

On a hierarchy of approximate models for flows of incompressible fluids through porous solids

Mechanics-based solution verification for porous media models

Manufacturing tolerant topology optimization

Calculus on Manifolds: A Modern Approach to Classical Theorems of Advanced Calculus

Inverse Problems in Groundwater Modeling

Design sensitivity analysis: overview and review

Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels

Visualizing the effectiveness of face masks in obstructing respiratory jets

Minimal repetition dynamic checkpointing algorithm for unsteady adjoint calculation

Adaptive checkpoint adjoint method for gradient estimation in neural ODE

The results-solutions to the adjoint problem and expressions for the sensitivities-and conclusions remain the same as the original form of Darcy-Brinkman equations. Since the derivations are a straightforward extension of the original case, we omitted the details for brevity.

An alternative form of Darcy-Brinkman equations, found in the literature, has the following mathematical form:where t p (x) is the prescribed traction, which need not be in the form of a pressure. Here Γ p is the part of the boundary on which traction boundary condition is prescribed. The solution to above forward problem is again denoted by v DB (x; k(x)) and p DB (x; k(x)).Using the adjoint method, the sensitivity of the total dissipation rate to the permeability field can be written as follows:By following a similar procedure as before, we get the following adjoint problem corresponding to the alternative form of Darcy-Brinkman equations:Even under this alternative form, the expression for the sensitivity of the total dissipation rate with respect to the permeability field is same as before:where Λ(x) satisfies the adjoint problem (A3a)-(A3d). Due to a change in one of the boundary conditions, the definitions for the classes of boundary value problems should be modified accordingly for the alternative form of Darcy-Brinkman equations.