Adjoint design sensitivity analysis and optimization of nonlinear structures using geometrical mapping approach

Computers and Structures 183 (2017) 1–13

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

Adjoint design sensitivity analysis and optimization of nonlinear structures using geometrical mapping approach Yintao Wei a,⇑, Chonglei Zhao a, Zhenhan Yao b, Patrice Hauret c, Xuebing Li a, Michael Kaliske d a

State Key Lab Automotive Safety & Energy, Tsinghua University, Beijing 100084, China Department of Engineering Mechanics, Tsinghua University, Beijing 100084, China c Mécanique Numérique - CTEU/PM/SIM/ER, Centre de Technologies Michelin, 63040 Clermont-Ferrand, Cedex 9, France d Institute for Structural Analysis, Technische Universität Dresden, Dresden, Germany b

a r t i c l e

i n f o

Article history: Received 20 February 2016 Accepted 12 January 2017

Keywords: Shape optimization Adjoint variable method Nonlinear structures with contact forces Sensitivity analysis Geometrical mapping approach

a b s t r a c t A new adjoint shape design sensitivity formulation for nonlinear structures subject to contact forces is developed. The method is based on a geometrical mapping approach where shape variation is regarded as a mapping characterized by the shape variation velocity field. An adjoint variable method is developed for performing sensitivity analysis of the average shear strain in the tire belt area, which profile is properly parameterized in function of arcs, allowing explicit design velocities calculated. It can be seen that the present method is much faster than the direct differentiation method and more accurate than the classical finite difference scheme. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction Sensitivity analysis in shape optimization problems is complicated by the fact that the profile of the structure may change sharply in the optimization process [1]. This issue becomes even more critical in highly nonlinear problems such as, for example, optimal design of automotive tires whose mechanical response, durability and dynamical forces/moments generated in the contact patch strongly depends on the outer profile of the structure. Shape design sensitivity analysis (DSA) is still a challenge for tire industry researchers due to the high level of nonlinearity in material, geometry and contact stresses [2]. In this case, finite element analysis obtains an optimal design at affordable computational cost. Global optimization algorithms that combine different types of response surface approximations are more efficient than local search methods that may get trapped in local minima. The Kriging method is more robust than standard response surface fitting because it does not require assuming the form of the response surface a priori. The gradient-based Kriging (GBK) method [3] can further reduce the number of samples. However, the issue of disposing of a stable and accurate sensitivity formulation becomes even more critical in this case [3,4]. Another way to reduce the computational cost of shape optimization is to use sample points

⇑ Corresponding author. E-mail address: [email protected] (Y. Wei). http://dx.doi.org/10.1016/j.compstruc.2017.01.004 0045-7949/Ó 2017 Elsevier Ltd. All rights reserved.

determined by means of the Latin hypercube design minimizing pairwise correlations and maximizing inter-site distances of sample points [4,5]. The material derivative approach (MDA) (see the textbook by Choi and Kim [6,7]) has been adopted in different engineering areas by [8–20]. The adjoint variable method (AVM) [8,21–34] allows improving computational efficiency of sensitivity evaluations. However, the theoretical derivation of AVM for shape DSA of highly nonlinear structures as tires is a challenging task in view of the complexity of kinematic relationships. Here, a new adjoint shape DSA formulation is proposed for nonlinear tire structures with contact forces. The new formulation entails two basic aspects: (i) derivation of design velocity using a geometrical mapping approach (GMA) and (ii) formulation of shape DSA for nonlinear tire structures starting from the parameterized tire profile. The final goal is to perform accurate and stable DSA at low computational cost. The idea behind GMA is that shape variations are regarded as a mapping characterized by the shape variation velocity field. Using the virtual work principle, sensitivity equations of state variables (i.e., displacements) with respect to the shape variation are formulated. The AVM for sensitivity analysis of average shear strain in the tire belt area is hence formulated. The tire profile is parameterized using parameterized arcs, which makes boundary design velocities explicitly calculable according to its design variables. A passenger radial tire is chosen as case study: remarkably, AVM is

2

Y. Wei et al. / Computers and Structures 183 (2017) 1–13

about 80% faster than direct differentiation method (DDM) in this case. This provides a useful framework for shape sensitivity analysis and optimization of highly nonlinear structures. The paper is structured as follows. Section 2 describes design variables and objective function of the tire shape optimization problem. Section 3 presents the geometrical mapping approach that will be used for deriving sensitivity formulation. Section 4 illustrates how tire shape is parameterized and presents the definition of shape design velocity. Section 5 explains how to formulate sensitivity of displacements and displacement gradients with respect to shape variables. Sensitivity derivation for the objective function considered in this study is presented in Sections 6 and 7, respectively for the DDM and AVM. Applications to a passenger radial tire contour optimization are discussed in Section 8. Finally, Section 9 summarizes the main findings of this study.

2. Shape variables and objective function This article addresses nonlinear structures with contact forces, with tires being the application example. The tire shape optimization problem in this study included the five design variables shown in Fig. 1: radius of the first crown arc R1 , radius of the second crown arc R2 , radius of the shoulder arc R3 , radius of the sidewall-shoulder arc R4 , meridian coordinate of the intersection point of the first and second arcs W 1 , tread contact width W 2 . Since R4 is not an independent variable as its value can be known from W 2 , the radius related to the shape of the inner boundary and parameters remain constant in the optimization.

The objective function to be minimized is the average effective shear strain which develops in the belt end area near the tireground contact region. The corresponding area in the meridian plane is shown in Fig. 2, and the 3D finite element model of tire is obtained from this 2D model by rotating around the center of tire. The circle of 360 degrees is divided into 50 sections on average and the footprint length of tire tread in the circumferential direction comes to 6 sections. Therefore the circumferential angle is taken as 2p
6=50 which is shown in Fig. 4. The objective function (also called performance function) can be expressed as t

1

sw

¼

1 sV b

Z

t
sVb

s e ðxÞds V

ð1Þ

where the left subscript s denotes the initial configuration corresponding to the shape parameter sa (a = 1, 2. . .5) mentioned in Fig. 1, the left superscript t denotes the load parameter, s V b is the total volume of the design area, and ts
e ðxÞ [35] is the effective shear strain:

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   ffi 2t t 2 t 1t 1 e e ¼ e  e dij ts eij  ts ell dij se ¼ 3 s ij s ij 3 s ij 3 s kk 3 t


ð2Þ

Fig. 2 shows the FE mesh adopted for the meridian plane of the tire. The model is comprised of different rubber materials and rebars elements as reinforcements. The FE model, homogeneously divided into 50 sections in the circumferential direction, includes in most part trilinear hexahedral elements (HEX8). The level of mesh refinement was such to have mesh independent solutions at reasonably low computational cost. The structure is subject to two independent load cases: (i) inflation pressure and (ii) contact with rigid road surface under the load carried. 3. Geometrical mapping approach

Fig. 1. Design variables in the shape optimization of a tire (units/mm).

In order to evaluate the sensitivity of cost function with respect to shape variations, strain sensitivities must be determined. For that purpose, sensitivity of displacements with respect to changes in shape must be obtained first. In this study, the GMA was utilized because the MDA may be misleading as the real material derivative is not suitable for formulating the variation in shape design. The partial derivative with respect to the shape parameter does not have any physical meaning, and the gradient is the derivative with respect to the Euler coordinates rather than to the Lagrange coordinates that appear in the displacement gradient of the solids. The GMA is presented with clear physical meaning, wherein the shape variation is regarded as a mapping characterized by the shape vari-

Fig. 2. FE model of the meridian plane of the tire and definition of the cost function in belt end area.

3

Y. Wei et al. / Computers and Structures 183 (2017) 1–13

ation velocity field, or a fictitious deformation of the continuum (i.e. not a real displacement) and there are no strains or stresses corresponding to such fictitious displacement fields. The shape variation is formulated as a geometrical mapping, characterized by the shape variation velocity or shape design velocity field 0 V ia . The mapping of coordinates can be written as 0

s xi

¼ 0 xi þ 0 V ia sa ¼ 0 xi þ s ui 0

i ¼ 1;    ; m; a ¼ 1;    ; k

ð3Þ

where 0 xi , 0s xi denote the Lagrange coordinates for the original shape and the shape design corresponding to parameters sa respectively, m is the number of independent degrees of freedom (DOFs), k denotes the number of design parameters, and 0s ui can be regarded as a virtual displacement, which is not a real displacement of the continuum. In addition, there are no strains or stresses corresponding to such a virtual displacement field. For a given loading condition, the relationship between the displacements of the corresponding mapping points can be expressed as t

s ui

¼ t0 ui þ t0 ui

a

ð4Þ

sa

where the right superscript a denotes @=@ sa , which means the change rate of the variable for the corresponding mapping point a

during the shape variation, and t0 ui at sa ¼ 0 is the sensitivity of displacements. The relationship between the displacement gradients of the corresponding mapping points can be expressed as s ui ;j t

a

¼ t0 ui ;j þ t0 ðui ;j Þ

ð5Þ

sa



where t0 ðui ;j Þ a denotes the sensitivity of displacement gradients, subscript j indicates the derivative of the displacement ui with respect to the jth coordinate, which can be formulated as:

  a a @ t t  ðu ; Þ , ¼ ð u ; Þ ¼ t0 ui ;j  t0 ui ;k 0 V ka ;j i i j j 0  @ sa s s¼0

ð6Þ

In which dds ½ts ðui ;j Þ consists of two parts: one is due to the differa

ence of ts ui and t0 ui , namely t0 ui ; the other one comes from the difference of ts ðui ;j Þ and t0 ðui ;j Þ in the case of identical ts ui and t0 ui . 4. Tire shape parameterization and shape design velocity calculation An important concept in the derivation of sensitivity is the shape design velocity. For the shape design, the boundary shape design velocity can be directly obtained from the definition used to map the coordinates (Eq. (3)): 0 V ia

¼

@ 0 ð x  0 xi Þ @ sa s i

ð7Þ

where the shape velocity field V ia represents the variation relative of the ith degree of freedom defining the shape with respect to the ath design parameter (variable). Provided ð0s xi  0 xi Þ can be formulated explicitly using the design variables sa . Otherwise, the finite difference approach can be followed, 0 V ia

¼ ð0s xi ðDsa Þ  0 xi Þ=Dsa

ð8Þ

where Dsa is a small increment of a design variable. As mentioned in ection 2, we use five design variables, i.e., R1 , R2 , R3 , R4 , W 1 , and W 2 , as shown in Fig. 1, to describe the outer profile for a series of 50–70 typical passenger car radial tires (PCR). Moreover, note that the specified data in Fig. 1 for each parameter are simply an example, and R4 can be determined by W 2 . The controlling equations of all the arcs can explicitly be described as follows:

ðx1  xo1 Þ2 þ ðy1  yo1 Þ2 ¼ R21 ; 2

2

ðW 1  xo1 Þ þ ðy2  yo1 Þ ¼
 yo2 Þ ¼ ðW 2  xo2 Þ þ ðy 2

yo4 ¼ D  H;

2

xo1 ¼ 0;

R21 ;

R22 ;

yo1 ¼ D  R1

ðW 1  xo2 Þ2 þ ðy2  yo2 Þ2 ¼ R22
 yo4 Þ2 ¼ R24 ; ðW 2  xo4 Þ2 þ ðy

xo4 ¼ B  R4

ðx3  xo2 Þ2 þ ðy3  yo2 Þ2 ¼ R22 ;

ðx3  xo3 Þ2 þ ðy3  yo3 Þ2 ¼ R23 ;

ðx4  xo3 Þ2 þ ðy4  yo3 Þ2 ¼ R23 ; ðx4  xo4 Þ2 þ ðy4  yo4 Þ2 ¼ R24 W 1  xo2 W 1  xo1 x3  xo2 x3  xo3 x4  xo3 x4  xo4 ¼ ; ¼ ; ¼ y2  yo2 y2  yo1 y3  yo2 y3  yo3 y4  yo3 y4  yo4 ð9Þ where B is the tire width; D is the diameter; xoi , yoi are the coordinates of the ith arc center; xi , yi are outer surface point coordinates belonging to the ith arc. This means that, giving the five design variables, the tire width B and the diameter D, all 4 arcs as well as the intersection points can be calculated; therefore, the tire shape can be determined. The design velocity calculation is a key to the shape optimization problem, which is closely related to shape parameterization. A shape change can be described by a mapping function T expressed by

xs  Tðx; sÞ;

Xs  TðX; sÞ

ð10Þ

where X represents the whole shape domain. The design velocity, physically, can be understood as the partial derivative of the shape coordinates with the shape parameter. If one imagines this partial derivative as ‘‘the design velocity”, the shape parameter s can be taken as the ‘‘time”, therefore the design velocity can be defined as

Vðxs ; sÞ 

dxs dTðx; sÞ @Tðx; sÞ ¼ ¼ ds @s ds

ð11Þ

It can be shown that the design velocity is the key parameter for mesh updating and to evaluate sensitivity of objective function with respect to shape variables. Because all the arcs can be explicitly formulated using Eq. (9), the tire contour can be parameterized using arcs. That is:



 px ðuÞ py ðuÞ ¼ FBðRÞ

ð12Þ

where u are the parametric coordinates (i.e. the length of each arc); F and B are coefficient matrices, and R is the vector of design variables. The real meaning of Eq. (12) is that one can explicitly calculate all outer surface point coordinates, from which the design velocity can be calculated. Then, the design velocity of each mesh point can be determined by

V ¼F

@B @R

ð13Þ

After the 2D design velocity on the meridian plane is obtained, the 3D design velocities can be calculated by (see Fig. 3)

V x ¼ V x;

V y ¼ V y  cosðhÞ;

V z ¼ V y  sinðhÞ

ð14Þ

The boundary design velocity can be easily calculated with Eq. (13). Various approaches have been presented in the literature [8,9], such as the boundary displacement method and the fictitious load method, for determining the domain shape design velocity field based on the boundary design velocity. The resulting domain velocity field is almost distributed across the domain, and the computation is time-consuming in general. Conversely, in the present approach, the domain velocity field is localized in one layer of finite elements, which is directly related to the design modified boundary. For HEX8 elements, the domain velocity field of the elements

4

Y. Wei et al. / Computers and Structures 183 (2017) 1–13

where sa and sb denote the ath and bth shape design variables. Note that the strain expressed in Eq. (17) is defined according to the nonlinear theory. For the virtual displacement ds ui , which is arbitrary and geo
including all metrically admissible for the design shape vector s shape variables, the corresponding virtual strain can be expressed as follows:

h i. ds eij ¼ ds ðui ;j Þ þ ds ðuj ;i Þ þ ds ðuk ;i Þts ðuk ;j Þ þ ts ðuk ;i Þds ðuk ;j Þ 2 1n ðd0 ui ;j þ d0 uj ;i þ d0 uk ;i t0 uk ;j þ t0 uk ;i d0 uk ;j Þ ¼ 2     þ d0 ðui ;j Þ a þ d0 ðuj ;i Þ a þ d0 ðuk ;i Þ a t0 uk ;j þ t0 ðuk ;i Þ a d0 uk ;j o   ð18Þ þd0 uk ;i t0 ðuk ;j Þ a þ t0 uk ;i d0 ðuk ;j Þ a sa Fig. 3. Sketch of the tire for 3D design velocities.

is directly determined by the boundary velocity field, and no additional computing cost is required. Furthermore, undesirable mesh distortions for composite structures can be naturally avoided.

where the small higher order terms are neglected. The full derivation of Eq. (18) is given in Appendix C. For the second Piola-Kirchhoff stress tensor, it can be written as ¼ t0 Sij þ t0 Sij

t

s Sij t 0 Sij

a

a

sa     sa t0 DTijkl ðt0 ðuk ;l Þ a þ t0 ðul ;k Þ a þ t0 ðum ;k Þ a t0 um ;l þ t0 um ;k t0 ðum ;l Þ a Þsa =2

5. Sensitivity of displacements and displacement gradients

ð19Þ

The system of equations for the sensitivity of displacements can be derived from the virtual work principle. For the shape-varied structure in the shape optimization problem, the virtual work equation can be formulated as

Z Z t t t t LM S d e d V ¼ t d u d S þ s qs f k d0 uk ds V þ s C c s ij s ij s s k 0 k s t sV sS sV Z t LM C ¼ ðts kn d0 un þ ðts un þ 0s g n Þd0 kn Þds S s c Z

t

ð15Þ

t st s Sc

In this equation ts C LM is the contact boundary term in Lagrange c multiplier approach, where Sst c denotes the actually contacted part of the possibly in contact boundary, 0 kn , ts kn denote the Lagrange multipliers while g n denotes the normal gap, ui ði ¼ 1; 2; 3Þ denotes the displacement, the left subscript s denotes the initial configuration corresponding to the shape parameter sa , and the left superscript t denotes the loading condition. The entire boundary is divided into a traction-free boundary ts Sf , a traction-given boundary s St , and the contacted boundary ts Sst . The traction-free boundary, including the part of the tire boundary which may come in contact with the load, is excluded from the virtual work equation. It should be noted that the division of the boundary is load and shape dependent. Here, for simplicity, frictionless contact is considered; thus, only a normal Lagrange multiplier is used. For the shape optimization problem, the configuration corresponding to s ¼ 0 is taken as reference. Hence, Eq. (15) can be rewritten as

Z

Z

t s s Sij ds eij j0 J jd0 V

t 0V

Z

¼

s t k ds uk ðds S=d0 SÞd0 S þ t

0S

t

Z

þ

t Sst 0 c

t

0V

t t s s qs f k ds uk j0 J jd0 V

ðts kn d0 un þ ðts un þ 0s g n Þd0 kn Þðds S=d0 SÞd0 S ð16Þ

where the left term represent the stress increment, which is obtained from the product between the fourth order tangent constitutive tensor DTijkl ; i; j; k; l ¼ 1; 2; 3 and the displacement gradient sensitivity tensor. As far as it concerns the other terms in Eq. (16), it can be mentioned that sq

¼ 0q

ð20Þ

which implies keeping the same mass density during shape changes, for a generic element e:

ds V e ¼ js0 J jd0 V e ds V e ¼ js0 J jd0 V e ¼ js0 J jj0 Jjdn10 dn20 dn30 ds V e ¼ js0 J jd0 V e ¼ js0 J jj0 Jjdn10 dn20 dn30 !   a js Jj  j0 Jj 0 jJj ¼ 1þ js0 J j ¼ 1 þ sa ¼ ð1 þ 0 V ia ;i sa Þ j0 Jj j0 Jj

ð21Þ

Introducing the above equations into the equation of virtual work corresponding to the shape variation - Eq. (16) - it follows: XZ

h

0Ve

e

t 0 Sij

þ t0 DTijkl ðt0 ðuk ;l Þ

a

a

þ t0 ðul ;k Þ

a

þ t0 ðum ;k Þ

t 0 um ;l



þ t0 um ;k t0 ðum ;l Þ a Þsa =2

i

h  ðd0 ui ;j þ d0 uj ;i þ d0 uk ;i t0 uk ;j þ t0 uk ;i d0 uk ;j Þ=2       þ d0 ðui ;j Þ a þ d0 ðuj ;i Þ a þd0 ðuk ;i Þ a t0 uk ;j þ t0 ðuk ;i Þ a d0 uk ;j þ d0 uk ;i t0 ðuk ;j Þ a i  þt0 uk ;i d0 ðuk ;j Þ a sa =2 ð1 þ 0 V ia ;i sa Þd0 V e Z X XZ a t t ¼ 0 q0 f k ð1 þ 0 V ia ;i sa Þd0 uk d0 V e 0 t k ð1 þ 0 Se sa Þd0 uk d0 Se þ e

þ

t 0 Se

XZ e

t Sst 0 c

0Ve

e

fðt0 kn a

þ ðt0 un þ t0 un

þ

t 0 kn

a

b

sa Þd0 un ð1 þ 0 Se sb Þ 

sa þ 0 0 g n Þd0 kn ð1 þ 0 Seb sb Þgd0 Se

ð22Þ

Recalling Eq. (5), the strain field can be expressed as follows:

. 1 nt t t t t t t 2¼ s ui ;j þ s uj ;i þ s uk ;i s uk ;j 0 ui ;j þ 0 uj ;i þ 0 uk ;i 0 uk ;j 2      þ t0 ðui ;j Þ a þ t0 ðuj ;i Þ a þ t0 ðuk ;i Þ a t0 uk ;j þ t0 uk ;i t0 ðuk ;j Þ a sa  o   þ t0 ðuk ;i Þ a t0 ðuk ;j Þ b sa sb ð17Þ

s eij ¼ t



t

Note that Se denotes the element surface. Developing this equation with respect to the shape variation factor sa , after having neglected the small higher order terms related to the stress field t0 Sij and the loads t0 t k , 0 qt0 f k , which are in equilibrium in the original shape, and introducing Eq. (6) for the sensitivity of the displacement gradient, Eq. (22) can be rewritten as

5

Y. Wei et al. / Computers and Structures 183 (2017) 1–13

XZ



0Ve

e

þ

The load vector in Eq. (24) can be formulated as follows:



½ðd0 ui ;j þ d0 uj ;i Þ=2t0 DTijkl ½ðt0 uka ;l þ t0 ul a ;k Þ=2d0 V e

XZ

½ðd0 ui ;j þ

0Ve

e

t a 0F

 d0 uj ;i Þ=2t0 DTijkl ½ðt0 um ;k t0 uma ;l



þ t0 um ;l t0 uma ;k Þ=2d0 V e þ

XZ 0Ve

e

0Ve



þ ¼

t t a 0 um ;l 0 um ;k Þ=2d0 V e

0Ve

e

XZ

 ðd0 uk Þðt0 tk Þ0 Sea d0 Se

t 0 Se

e

þ

XZ


ð0 qt0 f k Þ0 V ia ;i d0 V e 

0Ve

þ

XZ

0Ve

þ

XZ

XZ

d0 uj ;i Þ=2Þt0 DTijkl ððt0 uk ;m 0 V ma ;l ½ððd0 ui ;j þ

0Ve

ð0 K L þ t0 K N þ t0 K G Þt0 u

þ

þ

t a 0F F

e

½ððd0 un ;k 0 V ka ;i t0 un ;j þ ½ððd0 ui ;j

0Ve

þ

t 0 ul ;m 0 V ma ;k Þ=2Þd0 V e

XZ

XZ

0Ve

e a

0Ve

t t 0K L þ 0K N þ 0K G 0u a t t a ¼g 0 GC 0 u

¼ t0 F

a

þ

½ððd0 un ;i t0 un ;j

0KL

¼

e

t 0KN

¼

0V

XZ e

½ððd0 un ;i t0 un ;j

t Sst 0 c a

a

a

þ t0 F C  t0 C C t0 k

t 0KG

¼

XZ e

Tt 0 BL 0 D 0 BL d0 V

ð23Þ

ð24Þ

ð25Þ

T

0V

0V

T

ð0 BTL t0 D t0 BN þ t0 BN t0 D 0 BL þ t0 BN t0 D t0 BN Þd0 V

Tt 0 BG 0 S 0 BG d0 V

a

þ t0 K G Þt0 u

a

ð30Þ

0Ve

XZ

0Ve

XZ

0Ve

T t 0 BGVa 0 S d0 V e

þ

T t 0 BNVa 0 S d0 V e

þ

XZ e

0Ve

e

0Ve

T t 0 BLVa 0 Sd0 V e

XZ

Tt 0 BL 0 DENVa d0 V e

þ

Tt 0 BL 0 DELVa d0 V e

XZ e

0Ve

Tt 0 BN 0 DELVa d0 V e

Tt 0 BN 0 DENVa d0 V e

ð31Þ

where all the terms can be derived from Eq. (23). In Eqs. (25)–(27), stiffness matrix t0 K N , t0 K N , t0 K N including strain matrix BL , BN , BG preserve their classical meaning. However, the strain matrix parameters BLVa , BNVa , ELVa , ENVa related to design velocity in Eq. (31) are different from those defined in conventional finite element analysis as explained in Appendix D. The complete derivation is a routine process and omitted here for the sake of brevity. On the other hand, the change rates of the real loads can be calculated as

where GC is the kinematic matrix and C C denotes the contact matrix with regard to Lagrange multiplier, 0 K L , t0 K N , t0 K G are identical to those in the conventional nonlinear finite element analysis for the original design shape. That is:

XZ

XZ

e



a t 0F R

 t0 Q

ð26Þ

ð27Þ

a

¼

XZ e

Finally, the equations for evaluating sensitivity of displacements can be formulated in matrix form as follows: a

0Ve

e

þ d0 kn0 0 g n þ d0 kn t0 un þ ðd0 un t0 kn Þ0 Se Þd0 Se

t



a

ð0 BTL þ t0 BN Þt0 S 0 V ia ;i d0 V e

XZ

e

þ d0 un ;j t0 un ;i Þ=2Þt0 DTijkl ððt0 um ;k t0 um ;n 0 V na ;l XZ a þ t0 um ;l t0 um ;n 0 V na ;k Þ=2Þd0 V e þ ðd0 un t0 kn





T

0Ve

e

d0 un ;j t0 un ;i Þ=2Þt0 DTijkl ððt0 uk ;m 0 V ma ;l

a

XZ ¼ þ

XZ

e

e

a

¼ t0 F R  ð0 K La þ t0 K N a

d0 uj ;i Þ=2Þt0 DTijkl ððt0 um ;k t0 um ;n 0 V na ;l

þ t0 ul ;m 0 V ma ;k Þ=2Þd0 V e þ

a

¼ t0 F R þ t0 F Fa ¼ t0 F

½ððd0 ui ;k t0 V ka ;j

0Ve

e

0Ve

þ t0 um ;l t0 um ;n 0 V na ;k Þ=2Þd0 V e þ þ

ð29Þ

sa



d0 un ;k 0 V ka ;j t0 un ;i Þ=2Þt0 Sij d0 V e

e



In the sensitivity analysis, the system of equations to be solved

ðd0 uk Þ

e

þ



sa þ t0 K G þ t0 K G a sa Þðt0 u þ t0 u a sa Þ

where the fictitious load t0 F Fa can be formulated as

d0 un ;j t0 un ;k 0 V ka ;i Þ=2Þt0 Sij d0 V e

XZ

a

is

½ðd0 ui ;j þ d0 uj ;i þ d0 un ;i t0 un ;j

e

þ

¼ t0 F þ t0 F R

½ððd0 un ;i t0 un ;k 0 V ka ;j

þ d0 uj ;k t0 V ka ;i Þ=2Þt0 Sij d0 V e þ þ

a

d0 un ;j t0 un ;i Þ=2t0 Sij 0 V ka ;k d0 V e e

þ



0Ve

e

XZ e

þ

þ

ð28Þ

a

ð0 K L þ 0 K La sa þ t0 K N þ t0 K N

 ðd0 un ;i Þt0 Sij ðt0 una ;j Þd0 V e

XZ



þ t0 F Fa

where t0 Q , t0 F F denote the change rates of the real loads and the fictitious loads, respectively. For the shape-modified design, the system of equations in the nonlinear finite element analysis can be written as

þ d0 un ;j t0 un ;i Þ=2t0 DTijkl ½ðt0 uka ;l þ t0 ul a ;k Þ=2d0 V e XZ  þ ½ðd0 un ;i t0 un ;j þ d0 un ;j t0 un ;i Þ=2t0 DTijkl ½ðt0 um ;k t0 uma ;l e

a a

½ðd0 un ;i t0 un ;j



¼ t0 Q

0Ve

N T ð0 qt0 F Þ0 V ia ;i d0 V e þ

XZ e

0 Se



N T t0 t 0 Sea d0 Se ð32Þ

Sensitivity of displacements can be determined from the linear system (30). 6. Determination of cost function sensitivities with the direct differentiation method The sensitivities of cost function defined by Eq. (1) can be computed once sensitivities of displacements and the displacement gradients are known. Theaverage strain functional (1) defined for the tire shape optimization problem can be expressed as

Z t 0w

¼ 0V

f ðt0 ui ð0 xÞ; t0 ui ;j ð0 xÞÞd0 V

ð33Þ

If we assume that density remains unchanged during the shape variation, the cost function can be formulated as

6

Y. Wei et al. / Computers and Structures 183 (2017) 1–13

Z t

sw

¼ sV

Z

The adjoint equation system can hence be formulated as

f ðts ui ðs xÞ; ts ui ;j ðs xÞÞds V



a

¼ 0V

f ðt0 ui ð0 xÞ þ t0 ui ð0 xÞsa ; t0 ui ;j ð0 xÞ a

þ t0 ðui ;j Þ where

b

a

sa Þ½1 þ ð0 jJj =j0 JjÞsb d0 V t0 w þ t0 w sa

ð34Þ

#  8   8 @f X @f X m a m a m m t t ¼ N ðnÞ0 ui þ d0 V e 0 ðN ;j Þ0 ui @ui m¼1 @ui ;j m¼1 0Ve e ! ! 8 8 X X Z  @f  X m m n t n t  0 ðN ;k Þ0 ui 0 ðN ;j Þ0 V ka d0 V e @u ; i j 0Ve e m¼1 n¼1 8 XZ X þ f ðt0 ui ð0 xÞ; t0 ui ;j ð0 xÞÞ 0 ðN m ;i Þ0 V m ia d0 V e 0Ve

e

m¼1

ð35Þ Introducing

#  8   8 @f X @f X m a m a m m t t N ðnÞ0 ui þ d0 V e , 0 ðN ;j Þ0 ui @ui m¼1 @ui ;j m¼1

0Ve

e

¼



e

þ

@f @ui ;j

0Ve

XZ

0Ve

e

 X 8

! m t 0 ðN ;k Þ0 ui

m

!

8 X n¼1

ð37Þ

m¼1

@ t0 w t  a t
 a u þ 0w @ t0 u 0

¼

¼

a

ð38Þ 

Displacement sensitivities t0 u a are obtained by solving Eq. (30). Determination of cost function sensitivities with Eq. (38) and (30) is called direct sensitivity calculation (DDM). Since this method is computationally expensive as it requires solving the system (30), the AVM was developed in this study. 7. Determination of cost function sensitivities with the adjoint variable method



K

"

 ¼

a

þ t0 K N þ t0 K G

0KL

T t 0 GC

T t 0CC

0K L

þ t0 K N þ t0 K G

t 0CC

t 0 GC

0

( t

a

0u t a 0k

)

(

¼

a t 0F



þ t0 F Ca

g

)

ð39Þ

a

The sensitivity of cost function can be determined with the following expression t a 0w

¼

h

@ t0 w @ t0 u

@ t0 w @ t0 k

i

0K L

þ t0 K N þ t0 K G t 0 GC

t 0CC

1 ( t

0F

a

0

a

þ t0 F C

g

a

)

a


þ t0 w

ð40Þ where

@ t0 w @ t0 k

¼ 0 because the cost function does not directly depend on

the traction at the contact surface. Using the adjoint approach, let us put

h

@ t0 w @ t0 u

@ t0 w @ t0 k

i K þ t K þ t K 0 L 0 N 0 G t 0 GC

t 0CC

0

1

 ¼ KT

aT



@ t0 u

:

ð42Þ

;

0

#1 8  t T 9 < @0 w = :

0

@ t0 u

0

ð43Þ

;

  After the adjoint variables KT aT are solved from Eq. (43), and then introducing Eq. (41) into Eq. (40), the sensitivity of cost function is hence expressed as: t a 0w

¼ KT ft0 F

a

a


þ t0 F Ca g þ aT fg a g þ t0 w 



ð44Þ

Z 1 t
s e ðxÞds V sV b sVb rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   ffi 2t t 2 t 1t 1 t
e e ¼ e  e dij ts eij  ts ell dij se ¼ 3 s ij s ij 3 s ij 3 s kk 3

0

W1

 a

¼ KT

t

0F

a

  

  
a þ t0 F Ca þ aT g a þ t0 W 1

ð45Þ

ð41Þ

ð46Þ

If the design velocity is zero in a given region of design space, cost function will be insensitive to shape variations. In the adjoint system (43), the right side term is called adjoint load. In the present case, it can be computed as follows !T " #T @ t0 ðw1 Þ @ t0 ðw1 Þ @ t0 ðw1 Þ @ t0 ðw1 Þ @ t0 ðw1 Þ @ t0 ðw1 Þ ¼       @ t0 u @ t0 u1ð1;1Þ @ t0 u3I2ðI;1Þ @ t0 u3I1ðI;2Þ @ t0 u3IðI;3Þ @ t0 u3NðN;3Þ ( nb Z 1 Z 1 Z 1 X  @ t0 ðw1 Þ 1 1 t
e ðxÞ 1 ¼ P nb n @ t0 u 3 0 n¼1 0 V e n¼1 1 1 1 bðI; kÞ I2en

b ¼ 3ðI  1Þ þ k 


t 0 ij

e

2  X 8 8 X 1 t 6 p p  0 ehh dij 4 0 ðN ;j Þdik þ 0 ðN ;i Þdjk 3 p¼1 p¼1 ðn;pÞ¼I

The system (24) that allows sensitivity of displacements with respect to shape variables to be computed can be rewritten as follows



8 9 < @ t0 w T =

The matrix of coefficients of the adjoint equation system is the transpose of the original system of equations because the matrix of coefficients is not symmetric. Then,

t

8 X  t
a f ðt0 ui ð0 xÞ; t0 ui ;j ð0 xÞÞ 0 ðNm ;i Þ0 V m ia d0 V e , ¼ 0 w

The sensitivity of cost function can hence be expressed as t a 0w



Sensitivities with respect to shape variables can be determined as follows

n t n 0 ðN ;j Þ0 V ka d0 V e

m¼1

0

K

1 sw ¼

@w t  a u @ t0 u 0

XZ

T 

t

ð36Þ 

t 0CC

For example, for the average strain to be determined in the belt end area of the tire - Eq. (1) - the cost function is

"

XZ

þ t0 K N þ t0 K G t 0 GC

"

XZ

t a 0w

0KL

þ

8 X

0 ðN

p¼1 ðn;pÞ¼I

p

;i Þdkl

8 X q¼1

0 ðN

q

;j Þ0 t uql þ

ðn;pÞ¼I

8 X

0 ðN

p

p¼1

;i Þ0 t upl

8 X

0 ðN

q

;j Þdkl

q¼1 ðn;qÞ¼I

0 1 3 ) 8 8 8 X X 2 BX C 7 p p q t p d j0 Jjdn10 dn20 dn30  @ ðN ; Þd þ ðN ; Þ u ðN ; Þd A 5 0 0 0 ij l kl l 0 m l km 3 p¼1 q¼1 p¼1 ðn;pÞ¼I

ðn;qÞ¼I

ð47Þ

where the subscript b of the displacement u is the global number of the node displacement component for all nodes related to the cost function, which is denoted as bðI; kÞ for the I-th node k-th degree of freedom, and takes the value b ¼ 3ðI  1Þ þ k; n is the number of element related to the cost function, the summation on all related boundary elements, number 1 to nb , the summation performed only on the elements connected to the I-th node; q or p is the local number of the element node, the condition ðn; pÞ ¼ I under summation symbol means that, if the summation conducted only for the node in element n, coincides with the global I-th node; N p or Nq denote the shape functions of the p-th or q-th element node.

7

Y. Wei et al. / Computers and Structures 183 (2017) 1–13

(48)

In order to simplify Eq. (48), it can be written

(49)

The terms indicated in red are eliminated. In summary, the AVM requires to calculate the adjoint load and solve Eq. (43). Therefore, its computational cost is much lower than in the case of DDM if the number of optimization variables is larger than the number of objective functions.

8. Determination of sensitivities in the shape optimization of a passenger car radial tire The sensitivity analysis technique described in this study was integrated with a finite element program available to authors. Computational tasks are distributed as follows:

8

Y. Wei et al. / Computers and Structures 183 (2017) 1–13

Table 1 Partitions of design space defined for the 205/55R16 tire shape optimization problem.

Level Level Level Level Level Level Level Level

1 2 3 4 5 6 7 8

W1

W2

R1

R2

R3

40.6 42.02857 43.45714 44.88571 46.31429 47.74286 49.17143 50.6

90.3 91.01429 91.72857 92.44286 93.15714 93.87143 94.58571 95.3

600 685.7143 771.4286 857.1429 942.8571 1028.571 1114.286 1200

248 298.2857 348.5714 398.8571 449.1429 499.4286 549.7143 600

40 41.85714 43.71429 45.57143 47.42857 49.28571 51.14286 53

Fig. 4. 3D finite element model of 205/65R16 tire.

Table 2 Results of sensitivity analysis carried out for the 205/55R16 tire shape optimization problem. DOE sample

Design variables

Value of design variables

waA

waB

waC

w1

waAVM

waFDM

waDDM

#1

1 2 3 4 5

4.06E+01 9.46E+01 1.03E+03 3.49E+02 5.30E+01

1.46E05 9.89E07 6.21E07 6.98E09 1.64E04

4.61E07 2.63E07 1.30E07 7.80E08 1.21E04

9.23E05 7.55E07 3.09E06 3.92E06 4.54E04

4.99E02 4.99E02 4.99E02 4.99E02 4.99E02

7.72E05 2.77E08 2.34E06 3.99E06 4.10E804

7.44E05 1.40E06 2.33E06 3.22E06 1.18E04

7.72E05 2.77E08 2.34E06 3.99E06 4.10E04

#2

1 2 3 4 5

4.20E+01 9.32E+01 6.00E+02 4.49E+02 4.00E+01

3.10E06 86.96E05 1.31E06 1.37E09 2.95E04

4.92E07 2.99E05 6.50E07 1.88E07 7.20E05

5.82E05 1.01E04 1.71E05 8.12E06 1.67E04

4.84E02 4.84E02 4.84E02 4.84E02 4.84E02

5.46E05 1.41E04 1.52E05 8.31E06 3.89E804

5.89E05 2.69E04 1.72E05 8.33E06 2.88E04

5.46E05 1.40E04 1.52E05 8.31E06 3.89E804

#3

1 2 3 4 5

4.35E+01 9.03E+01 8.57E+02 4.99E+02 4.93E+01

6.03E06 6.86E07 9.31E07 7.11E08 1.93E04

81.02E806 4.84E06 3.19E07 1.04E08 1.06E04

83.87E05 1.11E05 4.92E06 1.76E06 3.77E04

5.04E02 5.04E02 5.04E02 5.04E02 5.04E02

3.17E05 1.53E05 3.67E06 1.68E06 2.90E04

3.27E05 2.03E05 4.04E06 1.45E06 2.00E04

3.17E05 1.53E05 3.67E06 1.68E06 2.90E04

#4

1 2 3 4 5

4.49E+01 9.17E+01 1.20E+03 3.99E+02 4.37E+01

8.71E06 1.25E05 3.92E07 5.07E08 2.60E04

1.47E06 2.14E05 1.72E07 9.13E08 7.96E05

8.98E05 3.54E05 2.91E06 3.62E06 2.56E04

5.25E02 5.25E02 5.25E02 5.25E02 5.25E02

7.96E05 4.43E05 2.35E06 3.66E06 7.64E05

8.71E05 7.87E05 2.80E06 3.43E06 1.98E04

7.96E05 4.43E05 2.35E06 3.66E06 7.64E05

#5

1 2 3 4 5

4.63E+01 9.10E+01 9.43E+02 2.48E+02 4.56E+01

1.72Ev05 1.33E05 6.96E07 2.36E07 2.33E04

2.33E06 1.23E05 1.58E07 5.45E07 6.59E05

1.63E04 1.23E05 4.98E06 9.33E06 1.76E04

5.07E02 5.07E02 5.07E02 5.07E02 5.07E02

1.49E04 3.79E05 4.13E06 9.64Ev06 8.56E06

1.54E04 7.19E05 4.62E06 8.56E06 8.98E05

1.49E04 3.79E05 4.13E06 9.64E06 8.56E06

#6

1 2 3 4

4.77E+01 9.24E+01 6.86E+02 5.50E+02

1.50E06 1.29E06 1.18E06 2.73E08

4.92E08 3.02E07 3.81E07 1.76E10

1.31E05 6.11E07 8.19E06 8.26E07

4.94E02 4.94E02 4.94E02 4.94E02

1.16E05 2.20E06 6.62E06 8.54E07

1.10E05 6.40E06 6.57E06 7.06E06

1.16E05 2.20E06 6.62E06 8.54E07

9

Y. Wei et al. / Computers and Structures 183 (2017) 1–13 Table 2 (continued) DOE sample

Design variables

Value of design variables

waA

waB

waC

w1

waAVM

waFDM

waDDM

5

5.11E+01

1.64E04

1.11E04

4.10E04

4.94E02

3.56E04

1.22E04

3.56E04

#7

1 2 3 4 5

4.92E+01 9.53E+01 7.71E+02 2.98E+02 4.74E+01

9.96E06 6.76E06 1.04E06 5.85E08 1.95E04

1.57E06 5.05E06 2.59E07 2.06E07 8.98E05

9.11E05 7.20E06 7.60E06 4.08E06 2.68E04

5.00E02 5.00E02 5.00E02 5.00E02 5.00E02

8.28E05 5.48E06 6.30E06 4.22E06 1.63E04

8.06E05 3.18E05 6.56E06 3.62E06 8.01E05

8.28E05 5.48E06 6.30E06 4.22E06 1.63E04

#8

1 2 3 4 5

5.06E+01 9.39E+01 1.11E+03 6.00E+02 4.19E+01

1.43E06 2.26E05 3.49E07 6.00E08 2.86E04

1.19E07 2.39E05 2.18E07 3.78E08 8.34E05

3.19E05 4.11E05 3.71E06 9.46E07 2.74E04

5.31E02 5.31E02 5.31E02 5.31E02 5.31E02

3.03E05 4.24E05 3.14E06 1.04E06 7.14E05

3.23E05 8.34E05 3.71E06 9.50E07 2.07E04

3.03E05 4.24E05 3.14E06 1.04E06 7.14E05

Table 3 Comparison of computing time of DDM and AVM (CPU time/s). DOE sample

DDM Total time

Contact solving time

DDM algorithm time

Total time

AVM Contact solving time

AVM algorithm time

#1 #2 #3 #4 #5 #6 #7 #8

430.50 436.67 444.28 440.23 436.19 440.42 440.59 437.98

341.05 346.59 355.19 346.27 346.05 348.83 351.38 348.38

89.45 90.08 89.09 93.97 90.14 91.59 89.22 89.61

359.22 364.97 373.56 364.17 363.75 366.13 369.03 366.42

339.64 345.80 354.27 344.89 344.52 347.11 349.72 346.70

19.58 19.17 19.30 19.28 19.23 19.02 19.31 19.72

Sensitivity calculations were carried out for a 205/50R16 passenger radial car tire using both DDM and AVM. It was found that the two methods give identical results. However, the pure computational time of AVM (determined by subtracting the computation time required by the contact problem solver from the total CPU time required by the whole process) is only about 1/5 of the CPU time required by DDM for evaluating sensitivities (see Table 3). DDM and AVM require almost the same time to solve the contact problem. However, because of the advantage of AVM mentioned in Section 7, the pure computational time of present approach can reduce to 1/n of DDM, where n is the number of design variables. Further discussion on this aspect is provided in Appendix B. Sensitivity of cost function was evaluated with respect to the five shape variables included in the optimization problem. For that purpose, eight levels were evenly selected by Latin Hypercube Design method [3] in the range of W 1 2 ½40:6; 50:6, W 2 2 ½90:3; 95:3, R1 2 ½600; 1200, R2 2 ½248; 600, and R3 2 ½40; 53 covering the whole design range for all variables as shown in Table 1. According to orthogonal-maximin Latin Hypercube Design algorithm [4,5], a number of samples greater than 1.5 times of the number of design variables is enough. Hence, the eight samples selected in the present analysis can efficiently cover the whole design range for five design variables. The 2D and 3D models of the selected 205/55R16 passenger car radial tire are shown in Fig. 2 and 4, respectively. Table 2 presents the results of the sensitivity analysis carried out for eight design points generated with the design of experiments technique.    waA ¼ KT t0 F a , waB ¼ KT t0 F Ca , waC ¼ aT fg a g, w1 is the objective function, a a a a wAVM ¼ wA þ wB  wC is the sensitivity value determined with the adjoint variable method, waFDM is the sensitivity value determined with the finite differentiation method, and waDDM is the sensitivity value determined with the direct differentiation method. It can be seen that AVM and DDM obtained identical values sensitivities while the FDM results show a slight error that may be caused by the fact that, unlike AVM and DDM, FDM computes sensitivity by

% reduction of CPU time passing from DDM to AVM

78.11% 78.72% 78.34% 79.48% 78.66% 79.24% 78.35% 78.00%

directly changing design variables. In this study, the interval for design variable change was set as one thousand of the original variable value.

9. Summary and conclusions This paper described an efficient method for performing sensitivity analysis in shape optimization problems of complex nonlinear structures such as automotive tires. Since computations for each design sample are time consuming, reducing the number of design samples is very important. In shape optimization problems, two issues must be addressed: (i) sensitivity calculation and (ii) disposing of an efficient optimizer. A possible strategy may be to combine AVM and GBK. While the latter was already studied in the literature [4], this research developed a novel AVM including a new adjoint shape design sensitivity formulation for nonlinear structures subject to contact forces. Sensitivity equations of displacements with respect to shape variations were formulated using the virtual work principle. Then, sensitivity analysis of average shear strains in the tire belt area was carried out. It was found that the proposed approach is more accurate than traditional finite difference sensitivity computations and much faster than the direct differentiation method as shown in Table 3. In view of this, the present approach has full potentiality for becoming a useful framework for shape optimization and sensitivity analysis of highly nonlinear structures including contact forces.

Acknowledgements Projects 11672148, 51275265 and 5117528 were supported by the National Science Foundation of China. This study was partially supported by the Manufacture Française des Pneumatiques Michelin and AvH Foundation. Valuable comments from three reviewers are greatly appreciated.

10

Y. Wei et al. / Computers and Structures 183 (2017) 1–13

Appendix A Notation t

The objective function

1

sw

s

t eij dij sV b t
s e ðxÞ 0 0 xi ; s xi m k 0 s ui

 a t 0 ui  a t 0 ðui ;j Þ

0 V ia

sa Dsa B,D

X

xoi ; yoi xi ; yi

T u F,B R t

LM

sCc

t 0 kn ,s kn t t st t f

s S ,s S ,s S t 0 Sij t 0 t k ,0

t 0f k

q

t t 0K L; 0K N; 0K G



a t a t 0Q ; 0FF t a 0FF  a t 0F R t a 0u

waAVM waFDM waDDM

The initial configuration corresponding to the shape parameter The load parameter The principle strains The Kronecker delta The total volume of the design area The effective shear strain The Lagrange coordinates for the original shape and the shape design corresponding to parameters sa The number of independent degrees of freedom (DOFs) The number of design parameters Virtual displacement The sensitivity of the displacements when sa ¼ 0 The sensitivity of displacement gradients Shape design velocity field The design variables A small increment of a design variable The tire width and diameter The whole shape domain The coordinates of the ith arc center The outer surface point coordinates belonging to the ith arc Mapping function The parametric coordinates The coefficient matrices The vector of design variables The contact boundary term in Lagrange multiplier approach The Lagrange multipliers The traction-free boundary, the traction-given boundary, and the contacted boundary The stress field The loads which are in equilibrium in the original shape Identical to those in the conventional nonlinear finite element analysis for the original design shape The change rates of the real loads and the fictitious loads The change rate of the fictitious load The change rate of the real loads Displacement sensitivities The sensitivity value determined with the adjoint variable method The sensitivity value determined with the finite differentiation method The sensitivity value determined with the direct differentiation method

Fig. A2.1. 3D finite element model and boundary conditions of the square rod.

in Fig. B.1. The top of rod is subject to uniform pressure while the bottom is put in frictionless contact with a stationary rigid surface. The rod is divided into 10 cubic HEX8 elements with length of 10 mm (see Fig. B.1). For the sake of simplicity, the side length lc of the first element is defined as design variable, closely related to nodal forces of top element. We assume that all four sides displace inwards or outwards together in a length d (see Fig. B.1). The objective function is the effective shear strain - Eq. (2) - developed in the ninth element. The shape design velocity is 0.5, i.e. the coordinate xc change rate with respect to the lc design variable. One can calculate the design sensitivities waFDM , waDDM , waAVM and the corresponding CPU times listed in Table B.1. Similar to Table 2, AVM and DDM get identical sensitivity values while FDM shows a small error. For this simple case, AVM is only 22% faster than DDM. Hence, the advantage of using AVM becomes much more significant for large scale problems.

Appendix C. Derivation of strain field sensitivity Eq. (18) According to the nonlinear theory, the strain field can be expressed as follows:

e ¼

t 0 ij

t

s eij



t 0 ui ;j

þ t0 uj ;i þ t0 uk ;i t0 uk ;j

. 2

ðC:1Þ

 . ¼ ts ui ;j þ ts uj ;i þ ts uk ;i st uk ;j 2  . ¼ t0 ui ;j þ t0 uj ;i þ t0 uk ;i t0 uk ;j 2  .     þ t0 ðui ;j Þ a þ t0 ðuj ;i Þ a þ t0 ðuk ;i Þ a t0 uk ;j þ t0 uk ;i t0 ðuk ;j Þ a sa 2  .   þ t0 ðuk ;i Þ a t0 ðuk ;j Þ b sa sb 2 ðC:2Þ

Table B.1 Comparison of FDM, DDM, and AVM results for the square rod problem.

Appendix B. Calibration case In order to generalize the findings presented in this study, FDM, DDM, and AVM are applied to the square rod problem schematized

CPU time (s) Sensitivity of objective function

FDM

DDM

AVM

2.80 2.66E02

3.25 2.61E02

2.55 2.61E02

11

Y. Wei et al. / Computers and Structures 183 (2017) 1–13

where ts ui ;j ¼ t0 ui ;j þ t0 ðui ;j Þ

a



sa , and t0 ðui ;j Þ a

  ¼ @@sa ½ts ðui ;j Þ

s¼0

, and in

which dds ½ts ðui ;j Þ consists of two parts: one is due to the difference a of ts ui and t0 ui , namely t0 ui , another one comes from the difference of ts ðui ;j Þ and t0 ðui ;j Þ in the case of identical ts ui and t0 ui . For the virtual displacement ds ui , which is arbitrary and geo
including all metrically admissible for the design shape vector s shape variables, the corresponding virtual strain can be expressed as follows i. ¼ ds ðui ;j Þ þ ds ðuj ;i Þ þ ds ðuk ;i Þts ðuk ;j Þ þ ts ðuk ;i Þds ðuk ;j Þ 2  .   ¼ d0 ui ;j þ d0 uj ;i þ d0 uk ;i t0 uk ;j þ t0 uk ;i d0 uk ;j 2 þ d0 ðui ;j Þ a

In order to simplify the present derivation, only frictionless contact is considered here. In Eq. (D.2), Sst c denotes the contact element. Developing equation (D.2) with respect to the shape variation factor sa , we have

XZ 0Ve

e

þ sa

h

ds eij







h

. i ðd0 ui ;j þ d0 uj ;i þ d0 un ;i t0 un ;j þ d0 un ;j t0 un ;i Þ 2 t0 Sij d0 V e

XZ 0Ve

e

h ðd0 ui ;j þ d0 uj ;i þ d0 un ;i t0 un ;j

      þd0 un ;j t0 un ;i Þ 2 t0 DTijkl t0 ðuk ;l Þ a þ t0 ðul ;k Þ a þt0 um ;k t0 ðum ;l Þ a



þd0 ðuj ;i Þ a þ d0 ðuk ;i Þ a t0 uk ;j þ t0 ðuk ;i Þ a d0 uk ;j þ d0 uk ;i t0 ðuk ;j Þ a .  .      þt0 uk ;i d0 ðuk ;j Þ a sa 2 þ d0 ðuk ;i Þ a t0 ðuk ;j Þ b þ t0 ðuk ;i Þ a d0 ðuk ;j Þ b sa sb 2

ðC:3Þ By omitting the small higher order term

 .     d0 ðuk ;i Þ a t0 ðuk ;j Þ b þ t0 ðuk ;i Þ a d0 ðuk ;j Þ b sa sb 2

a

þt0 um ;l t0 ðum ;k Þ

i XZ =2 d0 V e þ sa

0V

e a

h   d0 ðui ;j Þ a þ d0 ðuj ;i Þ a

a



þ d0 ðun ;j Þt0 ðun ;i Þ þd0 ðun ;i Þ a t0 ðun ;j Þ i XZ  ½ððd0 ui ;j þ d0 uj ;i þd0 ðun ;j Þ a t0 ðun ;i Þ =2 t0 Sij d0 V e þ sa

þd0 ðun ;i Þt0 ðun ;j Þ

0Ve

e

which implies that the shape variation is small terms, the strain sensitivities can be rewritten as

þ

 .   ds eij ¼ d0 ui ;j þ d0 uj ;i þ d0 uk ;i t0 uk ;j þ t0 uk ;i d0 uk ;j 2 d0 ðui ;j Þ a a

a

a

d0 un ;i t0 un ;j

þ sa sb



þ sa sb

Introducing Eqs (18)–(20) into the equation of virtual work corresponding to the shape variation, it follows

Z

t s s Sij ds eij j0 J jd0 V

XZ

0V

s t k ds uk ðds S=d0 SÞd0 S þ t t

0S

Z

þ

t

t t s s qs f k ds uk j0 J jd0 V

þ sa sb sc

t

t

t Sst 0 c

½ðts un þ 0 0 g n Þd0 kn þ ðts us  t0 us Þd0 ks

¼

t Ssl 0 c

Z

½ðs un þ 0 g n Þd0 kn ðds S=d0 SÞd0 S t

t Ssl 0 c

0

ðD:1Þ

h

0Ve

e

t 0 Sij

þ t0 DTijkl

 þt0 um ;k t0 ðum ;l Þ a



a t 0 ðuk ;l Þ

þ t0 ðul ;k Þ

a

a

þ t0 ðum ;k Þ

þ þ

t 0 Se

XZ e

t Sst 0 c

e

t Sst 0 c

XZ



t 0 kn

t

0 un

þ t0 kn

a

a

t Sst 0 c

e

þ sa

t Sst 0 c



d0 kn ðt0 un þ 0 0 g n Þd0 Se XZ 0Ve

e a

d0 un t0 kn d0 Se þ sa

XZ

d0 un t0 kn 0 Sea d0 Se þ sa

XZ t Sst 0 c

d0 un t0 kn

b 0 Se d0 Se

d0 kn t0 una d0 Se

XZ t Sst 0 c

e a

d0 uk0 qt0 f k 0 V sa ;s d0 V e 

t Sst 0 c

e 

t Sst 0 c

e

XZ

d0 uk t0 tk 0 Sea d0 Se þ sa

XZ

þ sa sb

d0 uk0 qt0 f k d0 V e

0Ve

e

t Sst 0 c

e



þ t0 um ;k t0 umb ;l Þ0 V sc ;s d0 V e

XZ

d0 un t0 kn d0 Se þ

XZ e

þ sa

a

t t t b 0 ðun ;j ÞÞ0 Dijkl ð0 uk ;l

0 Se



þ t0 um ;k t0 ðum ;l Þ b Þd0 V e



0Ve

XZ

a

½ðd0 ðui ;j Þ a þ d0 ðun ;i Þt0 ðun ;j Þ

e

XZ

t 0 um ;l

. i h . sa 2  d0 ui ;j þ d0 uj ;i þ d0 uk ;i t0 uk ;j þ t0 uk ;i d0 uk ;j 2



þ

þ sa

      þ d0 ðui ;j Þ a þ d0 ðuj ;i Þ a þ d0 ðuk ;i Þ a t0 uk ;j þ t0 ðuk ;i Þ a d0 uk ;j þ d0 uk ;i t0 ðuk ;j Þ a . i  þt0 uk ;i d0 ðuk ;j Þ a sa 2 ð1 þ 0 V ia ;i sa Þd0 V e XZ XZ a t t ¼ 0 q0 f k ð1 þ 0 V ia ;i sa Þd0 uk d0 V e 0 t k ð1 þ 0 Se sa Þd0 uk d0 Se þ e

b t t t 0 ðun ;j ÞÞ0 Dijkl ð0 ðuk ;l Þ

d0 uk t0 t k d0 Se þ

e

The above equation can be rewritten as XZ

st 0 Se

e

þ ðts ut  t0 ut Þd0 kt ðds S=d0 SÞd0 S Z þ ðts kn d0 un þ ts ts d0 us þ ts t t d0 ut Þðds S=d0 SÞd0 S þ

XZ

½ðd0 ðui ;j Þ a þ d0 ðun ;i Þt0 ðun ;j Þ

XZ

þ d0 ðun ;i Þ

Z þ

a

e

ðs kn d0 un þ s ks d0 us þ s kt d0 ut Þðds S=d0 SÞd0 S t

t Sst 0 c

0V





0Ve

þ d0 ðun ;i Þ

Z

¼



a

i  þd0 ðun ;i Þ a at0 ðun ;j ÞÞt0 Sij 0 V sb ;s d0 V e

e

t

ðd0 ui ;j þ d0 un ;i t0 un ;j Þt0 DTijkl ðt0 ðuk ;l Þ

þt0 um ;k t0 ðum ;l Þ a Þþðd0 ðui ;j Þ a a þ d0 ðun ;i Þt0 ðun ;j Þ a a

Appendix D. Derivation of displacement gradient sensitivity Eq. (23)

Z

h

0Ve

e

a

þd0 ðuj ;i Þ þ d0 ðuk ;i Þ t0 uk ;j þ t0 ðuk ;i Þ d0 uk ;j þ d0 uk ;i t0 ðuk ;j Þ .  þt0 uk ;i d0 ðuk ;j Þ a sa 2 ðC:4Þ

þ

XZ

d0 un ;j t0 un ;i Þ=2Þt0 Sij 0 V ka ;k d0 V e

þ

X

sa sb

e



d0 kn ðt0 un þ 0 0 g n Þ0 Sea d0 Se Z t Sst 0 c





d0 kn t0 una 0 Seb d0 Se ðD:3Þ

0Ve

e 

sa d0 un ð1 þ 0 Seb sb Þd0 Se

   þ t0 una sa þ 0 0 g n d0 kn ð1 þ 0 Seb sb Þd0 Se

ðD:2Þ

Neglecting small higher order terms, the terms related to the stress field t0 Sij and the loads t0 tk , 0 qt0 f k which are in equilibrium in the original shape, Eq. (D.3) can be rewritten as

12

Y. Wei et al. / Computers and Structures 183 (2017) 1–13

XZ 0Ve

e

XZ

þ

a t 0 ðul ;k Þ Þ=2Þd0 V e

þ

 d0 uj ;i Þ=2Þt0 Dijkl ððt0 um ;k t0 ðum ;l Þ a

þ

XZ

þ t0 ðul ;k Þ Þ=2Þd0 V e þ

XZ

½ððd0 un ;i t0 un ;j

þ d0 un ;j t0 un ;i Þ=2Þt0 Dijkl ððt0 um ;k t0 ðum ;l Þ 

0Ve

e

¼

 ðd0 uk Þðt0 tk Þ0 Sea d0 Se t

0 Se

e


ð0 qt0 f k Þ0 V ia ;i d0 V e 

þ

XZ

XZ

ðd0 uk Þ

½ððd0 ui ;j þ d0 uj ;i

þ d0 un ;i t0 un ;j þ d0 un ;j t0 un ;i Þ=2Þt0 Skl ð0 V ka ;k Þd0 V e XZ     ½ððd0 ðui ;j Þ a þ d0 ðuj ;i Þ a þ d0 ðun ;i Þ a t0 ðun ;j Þ e

0V

þ d0 ðun ;j Þ þ

XZ e

þ

a



t Sst 0 c

XZ e

t t 0 ðun ;i ÞÞ=2Þ0 Sij d0 V e

þ

X Z e

t Sst 0 c

d0 un t0 kn d0 Se

0

s xi

t Sst 0 c a

XZ e

t Sst 0 c

d0 kn ðt0 un

s xi ðnÞ

ðD:4Þ

When the structure is modelled using HEX8 element, one can obtain the following weak form control equations

Nm ¼

1 m m 0 0 0 ð1 þ nm 10 n1 Þð1 þ n20 n2 Þð1 þ n30 n3 Þ 8

0 xj ðnÞ

¼

ðD:5Þ

8 X Nm ðn10 ; n20 ; n30 Þ0 xm j m¼1

t 0 uj ðnÞ

¼

8 X

ðD:6Þ N

m

ðn10 ; n20 ; n30 Þt0 um j

0 J 220

0 J 320

0 J 230

0 J 330

3 5

ðD:9Þ

ð1 þ nm n0 n n0 Þ

ðD:10Þ

n0 ¼1 n0 –j0

 X 8 @N m t m 0 J k0 k u j0 Jj m¼1 @nk0 0 i

a

0 xi

 m 8 0 J k0 j X @N

j0 Jj m¼1 @nk0

8 X

ðD:12Þ

Nm ðnÞ0 xm i þ

8 X Nm ðnÞ0 V m ia sa

8 X N m ðnÞ0 V m ia

m a 0 xi

¼ 0V m ia

ðD:14Þ

Considering Eqs. (D.6)–(D.15), the virtual work equation (4.18) and its finite element form can further be rewritten as

XZ

      ðd0 ui ;j þ d0 uj ;i Þ=2 t0 DTijkl ðt0 uka ;l þ t0 ul a ;k Þ=2 d0 V e

0Ve

!

þ

m 0 V ka

þ

XZ



e

0Ve

e

0Ve

XZ XZ

0Ve

e

     ðd0 ui ;j þ d0 uj ;i Þ=2 t0 DTijkl ðt0 um ;k t0 uma ;l þ t0 um ;l t0 uma ;k Þ=2 d0 V e

h

i     ðd0 un ;i t0 un ;j þ d0 un ;j t0 un ;i Þ=2 t0 DTijkl ðt0 uka ;l þ t0 ul a ;k Þ=2 d0 V e

h

i   ðd0 un ;i t0 un ;j þ d0 un ;j t0 un ;i Þ=2 t0 DTijkl ðt0 um ;k t0 uma ;l

X   þt0 um ;l t0 uma ;k Þ=2 d0 V e þ ¼

ðD:8Þ

ðD:13Þ

m¼1

m¼1

ðD:7Þ  m 8 0 J k0 j X @N d 0 u i ;j ¼ d0 um i j0 Jj m¼1 @nk0

 @ d0 ðui ;j Þ a ¼ ds ðui ;j Þ ¼ d0 ui ;k 0 V ka ;j @ sa s¼0 ! !   X n 8 8 @Nm 0 J k0 j X @N 0 J k0 k n m ¼ d0 ui 0V j0 Jj m¼1 @nk0 j0 Jj n¼1 @nk0 ka

¼

ðnÞ ¼

þ

!

ðD:11Þ

t0 –k0

The change rate of the geometrical coordinates of the mapping points with respect to the shape variation parameter sa for a generic point and nodal points of an element can be formulated as

e

  a a @ t t ðs ui ;j Þ ¼ t0 ui ;j  t0 ui ;k 0 V ka ;j 0 ðui ;j Þ , ¼ @ sa s¼0

where

8

m¼1

m¼1

 m 8 0 J k0 j X @N t m a ¼ u  j0 Jj m¼1 @nk0 0 i

0 J 310

¼ 0 xi þ 0 V ia sa

0

þ 0 g n Þ0 Se d0 Se 0

0 J 210

For the shape variation, the geometrical coordinates for the initial configuration corresponding to the shape variation parameter sa can be formulated as

a





2 0 J 110 7 7 ¼ 4 0 J 120 5 0 J 130

s0 –j0

d0 kn ðt0 una þ 0 g na Þd0 Se d0 un t0 kn 0 Sea d0 Se þ

3

1 8 3 mY X 0 0x nm B C m 1 i j0 Jj ¼ ei0 j0 k0 J i0 1 J j0 2 J k0 3 ¼ ei0 j0 k0 @ ð1 þ nn0 nn0 ÞA 8 0 n ¼1 m¼1 n0 –i0 0 10 1 p p 3 q q 3 8 8 X X Y n 0 0 x2 Y n x 0 0 B CB C j 3 k ð1 þ nps0 ns0 ÞA@ ð1 þ nqt0 nt0 ÞA @ 8 s0 ¼1 8 t0 ¼1 p¼1 q¼1

0Ve

e

@ 0 x3 @n10 @ 0 x3 @n20 @ 0 x3 @n30

0

0Ve

e

@ 0 x2 @n10 @ 0 x2 @n20 @ 0 x2 @n30

P Q 0 where indices i; j ; n0 ; . . ., except those related to and , obey the rules of index symbol including the summation convention, ei0 j0 k0 is the permutation symbol. The determinant of the Jacobian matrix can be expressed compactly as

a

þ t0 um ;l t0 ðum ;k Þ a Þ=2Þd0 V e XZ  þ ½d0 uk ;i t0 Sij t0 ðuk ;j Þ a d0 V e XZ

mY 8 3 X nm j 0 0 xi

¼

0 J j0 i

m¼1

0Ve

e

0 x1 @n 0 6 @ 0 x11 6 @n 4 20 @ 0 x1 @n30

Using index symbol, it can be expressed compactly as



þ t0 um ;l t0 ðum ;k Þ a Þ=2Þd0 V e

½ððd0 un ;i t0 un ;j þ d0 un ;j t0 un ;i Þ=2Þt0 Dijkl ððt0 ðuk ;l Þ

a

¼

½ððd0 ui ;j

a

0Ve

e

0 Jðn10 ; n20 ; n30 Þ

0Ve

e

þ

2@

a

½ððd0 ui ;j þ d0 uj ;i Þ=2Þt0 Dijkl ððt0 ðuk ;l Þ

e

XZ e

 þ

a

t 0 Se

ðd0 uk Þðt0 t k Þ0 Se d0 Se þ

XZ e

0Ve

e

0Ve

XZ

Z 0Ve

ðd0 un ;i Þt0 Sij

XZ e

0Ve



t a 0 un ;j

d0 V e

ðd0 uk Þð0 qt0 f k Þ0 V ia ;i d0 V e

½ðd0 ui ;j þ d0 uj ;i þ d0 un ;i t0 un ;j þ d0 un ;j t0 un ;i Þ=2t0 Sij 0 V ka ;k d0 V e ½ððd0 un ;i t0 un ;k 0 V ka ;j þ d0 un ;j t0 un ;k 0 V ka ;i Þ=2Þt0 Sij d0 V e

Y. Wei et al. / Computers and Structures 183 (2017) 1–13

þ þ þ þ

XZ e

0Ve

e

0Ve

e

0Ve

e

0Ve

XZ XZ XZ

½ððd0 ui ;k t0 V ka ;j þ d0 uj ;k t0 V ka ;i Þ=2Þt0 Sij d0 V e ½ððd0 un ;k 0 V ka ;i t0 un ;j þ d0 un ;k 0 V ka ;j t0 un ;i Þ=2Þt0 Sij d0 V e ½ððd0 ui ;j þ d0 uj ;i Þ=2Þt0 DTijkl ððt0 uk ;m 0 V ma ;l þ t0 ul ;m 0 V ma ;k Þ=2Þd0 V e ½ððd0 ui ;j þ d0 uj ;i Þ=2Þt0 DTijkl ððt0 um ;k t0 um ;n 0 V na ;l

þ t0 um ;l t0 um ;n 0 V na ;k Þ=2Þd0 V e þ

XZ e

0Ve

½ððd0 un ;i t0 un ;j

þ d0 un ;j t0 un ;i Þ=2Þt0 DTijkl ððt0 uk ;m 0 V ma ;l þ t0 ul ;m 0 V ma ;k Þ=2Þd0 V e XZ þ ½ððd0 un ;i t0 un ;j þ d0 un ;j t0 un ;i Þ=2Þt0 DTijkl ððt0 um ;k t0 um ;n 0 V na ;l e

þ

0Ve

t t 0 um ;l 0 um ;n 0 V na ;k Þ=2Þd0 V e

þ

XZ e





þ d0 kn t0 una þ ðd0 un t0 kn Þ0 Sea Þd0 Se

t Sst 0 c

ðd0 un t0 kn

a



þ d0 kn 0 0 g na ðD:15Þ

Finally, sensitivity equations for state variables can be put in the matrix form of Eq. (24).

References [1] van Keulen F, Haftka RT, Kim NH. Review of options for structural design sensitivity analysis, 1: linear systems. Comput Methods Appl Mech Eng 2005;194(30–33):3213–43. [2] Wei Y, Rezgui A, Yao Z, Wang P. A comparative analysis of contact algorithms in contact shape optimization problems. Optimiz Eng 2012;13(4):595–623. [3] Xuan Y, Xiang J, Zhang W, Zhang Y. Gradient-based Kriging approximate model and its application research to optimization design. Sci China E 2009;52:1117–24. [4] Yao Z, Wei Y. Some ideas and progress on the shape optimization of nonlinear structures. Proc Eng 2012;31:600–7. [5] Joesph VR, Hung Y. Orthogonal-maximin Latin hypercube designs. Statist Sin 2008;18:171–86. [6] Choi KK, Kim NH. Structural sensitivity analysis and optimization, 1: linear systems. Berlin: Springer; 2005. [7] Choi KK, Kim NH. Structural sensitivity analysis and optimization, 2: nonlinear systems and applications. Berlin: Springer; 2005. [8] Burczynski T, Kane JH, Balakrishna C. Shape design sensitivity analysis via material derivative-adjoint variable technique for 3-D and 2-D curved boundary elements. Int J Numer Meth Eng 1995;38(17):2839–66. [9] Burczynski T, Kane JH, Balakrishna C. Comparison of shape design sensitivity analysis formulations via material derivative-adjoint variable and implicit differentiation techniques for 3-D and 2-D curved boundary element. Comput Methods Appl Mech Eng 1997;142(1–2):89–109. [10] Cho S, Choi KK. Design sensitivity analysis and optimization of non-linear transient dynamics, I: sizing design. Int J Numer Meth Eng 2000;48(3):351–73. [11] Choi KK, Duan W. Design sensitivity analysis and shape optimization of structural components with hyperelastic material. Comput Methods Appl Mech Eng 2000;187(1–2):219–43. [12] Grindeanu I, Chang KH, Choi KK, et al. Design sensitivity analysis of hyperelastic structures using a meshless method. AIAA J 1998;36(4):618–27.

13

[13] Hou GJW, Sheen JS, Chuang CH. Shape-sensitivity analysis and design optimization of linear, thermoelastic solids. AIAA J 1992;30(2):528–37. [14] Meric RA. Shape design sensitivity analysis and optimization for nonlinear heat and electric conduction problems. Numer Heat Transfer Part A-Appl 1998;34(2):185–203. [15] Park IH, Coulomb JL, Hahn SY. Design sensitivity analysis for nonlinear magnetostatic problems by continuum approach. J Phys III 1992;2 (11):2045–53. [16] Park IH, Lee HB, Kwak IG, et al. Design sensitivity analysis for steady-state eddy-current problems by continuum approach. IEEE Trans Magn 1994;30 (5):3411–4. [17] Park JS, Choi KK. Design sensitivity analysis and optimization of nonlinear structural systems with critical loads. J Mech Des 1992;114(2):305–12. [18] Ryu JS, Yao YY, Koh CS, et al. Optimal shape design of 3-D nonlinear electromagnetic devices using parameterized design sensitivity analysis. IEEE Trans Magn 2005;41(5):1792–5. [19] Santos JLT, Choi KK. Shape design sensitivity analysis of nonlinear structural systems. Struct Optim 1992;4(1):23–35. [20] Smith DE, Tortorelli DA, Tucker CL. Optimal design for polymer extrusion, II: sensitivity analysis for weakly-coupled nonlinear steady-state systems. Comput Methods Appl Mech Eng 1998;167(3–4):303–23. [21] Moon MY, Kim JH, Ha YD, et al. Adjoint design sensitivity analysis of dynamic crack propagation using peridynamic theory. Struct Multidisc Optimiz 2015;51(3):585–98. [22] Jang HL, Kim JH, Park Y, et al. Adjoint design sensitivity analysis of molecular dynamics in parallel computing environment. Int J Mech Mater Des 2014;10 (4):379–94. [23] Jang HL, Cho S. Adjoint shape design sensitivity analysis of fluid-solid interactions using concurrent mesh velocity in ALE formulation. Finite Elem Anal Des 2014;85:20–32. [24] Kim MG, Jang H, Kim H, et al. Multiscale adjoint design sensitivity analysis of atomistic-continuum dynamic systems using bridging scale decomposition. Modell Simul Mater Sci Eng 2013;21(3):035005 (25pp). [25] Matsumoto T, Yamada T, Takahashi T, et al. Acoustic design shape and topology sensitivity formulations based on adjoint method and BEM. CMESComp Model Eng Sci 2011;78(2):77–94. [26] Igarashi H, Watanabe K. Complex adjoint variable method for finite-element analysis of eddy current problems. IEEE Trans Magn 2010;46(8):2739–42. [27] Lee TH. Adjoint method for design sensitivity analysis of multiple eigenvalues and associated eigenvectors. AIAA J 2007;45(8):1998–2004. [28] Thomas JP, Hall KC, Dowell EH. Discrete adjoint approach for modeling unsteady aerodynamic design sensitivities. AIAA J 2005;43(9):1931–6. [29] Martins J, Alonso JJ, Reuther JJ. A coupled-adjoint sensitivity analysis method for high-fidelity aero-structural design. Optimiz Eng 2005;6(1):33–62. [30] Maute K, Nikbay M, Farhat C. Sensitivity analysis and design optimization of three-dimensional non-linear aeroelastic systems by the adjoint method. Int J Numer Meth Eng 2003;56(6):911–33. [31] Georgieva NK, Glavic S, Bakr MH, et al. Feasible adjoint sensitivity technique for EM design optimization. IEEE Trans Microw Theory Tech 2002;50 (12):2751–8. [32] Lee TH. An adjoint variable method for design sensitivity analysis of elastoplastic structures. KSME Int J 1999;13(3):246–52. [33] Park IH. Design sensitivity analysis for transient eddy current problems using finite element discretization and adjoint variable. IEEE Trans Magn 1996;32 (3):1242–5. [34] Belegundu AD, Arora JS. A sensitivity interpretation of adjoint variables in optimal-design. Comput Methods Appl Mech Eng 1985;48(1):81–9. [35] Guz AN. Fundamentals of the three-dimensional theory of stability of deformable bodies. Berlin Heidelberg: Springer-Verlag; 1999.