Introduction

In the last three decades, various topology optimization (TO) methods have been presented, and most have meanwhile attained a mature state. In addition, their popularity as design tools for achieving solutions to a wide variety of problems involving single/multi-physics is growing consistently. Among these, design problems involving fluidic pressure loads¹ pose several unique challenges, e.g., (i) identifying the structural boundary to apply such loads, (ii) determining the relationship between the pressure loads and the design variables, i.e., defining a design-dependent and continuous pressure field, and (iii) efficient calculation of the pressure load sensitivities. Such problems can be encountered in various applications such as air-, water- and/or snow-loaded civil and mechanical structures (aircraft, pumps, pressure containers, ships, turbomachinery), pneumatically or hydraulically actuated soft robotics or compliant mechanisms and pressure loaded mechanical metamaterials, e.g. , to name a few. Note, the shape or topology and performance of the optimized structures or compliant mechanisms are directly related to the magnitude, location, and direction of the pressure loads which vary with the design. In this paper, a novel approach addressing the aforementioned challenges to optimize and design pressure loaded structures and mechanisms is presented. Hereby we target a density-based TO framework.

In line with the outlined applications, we are not only interested in optimizing pressure-loaded stiff structures, but also in generating pressure-actuated compliant mechanisms (CMs). CMs are monolithic continua which transfer or transform energy, force or motion into desired work. Their performance relies on the motion obtained from the deformation of their flexible branches. The use of such mechanisms is on the rise in various applications as these mechanisms provide many advantages over their rigid-body counterparts. In addition, for a given input actuation, the output characteristic of a compliant mechanism can be customized, for instance, to achieve either output displacement in a certain desired fashion, e.g., path generation , shape morphing or maximum/minimum resulting (contact) force wherein grasping of an object is desired . and provide/mention various TO methods to synthesize structures and compliant mechanisms for the applications wherein input loads and constraints are considered invariant during the optimization. However, as mentioned above, a wide range of different applications with pressure loads can be found. A schematic diagram for a general problem with pressure loads is depicted in Fig. 1, whereas Fig. 2 is used to represent a schematic solution to the design problem with different optimized regions. A key problem characteristic is that the pressure-loaded surface is not defined a priori, but that it can be modified by the optimization process (Fig. 2) to maximize actuation or stiffness. Below, we review the proposed TO methods that involve pressure-loaded boundaries, for either structures or mechanism designs.

were first to present a TO method involving pressure loads. Thereafter, several approaches have been proposed to apply and provide a proper treatment of such loads in TO settings, which can be broadly classified into: (i) methods using boundary identification schemes , (ii) level set method based approaches , and (iii) approaches involving special methods, i.e. which avoid detecting the loading surface .

Boundary identification techniques, in general, are based on a priori chosen threshold density $`\rho_T`$, i.e., iso-density curves/surfaces are identified. used the iso-density approach to identify the pressure loading facets $`\mathrm{\Gamma}_\mathrm{p_b}`$ (Fig. 2) which they further interpolated via Bézier spline curves to apply the pressure loading. However, as per this iso-density (isolines) method may furnish isoline-islands and/or separated isolines. Consequently, valid loading facets may not be achieved. In addition, this method requires predefined starting and ending points for $`\mathrm{\Gamma_{p_b}}`$ . proposed a modified isolines technique to circumvent abnormalities associated with the isolines method. Refs. evaluated the sensitivities of the pressure load with respect to design variables using an efficient finite difference formulation. presented a method wherein one does not need to define starting and ending points a priori. In addition, they provided an analytical approach to calculate load sensitivities. Moreover, these studies considered sensitivities of the pressure loads, however they are confined to only those elements which are exposed to the pressure boundary loads $`\mathrm{\mathrm{\Gamma_{p_b}}}`$.

proposed a method wherein the evolving pressure loading boundary $`\mathrm{\Gamma_{p_b}}`$ is predefined using an additional set of variables, which are also optimized along with the design variables. proposed an element-based search method to locate the load surface. They used the actual boundary of the finite elements (FEs) to construct the load surface and thereafter, transferred pressure to corresponding element nodes directly. introduced an algorithm based on digital image processing and regional contour tracking to generate an appropriate pressure loading surface. They transferred pressure directly to nodes of the FEs. The methods presented in this paragraph do not account for load sensitivities within their TO setting.

As per , if the evolving pressure loaded boundary $`\mathrm{\Gamma_{p_b}}`$ coincides with the edges of the FEs then the load sensitivities with respect to design variables vanish or can be disregarded. Consequently, $`\mathrm{\Gamma_{p_b}}`$ no longer remains sensitive to infinitesimal alterations in the design variables (density fields) unless the threshold value $`\rho_T`$ is passed and thus, $`\mathrm{\Gamma_{p_b}}`$ jumps directly to the edges of a next set of FEs in the following TO iteration. Note that load sensitivities however may critically affect the optimal material layout of a given design problem, especially those pertaining to compliant mechanisms, as we will show in Sec. 14.5. Therefore, considering load sensitivities in problems involving pressure loads is highly desirable. In addition, ideally these sensitivities should be straightforward to compute, implement and computationally inexpensive.

In contrast to density-based TO, in level-set-based approaches an implicit boundary description is available that can be used to define the pressure load. On the other hand, being based on boundary motion, level-set methods tend to be more dependent on the initial design . employed a level set function (LSF) to represent the structural topology and overcame difficulties associated with the description of boundary curves in an efficient and robust way. employed two zero-level sets of two LSFs to represent the free boundary and the pressure boundary separately. employed the Distance Regularized Level Set Evolution (DRLSE) to locate the structural boundary. They used the zero level contour of an LSF to represent the loading boundary but did not regard load sensitivities. Recently, proposed a method wherein Laplace’s equation is employed to compute hydrostatic fluid pressure fields, in combination with interface tracking based on a flood fill procedure. Shape sensitivities in conjunction with Ersatz material interpolation approach are used within their approach.

Given the difficulties of identifying a discrete boundary within density-based TO and obtaining consistent sensitivity information, various researchers have employed special/alternative methods (without identifying pressure loading surfaces directly) to design structures experiencing pressure loading. presented an approach based on applying a fictitious thermal loading to solve pressure loaded problems. employed a mixed displacement-pressure formulation based finite element method in association with three-phase material (fluid/void/solid). Therein, an extra (compressible) void phase is introduced in the given design problem while limiting the volume fraction of the fluid phase and also, the mixed finite element methods have to fulfill the BB-condition which guarantees the stability of the element formulation . also used three-phase material to solve such problems. introduced a pseudo electric potential to model evolving structural boundaries. In their approach, pressure loads were directly applied upon the edges of FEs and thus, they did not account for load sensitivities. Additional physical fields or phases are typically introduced in these methods to handle the pressure loading. Our method follows a similar strategy based on Darcy’s law, which has not been reported before.

This paper presents a new approach to design both structures and compliant mechanisms loaded by design-dependent pressure loads using density-based topology optimization. The presented approach uses Darcy’s law in conjunction with a drainage term (Sec. 5.1.1) and standard FEs, for modeling and providing a suitable treatment of pressure loads. The drainage term is necessary to prevent pressure loads on structural boundaries that are not in contact with the pressure source, as explained in Sec. 5.1.1. Darcy’s law is adapted herein in a manner that the porosity of the FEs can be taken as design (density) dependent (Sec. 5.1) using a smooth Heaviside function facilitating smoothness and differentiability. Consequently, prescribed pressure loads are transferred into a design dependent pressure field using a PDE (Sec. 5.2.1) which is further solved using the finite element method. The determined pressure field is used to evaluate consistent nodal forces using the FE method (Sec. 5.2.2). This two step process offers a flexible and tunable method to apply the pressure loads and also, provides distributed load sensitivities, especially in the early stage of optimization. The latter is expected to enhance the exploratory characteristics of the TO process.

In addition, regarding applications most research on topology optimization involving pressure loads has thus far focused on compliance minimization problems and, a thorough search yielded only two research articles for designing pressure-actuated compliant mechanisms. employed the three-phase method proposed in to generate such mechanisms actuated via pressure loads whereas also used the three-phase method but in association with a displacement-based nonconforming FE method, which is not a standard FE approach. Herein, using the presented method, we not only design pressure-loaded structures but also pressure-actuated compliant mechanisms, which suggests the novel potentiality of the method.

In summary, we present the following new aspects:

• Darcy’s law is used with a drainage term to identify evolving pressure loading boundary which is performed by solving an associated PDE,

• the approach facilitates computationally inexpensive evaluation of the load sensitivities with respect to design variables using the adjoint-variable method,

• the load sensitivities are derived analytically and consistently considered within the presented approach while synthesizing structures and compliant mechanisms experiencing pressure loading,

• the importance of load sensitivity contributions, especially in the case of compliant mechanisms, is demonstrated,

• the method avoids explicit description of the pressure loading boundary (which proves cumbersome to extend to 3D),

• the robustness and efficacy of the approach is demonstrated via various standard design problems related to structures and compliant mechanisms,

• the method employs standard linear FEs, without the need for special FE formulations.

The remainder of the paper is organized as follows: Sec. 5 describes the modeling of pressure loading via Darcy’s law with a drainage term. Evaluation of consistent nodal forces from the obtained pressure field is presented therein. In Sec. 11, the topology optimization problem formulation for pressure loaded structures and small-deformation compliant mechanisms is presented with the associated sensitivity analysis. In addition, the presented method is verified using a pressure-loaded structure problem on a coarse mesh. Sec. 14 presents the solution of various benchmark design problems involving pressure loaded structures and small deformation compliant mechanisms. Lastly, conclusions are drawn in Sec. 8.

Modeling of Design Dependent Loading

The material boundary of a given design domain $`\mathrm{\Omega}`$ evolves as the TO progresses while forming an optimum material layout. Therefore, it is challenging especially in the initial stage of the optimization to locate an appropriate loading boundary $`\mathrm{\Gamma}_\mathrm{p_b}`$ for applying the pressure loads. In addition, while designing especially pressure-actuated compliant mechanisms, establishing a design dependent and continuous pressure field would aid to TO. Herein, Darcy’s law in conjunction with the drainage term, a volumetric material-dependent pressure loss, is employed to establish the pressure field as a function of material density vector $`\bm{\rho}`$.

Darcy’s law

Darcy’s law defines the ability of a fluid to flow through porous media such as rock, soil or sandstone. It states that fluid flow through a unit area is directly proportional to the pressure drop per unit length $`\nabla p`$ and inversely proportional to the resistance of the porous medium to the flow $`\mu`$ . Mathematically,

\begin{equation}
\label{sec2:eq1}
 \bm{q} = -\frac{\kappa}{\mu}\;\nabla p \quad = -K \;\nabla p,
\end{equation}

where $`\bm{q}, \,\kappa, \,\mu,\,\text{and},\,\nabla p`$ represent the flux ($`\si{\meter\per\second}`$), permeability ($`\si{\meter\squared}`$), fluid viscosity ($`\si{\newton\per\square\meter\second}`$) and pressure gradient ($`\si{\newton\per\cubic\meter}`$), respectively. Further, $`K`$ ($`\si{\meter\tothe{4}\per\newton\per\second}`$) is termed herein as a flow coefficient² which expresses the ability of a fluid to flow through a porous medium. The flow coefficient of each FE is assumed to be related to element density $`\rho_e`$. In order to differentiate between void ($`\rho_e=0`$) and solid ($`\rho_e=1`$) states of a FE, and at the same time ensuring a smooth and differentiable transition, $`K(\rho_e)`$ is modeled using a smooth Heaviside function as:

\begin{equation}
\label{sec2:eq2}
 K(\rho_e) = k_\mathrm{v} - k_\mathrm{vs} \frac{\tanh{\left(\beta_\mathrm{k}\eta_\mathrm{k}\right)}+\tanh{\left(\beta_\mathrm{k}(\rho_e - \eta_\mathrm{k})\right)}}{\tanh{\left(\beta_\mathrm{k} \eta_\mathrm{k}\right)}+\tanh{\left(\beta_\mathrm{k}(1 - \eta_\mathrm{k})\right)}},
\end{equation}

where $`k_\mathrm{vs}= (k_\mathrm{v}-k_\mathrm{s})`$, $`k_\mathrm{v}`$ and $`k_\mathrm{s}`$ are the flow coefficients for a void and solid FE, respectively. Further, $`\eta_\mathrm{k}`$ and $`\beta_\mathrm{k}`$ are two adjustable parameters which control the position of the step and the slope, respectively (Fig. 4). For sufficiently high $`\beta_\mathrm{k}`$, when $`\eta_\mathrm{k}>\rho_e`$, $`K(\rho_e) = k_\mathrm{v}`$ while when $`\eta_\mathrm{k}<\rho_e`$, $`K(\rho_e) = k_\mathrm{s}`$. In view of the permeability of an impervious material and viscosity of air, the flow coefficient of a solid element is chosen to be $`k_\mathrm{s} = \SI{e-10}{\meter\tothe{4}\per\newton\per\second}`$, whereas, $`k_\mathrm{v} = \SI{e-3}{\meter\tothe{4}\per\newton\per\second}`$ is taken to mimic a free flow with low resistance through the void regions.

Our intent is to smoothly and continuously distribute the pressure drop over a certain penetration depth of the solid facing the pressure source. To examine the interaction between structural features and applied pressure under Darcy’s law, consider Fig. 5. Darcy’s law renders a gradual pressure drop from the inner pressure boundary $`\mathrm{\Gamma_{p_b}}`$ to the outer pressure boundary $`\mathrm{\Gamma_{p_0}}`$ (Fig. 5). Consequently, equivalent nodal forces appear within the material as well as upon the associated boundaries. This penetrating pressure, originating because of Darcy’s law, is a smeared-out version of an applied pressure load on a sharp boundary or interface³. Note that, summing up the contributions of penetrating loads gives the resultant load. It is assumed that local differences in the load application have no significant effect on the global behaviour of the structure, in line with the Saint-Venant principle. The validity of this assumption will be checked later in a numerical example (Sec. 11.4).

Drainage term

Application of Darcy’s law alone introduces an undesired pressure distribution in the model when multiple walls are encountered between $`\mathrm{\Gamma_{p_b}} (p_\mathrm{in})`$ and $`\mathrm{\Gamma_{p_0}}(p_\mathrm{out})`$. That is, the pressure does not completely drop over the first boundary as illustrated in Fig. 6. To mitigate this issue, we introduce a drainage term, which is a volumetric density-dependent pressure loss, as

\begin{equation}
\label{sec2:eq3}
 {Q}_\mathrm{drain} = - H(\rho_e) (p - p_{\mathrm{out}}),
\end{equation}

where $`{Q}_\mathrm{drain}`$ denotes volumetric drainage per second in a unit volume ($`\si{\per\second}`$). $`H,\,p,\,p_\mathrm{out}`$ are drainage coefficient ($`\si{\meter\tothe{2}\per\newton\per\second}`$), continuous pressure field ($`\si{\newton\per\square\meter}`$), external pressure⁴ ($`\si{\newton\per\square\meter}`$), respectively. Conceptually, this term should drain/absorb the flow in the exterior structural boundary layer exposed to the pressure source, so that negligible flow (and pressure) acts on interior structural boundaries.

Similar to flow coefficient $`K (\rho_e)`$, the drainage coefficient $`H(\rho_e)`$ is also modeled using a smooth Heaviside function such that pressure drops to zero when $`\rho_e = 1`$ (Fig. 7). It is given by:

\begin{equation}
\label{sec2:eq4}
 H(\rho_e)    = h_{\mathrm{s}} \frac{\tanh{\left(\beta_\mathrm{h} \eta_\mathrm{h}\right)} + \tanh{\left(\beta_\mathrm{h} (\rho_e-\eta_\mathrm{h})\right)}}{\tanh{\left(\beta_\mathrm{h} \eta_\mathrm{h}\right)} + \tanh{\left(\beta_\mathrm{h} (1-\eta_\mathrm{h})\right)}},
\end{equation}

where, $`\beta_\mathrm{h}`$ and $`\eta_\mathrm{h}`$ are adjustable parameters similar to $`\beta_\mathrm{k}`$ and $`\eta_\mathrm{k}`$. $`h_\mathrm{s}`$ is the drainage coefficient of solid, which is used to control the thickness of the pressure-penetration layer. This formulation can effectively control the location and depth of penetration of the applied pressure. Note, $`h_\mathrm{s}`$ is related to $`k_\mathrm{s}`$ (Appendix 12) as:

\begin{equation}
\label{sec2:eq5}
 h_\mathrm{s} =\left(\frac{\ln{r}}{\Delta s}\right)^2 k_\mathrm{s},
\end{equation}

where $`r`$ is the ratio of input pressure at depth $`\Delta`$s, i.e., $`p|_{\Delta s} = rp_\mathrm{in}`$. Further, $`\Delta s`$ is the penetration depth of pressure, which can be set to the width or height of few FEs. Fig. 9 depicts a plot for the drainage coefficient $`H(\rho_e)`$ as a function of density. Note that the Heaviside parameters used in this plot are the same as those employed in Fig. 4.

Finite Element Formulation

This section presents the FE formulation of the proposed pressure load based on Darcy’s law, wherein the approach employs the standard FE method to solve the associated boundary value problems to determine the pressure and displacement fields. Standard 2D quadrilateral elements with bilinear shape functions are employed to parameterize the design domain. First, in addition to the Darcy equation (Eq. [sec2:eq1]), the equation of state using the law of conservation of mass in view of incompressible fluid is derived. Thereafter, the consistent nodal loads are determined from the derived pressure field.

State Equation

Fig. 10 shows in- and outflow through an infinitesimal volume element $`\mathrm{\Omega_e}`$. Now, using the conservation of mass for incompressible fluid one writes:

\begin{equation}
 \label{sec2:eq6}
\begin{aligned}
\left(q_{x}d y\;+\; q_{y}dx \;+\; {Q}_\mathrm{drain}d xd y \right)d z =&\\
\left(q_{x}d y \;+\; q_{y} d x \;+\; \left(\frac{\partial q_x}{\partial x}d x\right) d y \;+\;  \left(\frac{\partial q_y}{\partial y}d y\right) d x\right)d z, \\ \text{or,}\,\,
\frac{\partial q_x}{\partial x} + \frac{\partial q_y}{\partial y}-{Q}_\mathrm{drain}  = &0,\\ \text{or,}\,\,
\nabla\cdot\bm{q} -{Q}_\mathrm{drain} = &0.
\end{aligned}
\end{equation}

where $`q_x`$ and $`q_y`$ are the flux in $`x`$- and $`y`$-directions, respectively. In view of Eq. ([sec2:eq1]), Eq. ([sec2:eq6]) becomes:

\begin{equation}
\label{sec2:eq7}
\nabla\cdot \left(K\nabla p(\bm{x})\right) + {Q}_\mathrm{drain} = 0.
\end{equation}

Now, for the finite element formulation, we use the Galerkin approach to seek an approximate solution $`p (\bm{x})`$ such that:

\begin{equation}
\label{sec2:eq8}
\sum_{e=1}^{n_\mathrm{elem}}\left(\int_{\mathrm{\Omega}_e}\nabla\cdot \left(K\nabla p(\bm{x})\right)w(\bm{x}) d V + \int_{\mathrm{\Omega}_e}{Q}_\mathrm{drain}w(\bm{x}) d V\right) = 0,
\end{equation}

for every $`w (\bm{x})`$ constructed from the same basis functions as those employed for $`p(\bm{x})`$. The total number of elements is indicated via $`n_\mathrm{elem}`$. In the discrete setting, within each $`\mathrm{\Omega}_e{|_{e=1,\,2,\,3,\,\cdots,\,n_\mathrm{elem}}}`$, we have

\begin{equation}
\label{sec2:eq9}
p_e = \mathbf{N}_\text{p}\mathbf{p}_{e}, \qquad w = \mathbf{N}_\text{p} \mathrm{\bm{w}}_{e},
\end{equation}

where $`\mathbf{N}_\text{p} = [N_1,\, N_2,\,N_3,\,N_4]`$ are the bilinear shape functions in a physical element and $`\mathbf{p}_{e} = \tr{[p_1,\,p_2,\,p_3,\,p_4]}`$ is the nodal pressure. Now, with integration by parts and Greens’ theorem, Eq. ([sec2:eq8]) becomes on elemental level:

\begin{equation}
\label{sec2:eq10}
\begin{aligned}
\int_{\mathrm{\Omega}_e} K\left(\nabla w (\bm{x})\right)\cdot \left(\nabla p(\bm{x})\right)d V + \int_{\mathrm{\Omega}_e}{Q}_\mathrm{drain}w(\bm{x}) d V \\= -\int_{\mathrm{\Gamma}_e} w(\bm{x})\bm{q}_\mathrm{\Gamma}.\bm{n}_edA,
\end{aligned}
\end{equation}

where $`\bm{n}_e`$ is the boundary normal on surface $`\mathrm{\Gamma}_e`$ and therein, $`\bm{q}`$ changes to $`\bm{q}_\mathrm{\Gamma}`$. In view of Eq. ([sec2:eq3]) and Eq. ([sec2:eq9]), Eq. ([sec2:eq10]) gives:

\begin{equation}
 \label{sec2:eq11}
\begin{aligned}
\underbrace{\int_{\mathrm{\Omega}_e}\left( K~ \tr{\mathbf{B}}_\mathrm{p} \mathbf{B}_\mathrm{p}   + H ~\tr{\mathbf{N}}_\mathrm{p} \mathbf{N}_\mathrm{p} \right)d V}_{\mathbf{A}_e}~\mathbf{p}_e=\\
\underbrace{\int_{\mathrm{\Omega}_e}~H~\tr{\mathbf{N}}_\mathrm{p} p_\mathrm{out} ~~d V -
    \int_{\mathrm{\Gamma}_e}~ \tr{\mathbf{N}}_\mathrm{p} \bm{q}_\mathrm{\Gamma} \cdot \bm{n}_e~~d A}_{\mathbf{f}_e},
\end{aligned}
\end{equation}

where $`\mathbf{B}_\mathrm{p} =\nabla\mathbf{N}_\mathrm{p}`$ and $`\bm{q}_\mathrm{\Gamma}`$ is the Darcy flux through the boundary $`\mathrm{\Gamma}_e`$. In global sense, i.e., after assembly , Eq. ([sec2:eq11]) is written as

\begin{equation}
\label{sec2:eq12}
\mathbf{A}\mathbf{p} = \mathbf{f},
\end{equation}

where $`\mathbf{A}`$ is termed the global flow matrix, $`\mathbf{p}`$ and $`\mathbf{f}`$ are the global pressure vector and loading vector, respectively. Note, when $`p_\mathrm{out}= 0`$ and $`q_\mathrm{\Gamma}=0`$ then conveniently $`\mathbf{f} = 0`$ and therefore, the right hand side only contains the contribution from the prescribed pressure, which is the case we have considered while solving design problems in this paper.

Pressure field to consistent nodal loads

The force resulting from the pressure field is expressed as an equivalent body force. Fig. 11 depicts an infinitesimal volume element with pressure loads acting on it, which is used to relate the pressure field $`{p}(\bm{x})`$ and body force $`\bm{b}`$.

Writing the force equilibrium equations, one obtains:

\begin{equation}
\label{sec2:eq13}
\begin{bmatrix} 
pd z d y - pd z d y - \left(\pd{p}{x}d x\right)d z d y\\
pd z d x - pd z d x - \left(\pd{p}{y}d y\right)d z d x\\
pd x d y - pd x d y - \left(\pd{p}{z}d z\right)d x d y
\end{bmatrix} = \begin{bmatrix} b_x \\ b_y\\ b_z \end{bmatrix}d V,
\end{equation}

where, $`b_x,\,b_y,\,\text{and}\,b_z`$ are the components of the body force in $`x,\,y,\,\text{and}\,z`$ directions respectively. Eq. ([sec2:eq13]) can be written as⁵:

\begin{equation}
\label{sec2:eq14}
\quad \bm{b} dV = -\nabla p dV.
\end{equation}

In the discretized setting, $`-\nabla p dV = -\mathbf{B}_\mathrm{p} \mathbf{p}_e dV`$. In general, the external elemental force originating from the body force $`\bm{b}`$ and traction $`\bm{t}`$ in a FE setting , can be written as:

\begin{equation}
\label{sec2:eq15}
\mathbf{F}^e = \int_{\mathrm{\Gamma}_e} \tr{\mathbf{N}}_\mathbf{u}  \bm{t} \;dA \;+\int_{\mathrm{\Omega}_e} \tr{\mathbf{N}}_\mathbf{u} \bm{b} \;dV,
\end{equation}

where $`\mathbf{N}_\mathbf{u} = [N_1\mathbf{I},\, N_2\mathbf{I},\,N_3\mathbf{I},\,N_4\mathbf{I}]`$ with $`\mathbf{I}`$ as the identity matrix in $`\mathcal{R}^2`$ herein. In this work, we consider $`\bm{t} = 0`$. Thus, Eq. ([sec2:eq15]) gives the consistent nodal loads on elemental level as:

\begin{equation}
\label{sec2:eq16}
\mathbf{F}^e = - \int_{\mathrm{\Omega}_e} \tr{\mathbf{N}}_\mathbf{u} \nabla p d V = - \underbrace{\int_{\mathrm{\Omega}_e} \tr{\mathbf{N}}_\mathbf{u} \mathbf{B}_\mathrm{p}  d V}_{\mathbf{H}_e} \mathbf{p}_e.
\end{equation}

Next, in the global form, the consistent nodal loads $`\mathbf{F}`$ can be evaluated from the global pressure vector $`\mathbf{p}`$ (Eq. [sec2:eq12]) using the global conversion matrix $`\mathbf{H}`$ obtained by assembling all such $`\mathbf{H}_e`$ as:

\begin{equation}
\label{sec2:eq17}
\mathbf{F} = -\mathbf{H}\mathbf{p}.
\end{equation}

Note that H is independent of the design, the design-dependence of the loading enters through the pressure field obtained through Darcy’s law (Eq. [sec2:eq12]).

Conclusions

In this paper, a novel approach to perform topology optimization of design problems involving both pressure loaded structures and pressure-actuated compliant mechanisms is presented in a density-based setting. The approach permits use of standard finite element formulation and does not require explicit boudary description or tracking.

As pressure loads vary with the shape and location of the exposed structural boundary, a main challenge in such problems is to determine design dependent pressure field and its design sensitivity. In the proposed method, Darcy’s law in conjunction with a drainage term is used to define the design dependent pressure field by solving an associated PDE using the standard finite element method. The porosity of each FE is related to its material density via a smooth Heaviside function to ensure a smooth transition between void and solid elements. The drainage coefficient is also related to material density using a similar Heaviside function. The determined pressure field is further used to find the consistent nodal loads. In the early stage of the optimization, the obtained nodal loads are spread out within the design domain and thus, may enhance exploratory characteristics of the formulation and thereby the ability of the optimization process to find well-performing solutions.

The Darcy’s parameters, selected a priori to the optimization, affect the topologies of the final continua, and recommended values are provided based on the reported numerical experiments. The method facilitates analytical calculation of the load sensitivities with respect to the design variables using the computationally inexpensive adjoint-variable method. This availability of load sensitivities is an important advantage over various earlier approaches to handle pressure loads in topology optimization. In addition, it is noticed that consideration of load sensitivities within the approach does alter the final optimum designs, and that the load sensitivity terms are particularly important when designing compliant mechanisms. Moreover, in contrast to methods that use explicit boundary tracking, the proposed Darcy method offers the potential for relatively straightforward extension to 3D problems.

The effectiveness and robustness of the proposed method is verified by minimizing compliance and multi-criteria objectives for designing pressure-loaded structures and compliant mechanisms, respectively with given resource constraints. The method allows relocation of the pressure-loaded boundary during optimization, and smooth and steady convergence is observed. Extension to 3D structures and large displacement problems are prime directions for future research.

Acknowledgment

The authors are grateful to Krister Svanberg for providing the MATLAB implementation of his Method of Moving Asymptotes, which is used in this work.