Computer Modeling in Engineering & Sciences

DOI: 10.32604/cmes.2021.015894

ARTICLE

A Combined Shape and Topology Optimization Based on Isogeometric Boundary Element Method for 3D Acoustics

Jie Wang, Fuhang Jiang, Wenchang Zhao and Haibo Chen*

CAS Key Laboratory of Mechanical Behavior and Design of Materials, Department of Modern Mechanics, University of Science and Technology of China, Hefei, 230026, China
*Corresponding Author: Haibo Chen. Email: hbchen@ustc.edu.cn
Received: 21 January 2021; Accepted: 18 February 2021

1  Introduction

Isogeometric analysis (IGA) [1] has been proposed to eliminate the gap between computer-aided design (CAD) and computer-aided engineering (CAE), where non-uniform rational basis splines (NURBS) are widely used to represent geometry and approximate physical field variables, due to its advantages such as accurate description of geometry, high-order continuous fields, and flexible refinement. Compared with the finite element method (FEM) [2], boundary element method (BEM) [3] has exhibited superior performance in problems like exterior acoustics [4] because of its advantages including discretizing in boundary and fitting for infinite domain problems. The fully coupled FEM–BEM algorithm [5] was also proposed for acoustic-structural optimization design [6] to exploit the benefits of the two methods. Besides, various accelerated techniques such as adaptive cross approximation (ACA) [7], fast multipole method (FMM) [8–11], and wideband FMM [4,12,13], have also been adopted to improve BEM’s computational efficiency for large-scale problems. Considering that complex structures are usually described by NURBS patch in CAD software, which means that the geometry model are boundary-represented, the isogeometric boundary element method has also been successfully applied to potential problems [14–17], elastic problem [18–20], and acoustics [21–24]. FEM–BEM coupled analysis method was also combined with IGA in elasticity [25]. As the NURBS control points play a critical role in geometry construction, it is highly flexible to change the structural geometric shape. Refinement operators [1] such as h-refinement, p-refinement, and k-refinement have also been proposed to perfect IGA. Thus, conducting optimization design via IGA BEM is convenient.

The problem of sound emission of structures has attracted much attention in engineering practice. Structures such as sound barriers have also been investigated for their acoustic performance [26,27] in noise control. Baulac et al. [28] performed an optimization design of a sound barrier via BEM. Chen et al. [29] studied the distribution of absorbing material on a noise barrier by optimization design. Zhao et al. [30] applied topology optimization to the distribution of sound absorption material (SAM) that covers the structural surface with FMM BEM. Zhao et al. [31] also investigated the topology optimization of two different materials based on FEM–BEM coupled method. Liu et al. [32] and Chen et al. [33] conducted the sound barrier’s shape optimization by IGA BEM in 2D acoustics. Among these approaches, shape optimization and topology optimization are mainly used for designing structures to change their acoustic performance.

Shape optimization plays a great role in engineering problems to obtain the optimal shape under certain constraints with objectives. In this rational and automatic process, the geometry shape needs to be regenerated in each step, which causes limitation in conventional numerical analysis. In IGA, the structure surface is constructed by NURBS control points with weights and knot vectors; thus, the surface is easy to change by defining the control points as design parameters. Although IGA FEM has been firstly combined with shape optimization in fluid mechanics [34] and electromagnetism [35], IGA BEM is more suitable for the surface shape change in shape optimization, due to NURBS’ and BEM’s features. Shape optimization has also been conducted through IGA BEM in heat conduction [36], wave resistance [37], elasticity [38–40], and acoustics [13,32,41]. In engineering problems, sound structures may be covered with SAM or other materials. Different from design structure surface shape, the topology optimization focuses on finding the best topological distribution of SAM on the surface. Chen et al. [33] firstly conducted the topology optimization of the distribution of SAM on the sound barrier surface by IGA BEM with NURBS. Chen et al. [42] also studied the distribution of SAM via IGA BEM with subdivision surface method for 3D acoustics. Their findings demonstrated the frequency-dependent phenomenon in acoustic optimization design.

The previous works mainly conduct one type of optimization method to improve acoustic performance of structures. The present study focuses on combining the shape and topology optimization together efficiently to ensure better noise control for 3D acoustic problems. For combined optimization, researchers have done many works with the level set method (LSM) to conduct shape and topology update simultaneously by a uniform manner [43–45]. Matsumoto et al. [44] adopted the LSM in the acoustic scatter problem to design the shape and topology of the structure. In this work, the topology optimization we considered is the distribution of SAM on structural surface. The solid isotropic material with penalization (SIMP) [46] method is applied to the initialization of design variables by converting the discrete density values of artificial material into continuous variables and enabling the intermediate density to approach 0 or 1, where 0 means that no material is on the design area. Thus, we select another type of combined optimization, where the shape and topology of structure should be changed by constructing a connection between them. Researchers have applied the deformable simplicial complex (DSC) method to combine shape and topology optimization in elasticity problems [47–49], where the topology design variables were updated first, and then the shape update was implemented subsequently. This process was called total iteration step of the combined optimization. Lin et al. [50] performed shape and topology optimizations through a time-dependent adjoint approach for sensitivity analysis in 2D acoustic metamaterials and phononic crystals, where shape optimization followed the topology optimization subsequently to modify the structural geometries. Jiang et al. [51] compared four iteration schemes by the second kind of combined optimization and successfully applied them to the noise reduction on sound barriers based on IGA BEM in 2D acoustics. In this paper, we mainly expand the best iteration scheme proposed by Jiang et al. [51] to 3D acoustics and applied it to complicated structures, where the NURBS control points coordinates are set as design parameters for shape design, and the artificial densities on NURBS integral elements are selected as the design variables for topology analysis. As both the two types of parameters are related to IGA BEM, a bridge is constructed between the structural shape and the distribution of SAM on the surface, which helps the shape and topology change in each iteration step. Meanwhile, using IGA BEM also avoids mesh regeneration, which improves the efficiency in the optimization process.

To implement combined optimization, the gradient-based optimizer is generally used due to its high efficiency, for example, the method of moving asymptotes (MMA) [52]. This condition leads to the requirement of gradient information, which means additional step, i.e., sensitivity analysis. The present study consists of two parts: shape sensitivity analysis with respect to the control points, and topology sensitivity analysis with respect to the artificial densities. Finite difference method (FDM), direct differentiation method (DDM) [4,12,53], and adjoint variable method (AVM) [30,54,55] are the three main methods for both the two types of sensitivity analysis. The DDM and AVM provide the analytical expression of sensitivity, and the AVM shows higher efficiency for multiple design variables. As the shape design parameters and topology design variables may be large in engineering, AVM is more suitable for ordinary optimization problems to achieve high efficiency. Here, we adopt the AVM for both the two types of sensitivity analysis. Another difficulty in IGA BEM for 3D acoustics is the non-uniqueness problem [56,57], especially for exterior acoustics. The Burton–Miller method [56,58] is adopted to overcome this problem, where Marburg [59] has investigated the influence of different coupling parameters on the analysis result. However, using this method also results in the problem of hyper singular integrals [41]. Thus, the method based on subtraction of singularity technique [41,60] is also adopted to compute the singular integrals to achieve high accuracy. Moreover, to improve the efficiency of analysis conveniently, OpenMP tool in Fortran code is exploited in the matrix computation without changing the formulation of IGA BEM.

The present work is devoted to extending the combined optimization algorithm to 3D acoustics. The rest of this paper is organized as follows. In Section 2, the IGA BEM for 3D acoustics with the impedance boundary conditions are reviewed. Section 3 discusses the sensitivity analysis via AVM with respect to the NURBS control points and the artificial densities of SAM, respectively. Section 4 describes the three optimization procedure: shape optimization, topology optimization, and combined optimization. Section 5 introduces several numerical examples to validate the proposed approach. Finally, Section 6 provides the conclusions of this work.

2  Isogeometric BEM for 3D Acoustics

2.1 NURBS Surface for IGA

In IGA, we can define the knot vector Ξ={ξ1,ξ2,…,ξn+pg+1}, where ξi∈ℝ is the i-th knot, n is the number of basis functions, and pg is the polynomial order. Then, B-spline basis functions can be defined as follows:

Ni,0(ξ)={1,ifξi≤ξ<ξi+1,0,otherwise.(1)

For pg > 1, we have the following basis functions:

Ni,pg(ξ)=ξ-ξiξi+pg-ξiNi,pg-1(ξ)+ξi+pg+1-ξξi+pg+1-ξi+1Ni+1,pg-1(ξ).(2)

Then, a B-spline curve’s formulation can be described by

x(ξ)=∑i=1nNi,pg(ξ)Pi,(3)

where x(ξ) is the point located on the curve, and Pi is the i-th control point that corresponds to the knot.

For B-spline surface, we need two knot vectors Ξ1={ξ1,ξ2,…,ξn+pg1+1} and Ξ2={η1,η2,…,ηm+pg2+1}, where n and m are the numbers of basis functions for each parametric dimension. The formulation of B-spline surfaces can be described as follows:

x(ξ,η)=∑i=1n∑j=1mNi,pg1(ξ)Nj,pg2(η)Pi,j.(4)

By introducing a weight wi,j with each control point, NURBS basis functions for two dimensions can be obtained as

Ri,j(ξ,η)=Ni,pg1(ξ)Nj,pg2(η)wi,j∑a=1n∑b=1mNa,pg1(ξ)Nb,pg2(η)wa,b.(5)

Similarly, the NURBS surface formulation can be described as follows:

x(ξ,η)=∑i=1n∑j=1mRi,j(ξ,η)Pi,j=∑l=1NgRl(ξ,η)Pl,(6)

where the number of control points is represented by Ng=n×m. The derivative of the NURBS basis function was discussed by Simpson et al. [18].

2.2 Isogeometric BEM for Acoustics

Considering a domain Ω with a boundary Γ, the Helmholtz equation that governs the acoustic field can be described as

∇2p(x)+k2p(x)=0,∀x∈Ω,(7)

where p denotes the sound pressure, k=ω/c denotes the wave number, ω is the angular frequency, and c is the wave velocity in the acoustic medium. Here, we select e-iωt as the harmonic time dependence, where t denotes time.

For boundary element method, by applying Green’s second theorem and letting point x approach the boundary, we can obtain the conventional boundary integral equation (BIE), referred to as CBIE, and its outward normal BIE, referred to as HBIE. The formulations can be described as follows:

c(s)p(s)+∫ΓF(s,x)p(x)dΓ(x)=∫ΓG(s,x)∂p(x)∂n(x)dΓ(x)+pin(s), (8)

c(s)∂p(s)∂n(s)+∫ΓH(s,x)p(x)dΓ(x)=∫ΓK(s,x)∂p(x)∂n(x)dΓ(x)+∂pin(s)∂n(s), (9)

where s and x are the source point and field point, respectively; ∂p(x)/∂n(x)=q(x) denotes the flux of sound pressure; pin(s) indicates the sound pressure of the incident wave at point s; n(x) denotes the outward normal vector at point x, and jump term c(s) is equal to 1/2 if the boundary Γ is smooth around the source point s.

The Green’s function for 3D acoustics is given as follows [41]:

G(s,x)=eikr4πr, (10)

F(s,x)=-eikr4πr2(ikr-1)∂r∂n(x), (11)

K(s,x)=-eikr4πr2(ikr-1)∂r∂n(s), (12)

H(s,x)=eikr4πr3[(3-3ikr-k2r2)∂r∂n(s)∂r∂n(x)+(1-ikr)nj(s)nj(x)], (13)

where r=|s-x| denotes the Euclidean distance between points s and x.

In this study, we consider the impedance boundary condition

q(x)=ikβ(x)p(x),(14)

where β(x) denotes the normalized surface admittance at field point x.

In IGA BEM, the NURBS interpolations are applied to both the geometry and physical fields. In this study, we adopt different NURBS interpolation formulations to suit physical analysis, which means that the physical field is separated from the geometry. The knot vectors of the physical space can be represented as Ξd1 and Ξd2. Thus, the physical field variables can be expressed as follows:

p(ξ,η)=∑l=1NdRld(ξ,η)pl,q(ξ,η)=∑l=1NdRld(ξ,η)ql (15)

=ikβx∑l=1NdRld(ξ,η)pl, (16)

where pl and ql denote the sound pressure and flux coefficients associated with the l-th physical control point, respectively; Nd is the total number of physical control points; and βx represents the acoustic admittance at point x(ξ,η).

The location of collocation points in parametric space can be obtained by the Greville abscissa as follows:

ξc,i=ξi+1+ξi+2+…+ξi+pd1pd1,i=1,2,…,nd, (17)

ηc,j=ηj+1+ηj+2+…+ηj+pd2pd2,j=1,2,…,md, (18)

where pd1 and pd2 denote the orders of basis function of each dimension, respectively; nd and md are the number of collocation parameters in each dimension, and md×nd=Nd.

Similarly, after discretizing the boundary into Ne1×Ne2=Ne integral elements, Eqs. (8) and (9) can be rewritten as

c(s(ξc,i,ηc,j))∑l=1NdRld(ξc,i,ηc,j)pl+∑a=1Ne1∑b=1Ne2∑l=1Nd[∫ξaξa+1∫ηbηb+1F(s(ξc,i,ηc,j),x(ξ,η))Rld(ξ,η)J(ξ,η)dξdη]pl=∑a=1Ne1∑b=1Ne2∑l=1Ndikβx[∫ξaξa+1∫ηbηb+1G(s(ξc,i,ηc,j),x(ξ,η))Rld(ξ,η)J(ξ,η)dξdη]pl+pin(s)(19)

c(s(ξc,i,ηc,j))∑l=1NdikβsRld(ξc,i,ηc,j)pl+∑a=1Ne1∑b=1Ne2∑l=1Nd[∫ξaξa+1∫ηbηb+1H(s(ξc,i,ηc,j),x(ξ,η))Rld(ξ,η)J(ξ,η)dξdη]pl=∑a=1Ne1∑b=1Ne2∑l=1Ndikβx[∫ξaξa+1∫ηbηb+1K(s(ξc,i,ηc,j),x(ξ,η))Rld(ξ,η)J(ξ,η)dξdη]pl+∂pin(s)∂n(s),(20)

where [ξa,ξa+1]×[ηb,ηb+1] denotes the parameter space of an integral element; J(ξ,η) represents the Jacobian; (ξc,i,ηc,j) and (ξ,η) denote the two parameters of the source point s and field point x, respectively; βs denotes the acoustic admittance at the source point s(ξc,i,ηc,j), and we assume that the admittance keeps the same value over each element. The local support of NURBS basis functions, where a maximum of (pd1+1)×(pd2+1) basis functions are nonzero for values, can reduce the computation effort for the integral coefficient.

When the parameter (ξc,i,ηc,j) of the source point lies within the element parametric space [ξa,ξa+1]×[ηb,ηb+1], weak and hyper singularities appear in the kernel functions of Eqs. (19) and (20), respectively. Here, the detailed treatment of singular integrals is referenced in the use of Guiggiani method [60] by Chen et al. [41].

By adopting the Burton–Miller formulation [58] with coefficient α=-ik [59], Eqs. (19) and (20) can be combined and formulated as the matrix equation

(H-G)p=pin+αqin,(21)

where p is the complex vector of the physical values at the physical control points; and pin and qin denote the incident wave vectors of sound pressure and flux, respectively.

After Eq. (21) is solved, the sound pressure of points lying in the domain can be obtained by

pf=(Gf-Hf)p+pinf,(22)

where the computations of matrices Gf, Hf, and vector pinf are similar to those in Eq. (21).

3  Sensitivity Analysis through IGA BEM

Sensitivity analysis can obtain the derivatives of objection function with respect to different kinds of design variables. Thus, this type of analysis plays a critical role in the optimization process. The shape sensitivity and topology sensitivity analyses are presented in this section.

3.1 Shape Sensitivity Analysis

Control points usually control the configuration of structure surface in IGA, which means that they can be naturally set as the design parameters in shape optimization. In this study, we set certain control points of the NURBS surface as design parameters to change the shape, and the AVM proposed by Zheng et al. [54] is adopted for shape sensitivity analysis.

In DDM, Eqs. (8) and (9) can be differentiated with respect to the design parameter γe (usually represents the x, y, or z coordinate of a certain control point) as

c(s)∂p(s)∂γe+∫Γ∂F(s,x)∂γep(x)dΓ(x)+∫ΓF(s,x)∂p(x)∂γedΓ(x)+∫ΓF(s,x)p(x)d∂Γ(x)∂γe=∫Γ∂G(s,x)∂γeq(x)dΓ(x)+∫ΓG(s,x)∂q(x)∂γedΓ(x)+∫ΓG(s,x)q(x)d∂Γ(x)∂γe+∂pin(s)∂γe, (23)

c(s)∂q(s)∂γe+∫Γ∂H(s,x)∂γep(x)dΓ(x)+∫ΓH(s,x)∂p(x)∂γedΓ(x)+∫ΓH(s,x)p(x)d∂Γ(x)∂γe=∫Γ∂K(s,x)∂γeq(x)dΓ(x)+∫ΓK(s,x)∂q(x)∂γedΓ(x)+∫ΓK(s,x)q(x)d∂Γ(x)∂γe+∂∂γe(∂pin(s)∂n(s)), (24)

where the sensitivities of the kernel function are presented as

∂G(s,x)∂γe=eikr4πr2(ikr-1)∂r∂γe, (25)

∂F(s,x)∂γe=eikr4πr3[(2-2ikr-k2r2)∂r∂n(x)∂r∂γe+(ikr-1)r∂∂γe(∂r∂n(x))], (26)

∂K(s,x)∂γe=eikr4πr3[(2-2ikr-k2r2)∂r∂n(s)∂r∂γe+(ikr-1)r∂∂γe(∂r∂n(s))], (27)

∂H(s,x)∂γe=eikr4πr4(-9+9ikr+4k2r2-ik3r3)∂r∂n(s)∂r∂n(x)∂r∂γe-eikr4πr4(3-3ikr-k2r2)nj(s)nj(x)∂r∂γe+eikr4πr3(3-3ikr-k2r2)[∂∂γe(∂r∂n(s))∂r∂n(x)+∂r∂n(s)∂∂γe(∂r∂n(x))]+eikr4πr3(1-ikr)(∂nj(s)∂γenj(x)+nj(s)∂nj(x)∂γe).(28)

The sensitivity interpolation by NURBS basis functions are presented as

∂x(ξ,η)∂γe=∑l=1NgRl(ξ,η)∂Pl∂γe, (29)

∂p(ξ,η)∂γe=∑l=1NdRld(ξ,η)∂pl∂γe, (30)

∂q(ξ,η)∂γe=ikβx∑l=1NdRld(ξ,η)∂pl∂γe, (31)

where ∂pl∂γe denotes the sensitivity coefficient of the l-th control point. The other detailed sensitivity terms are discussed in the work of Chen et al. [41].

Similarly, using the same discretization as in Eq. (19), we can express Eqs. (23) and (24) as

c(s(ξc,i,ηc,j))∑l=1NdRld(ξc,i,ηc,j)∂pl∂γe=∑a=1Ne1∑b=1Ne2∑l=1Nd{∫ξaξa+1∫ηbηb+1[ikβx∂G(s(ξc,i,ηc,j),x(ξ,η))∂γe-∂F(s(ξc,i,ηc,j),x(ξ,η))∂γe]Rld(ξ,η)J(ξ,η)dξdη}pl+∑a=1Ne1∑b=1Ne2∑l=1Nd{∫ξaξa+1∫ηbηb+1[ikβxG(s(ξc,i,ηc,j),x(ξ,η))-F(s(ξc,i,ηc,j),x(ξ,η))]Rld(ξ,η)J(ξ,η)dξdη}∂pl∂γe+∑a=1Ne1∑b=1Ne2∑l=1Nd{∫ξaξa+1∫ηbηb+1[ikβxG(s(ξc,i,ηc,j),x(ξ,η))-F(s(ξc,i,ηc,j),x(ξ,η))]Rld(ξ,η)∂J(ξ,η)∂γedξdη}pl+∂pin(s)∂γe,(32)

c(s(ξc,i,ηc,j))∑l=1NdikβsRld(ξc,i,ηc,j)∂pl∂γe=∑a=1Ne1∑b=1Ne2∑l=1Nd{∫ξaξa+1∫ηbηb+1[ikβx∂K(s(ξc,i,ηc,j),x(ξ,η))∂γe-∂H(s(ξc,i,ηc,j),x(ξ,η))∂γe]Rld(ξ,η)J(ξ,η)dξdη}pl+∑a=1Ne1∑b=1Ne2∑l=1Nd{∫ξaξa+1∫ηbηb+1[ikβxK(s(ξc,i,ηc,j),x(ξ,η))-H(s(ξc,i,ηc,j),x(ξ,η))]Rld(ξ,η)J(ξ,η)dξdη}∂pl∂γe+∑a=1Ne1∑b=1Ne2∑l=1Nd{∫ξaξa+1∫ηbηb+1[ikβxK(s(ξc,i,ηc,j),x(ξ,η))-H(s(ξc,i,ηc,j),x(ξ,η))]Rld(ξ,η)∂J(ξ,η)∂γedξdη}pl+∂∂γe(∂pin(s)∂n(s)).(33)

Subsequently, still based on the Burton–Miller formulation, Eqs. (32) and (33) can be combined and expressed by matrix form as

(H-G)∂p∂γe=(∂G∂γe-∂H∂γe)p+∂pin∂γe+α∂qin∂γe,(34)

where matrices H and G have been computed by Eq. (21). ∂pin∂γe and ∂qin∂γe denote the vector of ∂pin∂γe and ∂∂γe(∂pin∂n), respectively. ∂G∂γe and ∂H∂γe are the derivatives of G and H. The detailed formulation of singular integrals in Eqs. (32) and (33) is presented in the work of Chen et al. [41].

For sensitivity analysis via DDM, after Eq. (34) is solved, the source point s is placed at the acoustic domain Ω, and their sensitivities can be calculated by

∂pf∂γe=(Gf-Hf)∂p∂γe+(∂Gf∂γe-∂Hf∂γe)p,(35)

where the computation of ∂Gf∂γe and ∂Hf∂γe is similar to those in ∂G∂γe and ∂H∂γe by Eq. (34).

However, some engineering problems may require more design parameters to guarantee the complexity of the structure, and the objective function in the optimization process may be merely related to physical values at field points. Thus, the AVM is adopted for ordinary optimization design. In accordance with Eqs. (34) and (35), the formulation can be rewritten as

∂pf∂γe=(Gf-Hf)(H-G)-1[(∂G∂γe-∂H∂γe)p+∂pin∂γe+α∂qin∂γe]+(∂Gf∂γe-∂Hf∂γe)p,(36)

where the calculation of (H-G)-1 usually requires extensive effort. To avoid this process, we introduce the following adjoint equation:

(H-G)TU=(Gf-Hf)T,(37)

and thus, the value of adjoint matrix U can be obtained by

UT=(Gf-Hf)(H-G)-1.(38)

Finally, by substituting Eq. (38) into Eq. (36), we can calculate the new formulation for the shape sensitivities at field points as follows:

∂pf∂γe=UT[(∂G∂γe-∂H∂γe)p+∂pin∂γe+α∂qin∂γe]+(∂Gf∂γe-∂Hf∂γe)Tp,(39)

Evidently, the adjoint matrix U is independent of all the design parameters, which means that Eq. (37) needs be solved once even if a large number of design parameters exist. Thus, the AVM can significantly improve the efficiency of shape sensitivity analysis.

3.2 Topology Sensitivity Analysis

In this study, we investigated the topology optimization of the distribution of SAM on structural surface, which is a problem with discrete values of 0 or 1, so that the mathematical computation in sensitivity analysis can be conducted. The SIMP method [46] is adopted to transform the discrete values into continuous values. The artificial density ρe of the e-th element is set as the design variable to conduct the topology optimization process. The formulation of admittance can be expressed as

βe=β0ρeμ,(40)

where μ denotes the penalization factor, and we define μ=3 in the analysis. β0 represents the normalized surface acoustic admittance. Here, we set its value as 0.1 for simplicity.

As the artificial density of the material of each element is set as a design variable, the AVM is also implemented for higher efficiency in topology sensitivity analysis. Usually, the objective function is related to the sound pressure of the field points as Π(pf), and in accordance with Eqs. (21) and (22), it can be rewritten as

Π̃=Π(pf)+ℜ{λ1T[(H-G)p-pin-αqin]+λ2T[pf-(Gf-Hf)p-pinf]},(41)

where ℜ represents the real part value of the complex vector. λ1T and λ2T are the adjoint vectors of the structure and field points, respectively, and they can be set to arbitrary values. Similarly, differentiating Eq. (41) with respect to the design variable ρe gives

∂Π̃∂ρe=∂Π∂ρe+ℜ{λ1T[-∂G∂ρep+(H-G)∂p∂ρe]+λ2T[∂pf∂ρe-∂Gf∂ρep-(Gf-Hf)∂p∂ρe]}.(42)

The derivatives of objective function Π(pf) can also be expressed as

∂Π∂ρe=ℜ{u1T∂p∂ρe+u2T∂pf∂ρe+u3},(43)

where the term u3 does not contain any term about ∂p∂ρe and ∂pf∂ρe.

Combining Eq. (42) with Eq. (43) yields

∂Π̃∂ρe=ℜ{[λ1T(H-G)+λ2T(Hf-Gf)+u1T]∂p∂ρe}+ℜ[(λ2T+u2T)∂pf∂ρe]+ℜ[-(λ1T∂G∂ρe+λ2T∂Gf∂ρe)p+u3].(44)