The application of the indirect boundary element method to optimum shape design

Pergamonl PII: soo45-7!M!q%)4lo439-7

Compurers & Strucrures Vol. 65. No. 5, pp. 651-6668,1997 0 1997 Elswier Science Ltd. All rights reserved Printed in Great Britain 004s7949/97 s17.00 + 0.00

THE APF’LICATION OF THE INDIRECT BOUNDARY METHOD TO OPTIMUM SHAPE DESIGN

ELEMENT

Osama K. Bedair Civil Engineering Department,

Concordia University, 1455 de Maisonneuve Blvd West, Montreal, Quebec, H3G lM8, Canada (Received 28 June 1996)

Abstract-This paper presents a strategy for optimum shape design using the indirect boundary element method. The optimization problem is solved iteratively by transforming the implicit non-linear problem to a sequence of explicit linear sub-problems. Beginning with initial geometry, the boundary is modified until convergence to the minimum stress distribution is achieved. The optimization is applicable to any loading condition and is adaptable to any boundary element code. The paper also presents three formulations for computing the stress sensitivity using the indirect boundary element method. The optimization algorithm is applied to optimize the shape of 90” V-notch under tensile stress and combined shear and tensile stresses. The optimum shapes and stress distributions for wide ranges of ~,~/a,, ratios are determined. The monotonic convergence of the peak boundary stress reflects the high accuracy the indirect boundary element method provides in the computation of the stress sensitivities. In all cases, convergence is attained at the fourth iteration. It is also shown that minimizing the boundary stress is a process that can be described as a curvature re-distribution process. The curvature is decreased in the highly stressed segment and is increased in the low stressed segments. Finally, it is shown that the constant stress distribution along the entire moving boundary is not always the optimum solution. The portion of the boundary, over which the optimum stress is constant, is a function of the applied loading. For uniform tensile loading, err, the optimum stress distribution is constant along the entire moving boundary. The constant stress region along the moving boundary then decreases as the shear loading is applied. 0 1997 Elsevier Science Ltd.

1. INTRODUCTION

Stress concentration problems are one of the most serious problems that are encountered in many structures. They normally result engineering wherever the curvature of the boundary is locally high. This could occur for example, when two surfaces with different angles need to be connected. Figure 1 shows some examples of stress concentration regions, denoted by the segments joining points AB. Other examples of such regions are also encountered in machined fillet of a crank shaft, a welded fillet region between a (chord and column member of an off-shore oil rig, or the root region of the threaded end closure cap of a high pressure vessel. Because of the potentially serious consequences of the development of fatigue crack in any of the stress concentration regions, considerable cost and effort are often expended by designers to determine the optimum shape for each of these regions, i.e. the shape which minimizes the peak stress. Because of the non-linear nature of the optimization problem, the computation cost is very expensive and extensive research is underway to improve the accuracy and reduce the computation cost of the optimization process. The finite element method has long been the predominant method of analysis. Recently, the

boundary element method, because of its superior efficiency for many problems, has attracted the interest of many researchers. In this paper, a boundary-element based technique is developed for shape optimization of stress concentration regions. The optimization problem is solved iteratively by transforming the implicit non-linear problem to a sequence of explicit linear subproblems. Beginning with initial geometry, the boundary is modified until convergence to the minimum stress distribution is achieved. The optimization algorithm is applicable to any loading conditions and is adaptable to any boundary element code. The optimization algorithm is applied to optimize the shape of a 90” V-notch. The optimum shapes for uniform tensile and wide range of combined tensile and shear stresses, crY,./en, are presented. 2. LITERATUREREVIEW Shape optimization has been a popular area of research in the last few years. Extensive literature has been published on shape optimization where the objective is to minimize the weight of the structure. A limited amount of literature is available where the objective is to minimize the maximum stress of the 651

0. K. Bedair

652

structure. This section represents a literature survey and review of the state of the art for optimum shape design of stress concentration problems. Because of the variety of methods in the treatment of the problem, the work done is classified into two major and experimental categories, namely analytical techniques and numerical techniques. 2.1. Analytical

and experimental

techniques

A constant stress distribution around the boundary of the stress concentration regions serves as the optimality criteria for the researchers using various analytical techniques. They use analytical procedures to find the geometry of the structure that correspond to this stress distribution. Cherepanov [l] treated the problem as an inverse plane elasticity problem. He derived analytical expressions to find, in general, the shape that corresponds to any prescribed stress distribution. By using the complex variable techniques of Mushelishvili and conformal mapping, he derived analytical expressions for infinite plates with one or many holes. Banichuk [2, 31 proved analytically, using the complex variables techniques, that a uniform tangential stress along the boundary of a hole in infinite plate loaded axially in tension is the optimum solution. Bjorkman and Richards [4] characterized the holes that do not perturb the Laplacian component of the tensor field in which they are placed as harmonic holes. Such holes are considered to be optimum. They employed complex variable techniques to solve the inverse problem and obtained the functional equations for these holes. They applied their formulations to a biaxially loaded plate and found that an ellipse with a major axis oriented parallel to the direction of the major principal axis is the harmonic hole of this problem. Durelli and Rajaiah [5] employed photoelastic techniques to achieve the optimum hole configuration in a finite plate loaded uniaxially in tension. The optimization process involves removal of material from lower stress regions around the hole by careful hand filing till an isochromatic fringe coincides with the boundary in the tensile and compression regions. The resulting stress distribution is uniform along the tensile and compression segment of the boundary. They concluded that the amount of reduction of stress concentration factors is proportional to the hole diameter to the plate width ratio. In other work [6], they applied the same optimality criteria to optimize the shape of the hole subject to uniaxial compressive stress. 2.2. Numerical

techniques

Two different numerical techniques has been used for minimizing design of stress concentration regions; the gradient-based techniques and non-gradientbased techniques. In the former, the stress gradients are used as a guide (or basis) to predict the geometric changes by which the maximum stress is reduced. The

search process is conducted as a part of the numerical optimization algorithm. The evaluation of the gradients, however, is processed outside the algorithm. For the non-gradient-based techniques, the constant stress distribution serves as the optimality criteria. Certain iterative procedures are often used to find the geometry that corresponds to the uniform stress distribution. The main advantage in the non-gradient-based techniques is that the stress gradients are not required to predict the geometric changes. Francavilla et al. [7] were amongst the first to approach the problem numerically. They approximated the moving boundary by a polynomial and treated its coefficients as the design variables. They employed the finite element method to evaluate the stress derivatives and employed the penalty function algorithm combined with a variable metric unconstrained minimization technique to get the numerical solution. They applied their formulation to optimize the shape of a fillet of a tension strip and to find the optimum configuration of a connecting rod. A second-order polynomial was used for the first example, while a fifth-order polynomial was used for the second example. Kristensen and Madsen [8] used different boundary representations to approximate the moveable boundary. They used a finite element model to evaluate the derivatives and modified the mesh at every iteration to improve the accuracy in the derivatives. They applied their algorithm to optimize a hole in an infinite plate under biaxial tension (aYT/crY_ = l/2) and to optimize the edge shape of a junction in a web frame under biaxial load or bending. Pedersen and Laursen [9] extended the work of Kristensen and Madsen [8] to three-dimensional axi-symmetric bodies. They applied their algorithm to optimize the shape of a shoulder fillet in a stepped bar under tension and bending, respectively. Middleton and Owen [lo] formulated the problem of optimizing axi-symmetric pressure vessels with torisphorical head intersected by a nozzle. Their optimization criterion is to minimize the maximum shearing stress occurring in the structure with prescribed geometrical constraints. They employed the penalty function procedure combined with a variable metric unconstrained minimization technique and the finite element method to compute the stress derivatives. Tvergaard [ 1l] formulated the problem of optimizing the shape of a fillet under pure tension or bending. To prevent sharp corners which result in an infinite local stress, he defined the boundary by a family of shape functions and treated the coefficients of these functions as the design variables. He used the variational analysis of stress field equations to evaluate the stress derivatives. Trompette and Marcelin [12, 131 minimized some integral stress functions that give rise to a uniform stress distribution along the moving boundary. They employed finite element methods to evaluate the stress gradients. They applied their formulation to

Application of the indirect boundary element method optimize the shape of a thrust bearing of a hydraulic jack loaded symmetrically and to optimize the shape of a shoulder fillet in a stepped bar loaded by pure bending. Schnack [14] used the constant stress distribution as an optimality criterion to develop a numerical procedure to find the shape associated with this distribution. He used the fundamental notch theory of Baud and Neuber as a basis to determine the feasible direction of the moving boundary. He employed the finite element method to evaluate the stress fields of the structure and applied his numerical procedure to optimize a hole in infinite plate loaded in tension biaxially (for o.V.y/oJJ = 1.5 and a,,/ a,, = 1.O) and to find the optimum transition in a fish plate. In other work, Schnack [15, 161 applied a similar procedure t.o three-dimensional axisymmetric structures. Demz and Morz [17] presented another iterative procedure for optimum stress distribution. They used a boundary perturbation analysis to derive optimality criteria and the finite element method to determine the optimum boundaries. However, their formulation is only valid for unloaded moving boundaries. In the area of shape sensitivity analysis, various numerical and analytical procedure using various methods of analysis was presented by Haug and Choi [18, 191 and Haug et al. [20]. 2.3. The application of boundary element for shape optimization

The boundary element method was first employed for shape optimization in the mid-80s. The method has proven to be more efficient and robust than its counterpart finite element method for many classes of problems. This is mainly due to the more accurate sensitivity information the method provides. Moreover, because discretization is confined only along the surface, the data processing is much simpler than the finite element method. In this section, a literature review of the various applications of the boundary element method in the context of shape optimization is presented. Mota Soares et al. [21] developed models for the shape optimal design of solids and hollow shafts. Their design objective is to maximize the torsional stiffness of the shaft for a given area. They compared the boundary element model with an equivalent finite element model and found that the former converges more rapidly. Zochowski and Mizukami [22] developed models for s:hape optimization of two-dimensional elasticity problems. They compared the boundary element model with an equivalent finite element model and reported that the boundary element model is more accurate, but less efficient in computational time. Burczynski and Adamczyk [23] developed models for two- and three-dimensional elasticity problems. Their design objective was to maximize the stiffness subject to a constant volume. They found that fewer analyses were required compared to the finite element method for the class

653

of problems they considered. Eizadian and Trompette [24] developed models for two-dimensional structures. They used substructuring techniques to represent the fixed and moving boundaries. Several applications were presented, including optimal shape design of a connecting rod and rotor. Numerical instabilities were reported in their studies. Defourny [25] employed boundary element and sequential quadratic programming for optimal design in potential problems. Mota Soares and Choi [26] presented models for two-dimensional elasticity and presented examples of minimizing the weight of a fillet. The boundary element method was also applied by Merit [27] to the optimal shape design of heat transfer problems. In the area of sensitivity analysis of elastic continua, efforts to date have been directed towards the direct boundary element method (DBEM). Mota Soares et al. [28], Mota Soares and Choi [26], Burczynski and Adamczyk [29], Choi and Ewak [30], Zhao and Adey [31] employed the variational approach to develop the sensitivity formulation based on DBEM. An improved formulation using the variational approach and including the effect of body forces was proposed by Saigal et al. [32]. Their formulation allows the evaluation of displacement and traction sensitivities at discrete locations of the domain. The stress sensitivities are then evaluated by direct differentiation of the boundary surface recovery equation. Kane and Saigal[33] presented an implicit differentiation formulation for sensitivity analysis using the DBEM. Their approach relies on a direct differentiation of the system of equations with respect to the shape design variables and then evaluation of the differentiated fundamental solution analytically. To evaluate the differentiated kernels which include the product of shape functions, fundamental solutions and their derivatives, they used a numerical integration techniques. Barone and Yang [34], developed a direct differentiation approach of the basic boundary element integral equation followed by a discretization of this equation. Sensitivity formulation for axisymmetric structures was presented by Saigal et al. [32]. In all of the above literature, the direct boundary element method (DBEM) was used as the numerical procedure for the optimization process. In this paper an optimization algorithm, using the indirect boundary element method (IBEM), is presented. As a first stage, a brief review is given about the indirect boundary element method. Then the mathematical formulation for the problem is given and the solution strategy. Finally some examples are given showing the efficiency of the algorithm. 3. THE INDIRECT BOUNDARY

ELEMENT

METHOD

(IBEM)

In this section, a brief description of the indirect boundary element method (IBEM), as applied to linearly elastic plane problems, is given to introduce

654

0. K. Bedair (a)

W

N

f-l

2/

.

\

‘i\ 4 \LV’ I

\

Fig. 3. Illustration of the numerical procedure of the IBEM for the plate with cavity problem: (a) the discretized cavity; (b) the fictitious model.

(a)

Let Gij,k and Uik denote the stress, crc and the displacement uir respectively, due to a unit point load at q acting in the k direction. Then, for a system of forces applied on S*, the stress induced at the field point can be computed by superposition of the stress and displacement fields due to the forces applied over increment dS* as

CL __-_-___-_-___-_-_-_---_-_---.

ue(x) =

s s

Gij>(4, x)P, (q)dS*

(1)

U,>(q,xV’dq)dS*.

(2)

s

Fig. 1. Examples of stress concentration regions; (a) crank hook; (b) a section of cylindrical part of jet engine fuel control housing.

the notation that will be used later when developing the sensitivity equations. Other versions of the method can also be found elsewhere [35-371. The philosophy of the IBEM is to assume that the body of interest, Q, is embedded in an infinite domain, SY, as shown in Fig. 2. A fictitious contour, S*, is then introduced and assumed to be separated by an infinitesimal distance, t, from the real body. The stresses and displacements at the field point, x, inside the body due to a point load P(q) applied on S* at a source point, q, can be calculated using the known solution for a point load in the infinite plane.

p

qd+s*

9 m

r /

‘,l

/‘,/;’

yy

,‘/ //’

IL,/’

X

/

/

/

s s*

3.

u,(x) =

s

Using the relation ti = aji(x) n,(x) one obtains for the tractions on S,

&(X) =

c

JS

7’dqr x)P&WS*

(3)

where T,k is the traction, t,, on S due to the afore-mentioned unit point load at q. If the stress intensities Pk(q) are known, the stress and displacement inside S can be computed directly from eqns (1) and (2). In this case, however, Pk(q) is not known and needs to be evaluated from the prescribed traction or displacement information on S using eqns (2) and (3) in an inverse manner. As an illustration, suppose that the traction boundary conditions, t, and ty are prescribed everywhere on the boundary, S. Therefore, P, and P, are the solutions to the coupled integral equations

d

/’

x

.E

&(x) =

no

Fig. 2. Analysis of the stress boundary value problem by the indirect boundary element method.

[TXq, x)P.&) + T&j, s .V

t,.(x) =

[T,:.&, x)P.r(q) + T,&, ss

x)P,ldS*

(4)

x)P,.]dS*. (5)

Application of the indirect boundary element method

655

Fig. 6. Schematic of the disolacements involved in the convexity constraint of node j.

u:(x) = i j=l

W,.&,

xX(q)

+ U.&j,

xV’:.ldS*

(8)

s .P

Typical discretization Fig. 4. Schematic of a typical stress concentration region (upper); and the segment to be optimized (lower).

It is generally not possible to find the analytical solution of these equations. However, the integral may be evaluated approximately by discretizing the fictitious boundary S* into N elements with assumed traction distributirons P, and Py. For a prescribed boundary conditions on S, the above integral equations then reduce to the following algebraic equations

t:(x) = -f ,=I

t;.(x) = f j=l

Similarly, prescribed, obtained

s s

Px,x(q,xP&)

Mixed formulations involving displacement and traction boundary conditions are handled by selecting the appropriate equations from eqns (6) to (9). The resulting set of equations may be written in a general matrix form as b’ = [#]pj

where {b’} is a vector contains the prescribed tractions or displacements boundary conditions, [A’i] is the matrix of the corresponding influence coetkients, and (pi) is the vector of the fictitious tractions. For brevity, eqn (10) will be presented hereafter as

+ Tw(q, x)P;,ldS* (6)

b = Ax.

9’

K,.&, xE(q) + T,.,.(q,xP;.l.ldS*. (7)

p

if the displacements u, and U, are the following algebraic equations are

(11)

Solving the above system of equations yields traction intensities P,,P?.Once these intensities the stresses or displacements at known, interior/boundary points within the domain, R, be directly evaluated from

G(x) =

5

j= I

+

. +

G.n,.&, xPXqW*

S*’

s

XV’!.dS*

s

XV’+.*

G.uLr (4,

G.n,r(q,

S”

the are any can

[S

1

9

Fig. 5. Schematic of a boundary segment at the kth and k + 1 stage.

(10)

1

(12)

(13)

0. K. Bedair

656

Similarly, eqns (15) and (16) can be written in a matrix form as

G,w(q, x)R(qW

{u’} = [uB]{Pj}.

j/,dq,x)P+*

+

(14)

For convenience, eqns (17) and (18) will be presented hereafter as d = cx.

z&(x) =

+

2 j=l

U.v(q, x)P!&)dS*

s

Udq, x)P:dS*

P’

u.:.(x) = t

j= I

For traction-free boundaries, such as the cavity shown in Fig. 3(a) t, and tYof eqns (4) and (5) are zero. Therefore, it is convenient to express the stress field as a summation of two components. The first component rY, is the free-field stress; the second component, &, is the perturbation of the free field stress due to the presence of the cavity, i.e.

U,;.&, xV’!JqW*

Y’

UJq,

x)P;.dS*

1.

(16)

Equations (12)-(14) can be written in matrix form as

{&} = [G"]{Pq.

(17)

where crfi is obtained by solving the problem without the cavity. The magnitude of 0; is determined by means of the fictitious model depicted in Fig. 3(b), which has identical discretization and is assumed to be at infinitesimal distance, E, from the real cavity S. Each element on the model has a prescribed shear and normal traction distributions Pi and pi of unknown intensities. In order to satisfy the boundary conditions on S, every element on the fictitious boundary S* should have a resultant shear and

Zero entries -

(19)

[S

s9’

+

(15)

1

(18)

1 0

*** **t **** l l

‘*-FIT

____I*

q

I*** I***

*

*, *,

l

*,-

**c************,**

**k***********,** l *it***********,** l *k_t&h********l+ p r r ry----

I

1

4 0

Active region

0

--

*1 * I l

l

0

I

Zero entries AA

AC Zero entries

Ab

Zero entries Fig. 7. Properties of AA, Ab, AC due to a perturbation of node j.

Application of the indirect boundary element method

I

6.5 Input data 6.0

E0

r

651

l LP ABE

Sensitivity processor

Optimizer

3.51 18

1

20

I

I

22

24

I

I

26 28 Element

I

I

30

32

I 34

I

36

Fig. 10. LP vs BE prediction after the first iteration for b,, = 1.

Analysis

Yes


Fig. 8. Flow chart of the optimization algorithm. normal tractions equal and opposite to the free field stress. In other words, by substituting t; = -aI and t; = -0: in eqns (6) and (7), and solving the system of eqn (lo), Pi and P:’ are determined. The magnitude of 0: at any interior/boundary points can then be calculated from eqns (12 j( 14). 4. FORMULATION OF THE OPTIMIZATION PROBLEM

In optimum shape design of stress concentration problems, the major concern is the region of the structure containing a stress concentration. The aim is to find the “best” shape for this part structure by varying certain shape parameters in order to find the best stress distribution. As will be shown later, the

optimization process can be described as a stress re-distribution process. By modifying the curvature of this region, stress from the high stress region is “shifted” to the surrounding low stress regions. It will be also shown that modifying this region affects only the local stress field near to the part of the boundary being modified. Consider a typical stress concentration region shown in Fig. 4, which shows two fixed boundaries I-, and I2 separated by angle B, and a transition region between A and B. This transition region AB could represent any of the stress concentration regions shown in Fig. 1 or it may also represent the welded fillet region between a chord and a column member of an off-shore oil rig, the machined fillet of a crank shaft or turbine blade, or the root region of the threaded end closure cap of a high pressure vessel. The design problem can be stated as follows. What is the shape of transition connecting points A and B which minimizes the maximum stress. Optimizing the shape of the segment AB becomes the problem of determining a unique set of locations of these N variable nodes which minimizes the maximum stress umarat any of M stress points along the boundary, i.e. min{maxcr,}

i= 1,2,. . .,M

(21)

where M is the total number of the stress points. In order to control the stress along the whole boundary, it is convenient to recast this problem as the minimization of a fictitious stress Q that must be less than or equal the original peak stress and is greater than or equal to the boundary stress at each of the 1 < i < A4 preselected locations which may, but need not, coincide with the N nodes, i.e. minimize 0 subject to 2N + 1 constraints a-Ui>O

Fig. 9. Geometric tdetails of the non-optimized notch, (t/R = 6, ye = YR= 45”).

6, -

UL

2 0

(22)

0. K. Bedair

658 30

Iteration

20

I 20

I 21

I 22

I

23

I

I

24

25

I

26

I

21

I

I

28

29

I

30

I

31

I

32

I

I

33

34

Elements Fig. 11. Stress change each iteration for G = 1

where cr,.is any specified lower bound and omaris the original maximum boundary stress. Since we do not known the location (or locations) of the maximum stress emaX,the fictitious stress o is treated as an unknown in additional to the design variables. Note that changing the loading boundary conditions changes only the upper limit of the fictitious stress d in the above formulation. The corresponding change in the non-linearity of the stress function affects the convergence of the optimization process. 4.1. Design variables Figure 5 depicts, a typical three noded segment of the discretized boundary. The open and solid circles depict the respective locations of nodes j - 1, j, j + 1 at the kth and k + 1st iterative stages. The S;- ,, ST and S:+ , denote the movements of these nodes from their positions at the kth stage. The 6.0-

movement of each node can be made independently, by taking the individual nodal coordinates as the design variables, or be dictated by a shape function, e.g. spline function, that controls the movement of all (or part of) the boundary. In the second approach, the moving boundary is approximated by one or a set of functions. Preselected parameter(s) of this or these functions are taken as the design variables. Increasing the number of design variables gives a better control of the curvature of the moving boundary. The major advantage of this approach is that the number of variables is not very large compared to the number of nodal coordinates. In addition, the smoothness of the boundary is better controlled. However, tangency and continuity, bounds etc., constraints are sometimes required. The major drawback of this approach is that the non-linearity of the stress function depends on the type of the selected shape function. Convergence of the optimization process is generally slow. In 6.0 r

I z E m 5.0i? d s 4.5 $ I

.LP

ABE

I I I I

.LP

ABE

4.0 4 3.5

I

I

I

19 20 21

I

I

I I I I I I I I I I I I

22 23 24 25 26 27 28 29 30 31

32 33 34 35

Elements Fig. 12. LP v BE prediction after the second iteration for

&V= 1.

3.511

I9

1 I I 1 I 1 I I 1 I I I ’ I ’ I

20 21

22 23 24 25 26 21

28

29 30 31

32 33

34 35

Elements Fig. 13. LP vs BE prediction after the third iteration for u.w= 1.

Application of the indirect boundary element method

---

Initial

--

Second iteration

-

First iteration Third iteration

This can be achieved by imposing additional geometrical constraints which assure that the movement of each of 1 <O

-6*%kih

-0.2 I

0.2I

0I

0.4I

0.6I

X-Coordinate

Fig. 14. Optimum shape evolution for o,.~= 1.

addition, the resulting stress distribution is also not uniform. By treating the nodal points of the boundary element discretization as the design variables, the boundary is defined by the location of a finite number of nodes. The coordinates of these nodes are treated as the design variables. Therefore, the number of the design variables is larger than the previous case. In addition, infeasible shapes from manufacturing or stress points of view may result. However, better control of the curvature is achieved, since each node can move independently. Therefore, uniform stress distributions and rapid convergence can be achieved. In this study, the nodal coordinates of the BE mesh are taken as the design variables. Instead of having two design variables (x,~) to determine the nodal position, the nodal movements (S,) along the bisector line are treated as the design variable. The angle of the bisector line is calculated by averaging the angles of the two adjacent linear elements. 4.2. Convexity requirement In order to prevent infeasible or unsmooth geometries, the curvature of the boundary is required to remain convex during the optimization process. - - - First iteration -Third iteration

Initial

I

---

Second iteration

-0.2

I

I

I

0

0.2

I

0.4

0.6

X-Coordinate Fig. 15. Optimum boundary stress evolution for u, = 1.

(23)

where each Sj is inherently a non-linear function of movements Sj_ 1 and Sj+ 1. In addition to the smoothness of the moving boundary, the convexity constraints eliminate the necessity of imposing additional geometrical constraints that increase the size of the problem. For example, they prevent reversal of curvature of the moving boundary which is undesirable from a manufacturing point of view. Having determined the design variables of the problem brings us to the major difficulty in solving the optimization problem which lies in the implicit character of the stress constraints of eqn (22). In other word, these functions are not explicitly known in terms of the design variables which force us to perform some approximations to overcome this difficulty as will be discussed in the next section. 4.3. Sequential linear programming (SLP) As mentioned in the previous section, the implicit character of the stress constraints constitutes the major difficulty in solving the optimization problem. One of the attractive numerical optimization algorithms for such problem is sequential linear programming. The idea of sequential linear programming is to replace the original non-linear problem with a sequence of linear subproblems. The approximate subproblem is obtained by writing the linear series expansion for the non-linear functions (whether the cost or constraints functions). Since linear approximations of the functions are used, limits must be imposed on changes in the design variables to the region of applicability of the linearization and to prevent unbounded solutions. These limits are usually called move limits in optimization literature and expressed as AS% I’ < AS!1.< AS?” I

3kk+rI

659

(24)

where ASP, AS:” are the lower and upper limits in the kth iteration, respectively. It should be emphasised that the selection of proper move limits is usually one of the most difficult tasks in using SLP. Large move limits invalidate the linear approximation while small ones slow the convergence of the algorithm. The selection of these limits is highly dependent on the degree of non-linearity of the linearized functions and experience with SLP. Although there have been recommendations concerning the selection of move limits in

0. K. Bedair

660

Notch region

t

0

-1

1

X-Coordinate Fig. 16. Stress change along the boundary due to shape variation for G = 1.

literature, they may only be efficient for a certain class of problems and cannot be generalized. Returning to the optimization problem previously formulated eqns (22) and (23), the stress and convexity constraints can be linearized for each of the 1 < j < N nodes from their locations in the kth state as a;+’ = a: + /$ fJf;.s; ,=I

(25)

$ = A; + S:,_,.S;_, + S:,+,.S;+,.

(26)

In eqn (25) 0: = %:/a$, the sensitivity of the boundary at point i to movement of node j. As can be seen from Fig. 6 ATis the maximum distance node j could move if nodes j - 1 and j + 1 did not move, while Si_, and Sji + 1 reflect the sensitivity of $ to movements of nodes j - 1 and j + 1. The Sj of eqns (25) and (26) are now the design variables for the linearized problem. The explicit linear problem can now be solved by any of the -

well-established linear programming algorithms. Since the stress functions are not linear with respect to the design variables, the optimum solution of the subproblem may not correspond to the primal (non-linear) problem. Therefore, the process is repeated (with updating of the move limits) until no further improvement of the objective function is attained. In other words, when the change in the design variables becomes negligible, i.e. IS,1 < C, where 6 is any small number giving the tolerance for the change in the design variables Sj. 5. SENSITIVITY

FORMULATION

In order to find the stress derivatives or sensitivities required for the SLP problem as defined by eqn (29, three different strategies can be used using the indirect boundary element method (IBEM); the implicit differentiation approach (IDA) (sometimes this approach is referred to as the direct approach), finite difference approach (FDA) and first order approximation approach (FOA). The IDA is based

Original .........Optimized

-6.0 -

-

Original .........Optimized

-6.2 -6.4 e *

-6.6 -

-6.8 -7.0

9 'ii : "b VI

0

2-

4

IJ xx

6

8

Fig. 17. Comparison of the non-optimized and optimized u.~.~ Fig. 18. Comparison of the non-optimized and optimized a, along the centerline of the notch. along the centerline of the notch.

Application of the indirect boundary element method

the influence matrix C, and dx/asj are the fictitious traction sensitivities. The fictitious traction sensitivities are obtained by differentiating eqn (11) with respect to the shape design variable Sj to get

\;_Iti.p...<

661

A$=[$-$]

(28)

or

Influenced region

Fig. 19.Schematic of the influenced region around the notch due to shape variation process.

A;=jI

on direct different.iation of the system of equations that evaluate the traction sensitivities of the fictitious model. By solving the resulting system of equations, the fictitious traction sensitivities are determined. The stress or displacement sensitivities are then evaluated by direct substitution of the fictitious traction sensitivities in the differentiated stress or displacement equations. The differentiated fundamental solution is evaluated using a finite difference scheme. In the FDA, a simple finite difference scheme is used to evaluate the stress or displacement sensitivities by perturbing the BE mesh and solving N system of equations (where N is the number of design variables). In the P3A, a first-order approximation is performed on the system of equations that evaluate the fictitious traction sensitivities, leading to a single system of equations solved for N sets of boundary conditions, 5.1. Implicit dtferentiation approach (IDA) In this case, the stress or displacement sensitivities are determined by direct differentiation of eqn (19) to obtain ad

ac

ax

s,=&x+c&

(27)

where x and C, are obtained from the primary analysis; aC/as, is obtained by direct differentiation of

oxy’“xx 0 ..‘.... 0.18 ----0.6-1.2 ---

where P can be interpreted as a fictitious boundary condition applied to the unperturbed geometry. Solving the above system of linear equations gives the required fictitious traction sensitivities, ax/as,. The evaluation of aA/&,, ab/asj and aC/asj in eqns (21) and (22) can be performed analytically or by a finite difference scheme. With the former method, additional kernels need to be added to the main BE code. Moreover, numerical difficulties arise in the treatment of the singularities in the differentiated kernels. The additional kernels and the special treatment of the singularities require substantial programming effort and computer time. With a finite difference scheme, no additional kernels need to be added to the main BE code. In the case of traction-free boundaries, the free field stress determines the boundary condition for the fictitious model. Therefore, the boundary stress sensitivity is the sum of the sensitivities of both crBand cc with respect to the shape design variable s,, i.e.

where &$/as, is obtained from aaR tt

aa:dxj

as, -

axj dsj + z

aG dYj

Note that at interior points a+%, 0.36 SW lfg~d&,,

ds,

(31)

is zero.

5.2. Finite diJference approach (FDA) Alternatively, the stress or displacement sensitivities can be calculated by a simple finite difference approximation, i.e. ad Ad -_N-= as, Aq

I

I

0.5

1.0

X-Coordinate Fig. 20. Boundary stress distributions for u~,~/G = 0, 0.18, 0.36, 0.6, 1.2 and 1.8.

CjXj - CX

Asj

where x, is obtained by perturbing the BE mesh and solving the system of equation N times (where N is the number of design variables), i.e. b, = Ajxi.

(33)

0. K. Bedair

662

The major drawback of this approach is the computation cost involved in every perturbation. The system of eqn (33) needs to be generated and solved N time for each iteration. However, this approach was found to yield good accuracy with the proper selection of the finite difference interval. Moreover, it does not require access to, nor modification of the BE source code.

Direction 40-

of tl

uxy’“xx-m

5.3. First-order approximation (FOA) approach Any approximation that avoids the necessity of having to solve the N + 1 sets of equations defined by eqn (33) without loss of accuracy is highly desirable. To achieve this, it is convenient to note that each of eqn (33) can be rewritten as AA,)(x

(A +

+

AXj) = (b + Ab,)

(34)

where AA,, Axj and Ab, denote the respective changes in A, x and b of eqn (9) as a result of the perturbation of the shape design variable Sj. If the second order product, AAj Ax,, is neglected, a first-order approximation of Xj is given by A$ =

b;”

(35)

where b,* = (Ab, - AA,x). Thus, instead of generating and solving N + 1 sets of equations, one has the much less costly task of solving only a single set of equations for N + 1 sets of boundary conditions, i.e. A{x, Z,,

. ) &} = {b, bS, . . ) bff}.

The stress sensitivities before,

are then approximated

e N & = CjZj - CX ASj as, - As,

(36) as

(37)

This approach reduces the cost of every iteration by a factor of (1N + l)/(N + l), where 1 is a parameter that depends on the efficiency of the programme, but is nevertheless very much less than one. The disadvantage of this approach is that it is necessary to have access to the BE source code. It can be observed that neither eqn (29) nor eqn (35) requires the expensive computation of a new A matrix for every perturbation. Instead, we are able to use the same kernels for the original and perturbed geometries by adding fictitious boundary conditions to the fictitious model. Therefore, the P and b;” of eqns (29) and (35), respectively, can be interpreted as a fictitious boundary conditions on the unperturbed geometry. 5.4. Properties of AAj AC) and Abj In a typical BE analysis, the evaluation of the kernels and the solution of the system of equations normally occupies a significant amount of computational time in the whole operation. By applying a

l

Angle to centroid of element with urnax

0-

Fig. 2 1.Variation of the peak stress location as u.~,./u,..ratio increases. finite difference scheme, we are actually evaluating the changes in the kernels of A and C due to the perturbation of the original BE mesh. If the nodal coordinates are the design variables, by perturbing only the jth node, only the adjacent two elements change from their original positions. Therefore, only the kernels involving these elements need to be recalculated. The remaining kernels remain unchanged. Thus, there are non-zero entries in AA,, AC, and Ab, only in the few rows and columns shown in Fig. 7. Therefore, the computation can be optimized by careful programming to avoid the computation of these zero entries. To have a useful insight to the time saving by doing that, if we assume that the geometry being optimized is discretized into 100 boundary element, with 20 design variables. We further assume that 10 iterations are required to reach the optimum solution, we can avoid the computation of 39, 192 kernels of the AA, matrices, 19,396 kernels in each of the AC, matrices and 196 entries in each of the A$ vector in each iteration, a total of 11,756,800 kernels in the whole operation! The process of solving the system of equations can be optimized by saving the A matrix generated in the original analysis in triangular form, and computing the triangularization operation to the fictitious load vector every perturbation. 6. COMPUTER IMPLEMENTATION

The major components of the optimization process are depicted in Fig. 8. l The first step is to define the initial boundary geometry and the loading conditions. The boundary consists of the transition region whose shape will be subsequently be modified and the fixed region. l The second step is to analyze the geometry in order to determine the maximum stress value urnax.

Application of the indirect boundary element method -5.8 -

0 2 a $ -5.9 3 &

-0

A

--

of the indirect boundary element method (IBEM), the initially circular segment AB and each of the flanking sides were approximated by 15 and 19 linear elements, respectively. Nodes 21-34 were allowed to move and stresses at the centroids of elements 19-36 were required to be less than the unknown upper bound stress cr. Therefore the explicit LP problem is stated as; minimize d Subject to

Q,y%x

_ _ _ 0.06 ..... .. 0.12

0.18---

0.24--*0.36

\

-

1 pi

1

-6.0 -

-0.4

%-

1

I

I

-0.2

0

0.2

663

I

I

0.4

0.6

X-Coordinate

(38)

Fig. 22. Optimum shapes for u,~~/u,,= 0, 0.06, 0.12, 0.18, 0.24 and 0.36. where

The third step is to evaluate the stress sensitivities c1; through the sensitivity processor. Accuracy of the c>, depends on the method of analysis and is very important to converge rapidly. l The fourth step is to find the optimum solution of the approximate problem. This problem involves formulating the LP subproblem, proper selection of the move limits and solving the LP problem. Geometrical constraints are imposed at this stage. l The fifth step is to update the current geometry and reanalyse the structure to determine the reduced Qmal.

-CL, d c G &.

(39)

l

7. NUMERICAL EXAMPLES

Consider the problem of optimizing the shape of the segment AB at the root of a symmetric 90” V-notch (i.e. yL = yR = 45”) in the traction free edge of the semi-infinite plate shown in Fig. 9. All nodes along the moving boundary AB were required to satisfy the convexity requirements described previously. In this section the optimum shapes are determined for two loading conditions: (1) uniform tensile stress; (2) combined uniform tensile stress and shear stresses at infinity. For the analysis by means

-0.6

-0.4

I

I

-0.2 0 0.2 X-Coordinate

I

0.4

I

0.6

Fig. 23. Optimum stress distributions for ~.~.~/a.~.~ = 0, 0.06, 0.18, 0.24, 0.3 and 0.36. CAS 65/S-B

AS:+BS:_,-m(xik-x:_,)-y,k_,+y~~O A = sin yj - m

cos l/i

B=mCOSyj-sinyj_,

(40) (41)

(42)

(.Yf+,+S;+,Sinyj+~)-(yk+$_,sinyj_,) m = (xf- , + Sk ,+I COSYj-I)- (X:+1 + Sj-1 Sitlyj-1) (43)

where a:, c: are computed for the kth state and Sk denote the movement from the kth state to the k + 1st state. The centroid of element i (where the boundary stress and sensitivity coefficients aij are calculated) lies midway between node j = i and j = i + 1. Moreover, in order to validate the linear approximation, the nodal movements were limited by the following side constraints;

where AS,? and ASP are the lower and upper limits in the kth iteration, respectively. The move limit strategy was to allow for the largest possible nodal movements in the first iterations. Then they were reduced, when necessary, as the shape approaches the optimum. This was achieved, iteratively, by testing different values of move limits until no further change in the LP solution is obtained, or in other word, when the side constraints become redundant. Move limits of 0.1 unit satisfied this criteria for the first two iterations. They were then reduced to 0.05 units for subsequent iterations. I. 1. Uniform tensile stress

In this case, the plate was subjected to a uniaxial tensile stress of unit magnitude at infinity. The initial peak stress is 6.34 at the root of the notch (element

664

0.

K. Bedair

6.4r

I 5.41 01234567

I

I

I

I

I

I

I 5.41 01234567

Iteration 6.4

I

I

I

I

I

I

Iteration 6.5

r

r

0.18

. 5.41

I 01234567

I

I

I

I

I

I

Iteration

Iteration

Fig. 24. Convergence histories for un /OVV = 0, 0.06, 0.18 and 0.24.

27). This value is within 0.94% of the value given by

Peterson [38]. The LP and BE stress distributions after the first iteration, k = 1, are shown in Fig. 10. The LP prediction has a value of 5.85 over the entire region between elements 21 and 33. The actual stress distribution varies between 5.94 at each end and 5.8 at the root of the notch. As shown in Fig. 11, this decrease in the initially high stress region at elements 24-29 is achieved by increasing the initially low stress at elements 21-23. The LP and BE stress distributions after the second iteration, k = 2, are shown in Fig. 12. The LP has a constant stress value of 5.7 1 at elements 20-34. The BE stress distribution has a reduced peak stress value of 5.62 at elements 25-28. The constant stress value is decreased at this stage by, increasing the stresses at elements 20 and 34 in order to decrease the stresses at elements 21-33 as shown in Fig. 11. The large increase in the stress values at these two elements can be thought as being composed of two parts. The first part is what they gained in the first iteration. The second part is what they distribute to elements 2&34 to decrease their stresses. The third iteration, k = 3, produces only a slight change in the boundary stress. Stresses at elements 20, 21, 33, 33 increase slightly to expand the constant stress region. The peak stress is reduced to a constant value of 5.56 at elements 2&34 and the LP has a constant value of 5.57 at elements 20-34 as shown in Fig. 13. Subsequent iterations produced a negligible change in the boundary stress. Therefore, convergence has been achieved.

The shape evolution during the optimization process is shown in Fig. 14 and can be summarized as follow; -k=l l l l

Flattening in the central region. Increase in curvature of the outer regions. Almost no change in the depth of the notch.

-k=2

Increase in curvature of the outer regions. a Upward movement of the notch root.

l

Evolution of stress during the optimization is shown in Fig. 15. Note that widening of the constant stress region is largely due to the major increase in stress at the outer nodes. From the above description, one can conclude that the optimization philosophy for this class of problems requires that the stress at the centre part of the notch (where the stress is highest) be reduced while at the same time allowing the stress at the ends of the notch to increase as shown in Fig. 16. It interesting to note that the average non-optimized stress for elements 20-34 (over which the optimum value stress is constant) is 5.61 which is very close to the optimum value. It can also be observed from Fig. 16 that the change in the tangential stress along the traction-free boundary due to the shape variation

Application of the indirect boundary element method

-5.8

1.2 .--*c--.

1.8 .--*--5.9

2.4

2

2 ._ 2

-6

8

v > -6.1

-6.2

-0.6

-0.6

a4

-0.2

0

0.2

0.4

0.6

0.6

X-Coordinate Fig. 25. Optimum shapes for u.~~/u, = 0.6, 1.2, 1.8 and 2.4.

process confined ,within distance f R measured from the centerline of the notch. Stresses at greater distances remain almost unchanged from their original values during the optimization process. The changes of CT,.~ and a,, along the centerline of the notch are also depicted in Figs 17 and 18, respectively. A reduction in the stress values is attained up to distance R measured from the deepest point of the original notch root. Beyond that distance, stresses remain almost unperturbed. Based on the observations of Figs 16-18, one can conclude that optimizing th.e shape of the notch root affects the

stresses only in a region within a semi-circle of radius R, centred at the deepest point of the original notch

root as shown in Fig. 19. 1.2. Combined shear and tensile stresses In this case, the magnitudes of the symmetric and anti-symmetric components of the high gradient stress fields around the root of the notch were varied by means of the prescribed values of a uniform a,, and linearly increasing u., . The non-optimized boundary stress distributions for ~.~~/cr.~.~ = 0.00, 0.18, 0.36, 0.60, 1.20, 1.80 are shown in Fig. 20. As a.V,./o.r.V ratio increases, the location of the peak stress, a,,, moves to the right by an angle 0 measured from the

#J-5.9 3 T-60 8. u $ -6.1 -2 -0.8 -0.6 -0.1 -0.2

0

0.2

0.4

0.6

0.8

X-Coordinate Fig. 26. Optimum stress distributions for a.~~/~.~., = 0.6, 1.2, 1.8 and 2.4.

-0.8 -0.6 -0.4 -0.2

0

0.2

0.4

0.6

0.8

X-Coordinate Fig. 27. Optimum shapes for ~.~,/a,~.~ = 3 and 6.

0. K. Bedair

666 2015-

D-A-A-A&A-h-A

=xy’=xx -

3

P’

8. CONCLUSIONS

Shape optimization has received the attention of many researchers [l-26, 28-34, 39-631. Several numerical and analytical procedures have been presented in the literature. The paper presented an alternative numerical approach for optimum shape design using the indirect boundary element method. The philosophy of the method is to transform the implicit non-linear problem into a sequence of explicit linear sub-problems. The explicit subproblem I I I I I I I -20L 1 can then be solved by any linear programming -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 technique. Beginning with the initial geometry, the X-Coordinate boundary is modified until convergence to the Fig. 28. Optimum stress distributions for o~,./u.~~ = 3 and 6. minimum stress distribution is achieved. The optimization algorithm is applicable to any loading conditions and is adaptable to any existing BE-Code. centerline of the notch as shown in Fig. 21. Note that The various sensitivity formulations using the the direction of 0 depends on the sign of the a.r,./o,r.r indirect boundary element method has been preratio. The peak stress values for these ratios are 6.34 sented. A major advantage of these formulations over (6’= O.O’), 6.344 (0 = 4.8”), 6.52 (0 = 14.2”), 6.95 other boundary element based sensitivity analysis is (6 = 19.0”), 8.37 (0 = 24.0”) and 10.02 (0 = 28.5% that the stress sensitivities are not obtained by respectively. For cr,J/u,, ratios of 2.4, 3.0 and 6.0 the differentiating the approximated displacement field location of the peak stress is fixed at 0 = 33.5”. The twice as in papers by Kane and Saigal[33] and peak stress values for these ratios are 11.8, 13.6 and Barone and Yang [34]. Moreover, the derivation of 22.5, respectively. The location of the peak stress higher order kernels as done by Kane and Saigal[33] remains unchanged at 33.5” for any higher Q/O, and Barone and Yang [34] and the associated special ratio. integration techniques which require a substantial The optimized boundaries and stress distributions programming effort and expensive computer time are for Q,/u.~.~= 0.06, 0.18, 0.24, 0.30, 0.36 are shown in avoided. Figs 22 and 23. The initial peak stresses are 6.34,6.34, Minimizing the boundary stress proved to be 6.38, 6.43 and 6.52, respectively. The optimized closely linked to a re-distribution of the curvature. values are 5.74, 5.75, 5.82, 5.92 and 6.02, respectively. The curvature is decreased in the originally most In each case the boundary shape and stress highly stressed segment and is increased in the distributions are skewed, with the regions of originally least stressed segments. By modifying the minimum curvature and maximum stress occurring curvature in this way, the stress is lowered at the on the side of the notch where the boundary stresses segment with low curvature and is increased at the from Q.~.~ and 0,. are additive. Note also that with segment with high curvature. It also has been shown increasing a,ry/a, values there is a slight decrease in that modifying the curvature at the stress concenthe portion of the notch boundary over which the tration region has a local inhuence on the structure. peak stress is maximum, and an increasing portion of Stresses within r < R contour around the moving the boundary which has reached the tangent to one boundary only changes due to the curvature of the flanking side of the V-notch. As u,/u,, re-distribution process. increases, the notch depth increases until it reaches The investigation has also shown that the constant the initial depth of the original (non-optimized) stress distribution along the entire moving boundary shape. This is achieved at u,/u, = 0.36. is not a permanent optimum solution as described The convergence history of some of these load by Chereponov [ 11, Banichuk [3], Durelli and cases is shown in Fig. 24. As can be seen that the Rajaiah [5], Trompette and Marcelin [12] and optimum solution is attained at the fourth iteration Demz and Morz [ 171. The portion of the boundary for all of these load cases, reflecting the high over which the optimum stress is constant is only efficiency of the optimization algorithm. Figures a function of the type and magnitude of the 25-28 show the optimum boundaries and stress applied load on the structure. For uniform tensile distributions for ~.,./a, = 0.6, 1.2, 1.8, 2.4, 3.0, stress, ulr, the optimum stress distribution is constant 6.0. The initial (non-optimized) peak stresses for along the entire moving boundary. The constant these ratios are 6.88, 8.37, 10.02, 11.80, 13.60 and stress region decreases as the magnitude of u.&,, 22.50, respectively. The optimized peak stresses are ratio increases. 6.25, 7.36, 8.50, 9.80, 11.26 and 18.20, respectively. Acknowledgemenrs-The financial support provided from As can be seen, the depth of the notch continues the Natural Sciences and Engineering Research Council to increase as the magnitude of ux~/uxx ratio (NSERC) of Canada is greatly acknowledged. The increases with a decrease in the constant stress discussion and comments of Dr J. C. Thompson about the subject are greatly appreciated. region.

Application of the indirect boundary element method REFERENCES

1. Chereponov, G. P., Inverse Problems Of The Plane Theory Of Elasticity, PMM, 1974, Vol. 38, No. 6,

pp. 963-979. 2. Banichuk, N. V., Optimality Conditions In The Problem Of Seeking The Hole Shapes In Elastic Bodies, Vol. 41. PMM, 1977, pp. 92&925. 3. Banichuk, N. V., Problems and Methocis of Optimal Structure Design. Plenum Press, New York, 1983. 4. Bjorkman, G. S. and Richards, R., Harmonic holes-an inverse problem in elasticity. Journal of Applied Mechanics,

1976, 43, 414418.

5. Durelli, A. J. and Rajaiah, K., Optimum hole shape in finite plates under uniaxial load. Journal of Applied Mechanics,

1979, 46, 691-695.

6. Durelli, A. J., Erickson, M. and Rajaiah, K., Optimum shapes of central holes in square plates subjected to uniaxial unifomi loads. International Journal of Solids Structures,

1981, 17, 787-792.

I. Francavilla, A., Ramakrishnan, C. V. and Zienkiewicz, 0. C., Optimization of shape to minimize stress concentration. Journal of Strain Analysis, 1975, 12, 63-70.

8. Kristensen, E. S. and Madsen, N. F., On the optimum shape of fillets in plates subjected to multiple in plane loading cases. International Journal of Numerical Methods in Eng!neering, 1976, 10 ,1007-1019. 9. Pedersen, P. and Laursen, C. L., Design for minimum stress concentration by finite elements and linear programming. Journal of Structure Mechanics, 1983,10, 375-39 1.

10. Middleton, J. and Owen, D. R. J., Automated design optimization to minimize shearing stress in axisymmetric pressure vesr.els. Nuclear Engineering Design, 1977, 44, 357-366. 11. Tvergaard, V., On the optimum shape of a fillet in a flat bar with restrictions. Proceedings of IUTAM Symposium on Optimization in Structural Design, eds A. Sawczuk and Z. Morz. Springer, Warsaw, 1975, pp. 181-195. 12. Trompette, P. and Marcelin, J. C., On the choice of objectives in shape optimization. Engineering Optics, 1987, 11, 89-102.

13. Trompette, P., Marcelin, J. L. and Lallemand, C., Optimal shape design of axisymmetric structures. In Automated Structural Design, 1986, pp. 283-295. 14. Schnack, E., An optimization procedure for stress concentrations by the finite element technique. International Journal of Numerical Methods in Engineering, 1979, 14, 115-124.

15. Schnack, E. and Inacu, G., Shape design of elastostatic structures based on local perturbation analysis. Structural Optics, 1989, 1, 117-125. 16. Schnack, E. and Inacu, G., Gradientless computer methods in shape designing. Proceedings Of Second International Conference On Computer Aided Optimum Design Of Stru(:ture, OPTI 91, Boston, USA, 1991,

pp. 363-376. 17. Dems, K. and Morz, Z., Multiparameter structural shape optimization by finite element method. International Journal of Numerical Methods in Engineering, 1978, 13, 247-263. 18. Haug, E. J. and Choi, K. K., A variational method for shape optimal design of elastic structures. In New Directions in Optimum Structural Design. Wiley, 1984. 19. Haug, E. J. and Choi, K. K., Material derivative methods for shape design sensitivity analysis. In Automated Structural Design. Plenum Press, 1986, pp. 29-55. 20. Haug, E. J., Choi, K. K. and Komkov, V., Design Sensitivity Anal;jsis Of Structural Systems. Academic Press, New Yor.k, 1986.

661

21. Mota Soares, C. A., Rodrigues, H. C., Oliveira, L. M. and Haug, E. J., Optimization of the geometry of shafts using boundary elements. Journal of Mechanism, Transmissions and Automation In Design, 1984a, 106, 199-202. 22. Zochowski, A. and Mizukami, K., A Comparison Of BEM and FEM in Minimum Weight De&m. in Boundary Elements. Springer, West Germany,- 1985, pp. 901-911. 23. Burczynski, T. and Adamczyk, T., The boundary element formulation for multiparameter structure shape optimization. Applied Mathematics Modelling, 1985a, 9, 195-200. 24. Eizadian, D. and Trompette, P., Shape optimization of bidimensional structures by boundary element method. Conference on CAD/Cam, Design, Tuscan, Arizona.

Robotics and Automation in

25. Defourny, M., Optimization techniques and boundary element method. In Boundary Elements, Vol. 3, Springer, Berlin, 1988, pp. 479490. 26. Mota Soares, C. A. and Choi, K. K., Boundary elements in shape optimal design of structures. Computer Aided Optimal Design, NATO/NASA/NSF/ USAF Advanced Study, Vol. 2, Portugal, 1986,

pp. 145-185. 27. Merit, R. A. Boundary element methods for static optimal heating of solids. Journal of Heat Transfer, 1984, 106, 876880. 28. Mota Soares, C. A., Rodrigues, H. C. and Choi, K. K.,

Shape optimal structural design using boundary elements and minimum compliance techniques. Journal of Mechanism, Transmissions and Automation In Design,

1984b, 106, 518-523. 29. Burczynski, T. and Adamczyk, T., The boundary element method for shape design synthesis of elastic structures. 7th International Conference On Boundary Element Methods. Springer, Berlin, 1985b, pp. 408421. 30. Choi, J. H. and Kwak, B. M., Boundary integral equation method for shape design sensitivity analysis. In Computer Aided Optimal Design, Structural and Mechanical Systems, Vol. 2, ed. C. A. Mota Soares.

NATO/NASA/NSF/USAF, 1986, pp. 186-214. 31. Zhao, Z. and Adey, R. A., Shape optimization-a numerical consideration. Proceedings of the 11th International Conference on BEM, Cambridge, USA. 32. Saigal, S. J., Borggard, J. T. and Kane, J. H. Boundary element implicit differentiation equations for design sensitivities of axisymmetric structures. International Journal of Solids Structures, 1989, 25(5), 527-538. 33. Kane, J. H. and Saigal, S., Design sensitivity analysis of solids using BEM. Journal of Engineering Mechanics

ASCE, 1988, 114, 1703-1722. 34. Barone, M. R. and Yang, R. J., Boundary integral equations for recovery of design sensitivities in shape optimization. AIAA Journal, 1988, 26, 589594. 35. Banerjee, P. K. and Butterfield, R., Boundary Element Methods In Engineering Science. McGraw-Hill, U.K., 1981. 36. Cruse, T. A., Application of the boundary integral equation method to three-dimensional stress analysis. Journal of Computers and Structures, 1913, 3, 507-521. 37. Massonnet, C. E., Numerical use of integral procedures. In Stress Analysis, eds 0. C. Zienkiewicz and G. S.

Hoister. Wiley, New York, 1965, pp. 198-235. 38. Peterson. R. E., Stress Concentration Factors. Wilev. .., New York, 1974. 39. Bennett, J. A. and Botkin, M. E., Structural shape optimization with geometric description and adaptive mesh refinement. AIAA Journal, 1985, 33, 458-464. 40. Brabiant, V. and Fleury, Shape optimal design using B-splines. Computer Methods in Applied Mechanics And Engineering,

1984, 44, 247-261.

668

0. K. Bedair

41. Braibant,

V., Shape sensitivity by finite elements.

Journal of Structural Mechanics, 1986, 14, 209-228.

42. Cea,

J., Problems

Optimization Of

of shape

Distributed

optimal Parameter

design.

In

Structures.

Sijthoff Noordhoff, Netherlands, 1981, pp. 1981. 43. Choi, K. K. and Haug, E. J., Shape design sensitivity analysis of elastic structures. Journal of Structural Mechanics, 1983, 11, 231-269. 44. Choi, K. K. and Frederich, M. C., Implementation of design sensitivity analysis with existing finite element codes. ASME Journal of Mechanisms, Transmissions and Automation in Design, 1987, 109, 385-391. 45. Defourny, M., Boundary element method and design optimization. In Boundary Elements, Vol. 9, Springer, Berlin, 1986, pp. 4633471. 46. Dems, K. and Morz, Z., Variational approach by means of adjoint systems to structural optimization and sensitivity analysis, II structural shape variation. International Journal of Solids and Structures, 1984, 20, 527-552.

41. Gulperi, A. and Mukherjee, S., Design sensitivity analysis for potential problems by the derivative boundary element method. In Boundary Element Methods in Engineering, 1989, pp. 412-419. 48. Haftka, R. T. and Grandhi, R. V., Structural shape optimization-a survey. Computer Methods of Applied Mechanics in Engineering,

1986, 57, 91-106.

49. Haftka, R. T. and Adelman, H. M., Recent developments in structural sensitivity analysis. In Structural Optimization 1989, pp. 137-151. 50. Haftka, R. T. and Mroz, Z., First- and second-order sensitivity derivatives of linear and non-linear structures. AIAA Journal, 1986, 24, 1187-1192. 51. Immam, M. H., Three-dimensional shape optimization. International Journal of Numerical Methods in Engineering, 1982, 18, 661-673.

52. Miranda, C. and El-Yafi, F., 3-D optimum design using BEM technique. In Boundary Element Methods, Vol. 3. Springer, Berlin, 1988, pp. 449-462.

53. Pironneau, O., Optimal Shape Design for Elliptical Systems. Springer, Berlin, 1984. 54. Queau, J. P. and Trompette, P., Two-dimensional shape optimal design by finite element method. International Journal of Numerical Methods in Engineering, 1980, 15,

1603-1612. 55. Ramakrishnan, C. V. and Francavilla, A., Structural shape optimization using penalty functions. Journal of Structural Mechanics, 1975, 3, 403422.

56. Richetts, R. E. and Zienkiewicz, 0. C., Shape optimization of continuum structures. In New Directions in Structural Design, eds R. H. Atrek, R. H. Gallagher, K. M. Ragsdell and 0. C. Zienkiewicz, Wiley, New York, 1984. 57. Saigal, S. and Aithal, R., A variational approach for the sensitivity of stress constraints using boundary elements. In Boundary Element Methods in Engineering, eds B. S. Annigeri and K. Tseng. USA, 1989, pp. 435-441. 58. Wang, S. Y., Sun, Y. and Gallagher, R. H., Sensitivity analysis in shape optimization of continuum structures. Computers and Structures, 1985, 20, 85s867. 59. Wech, M. and Steinke, P., An efficient technique in shape optimization. Journal of Structural Mechanics, 1984, 11, 433449. 60. Yang, R. J. and Botkin, M. E., Comparison between the

variational and implicit differentiation approaches to shape design sensitivities. AZAA Journal, 1986, 25, 1027-1032.

61. Yang, R. J. and Botkin, M. E., Accuracy of the domain material derivative approach to shape design sensitivity. AIAA Journal, 1986, 25, 16061610.

62. Yang, R. J. and Choi, K. K., Accuracy of finite element based shape design sensitivity analysis. Journal of Structural Mechanics,

1985, 131 xx-xi.

63. Zienkiewicz. 0. C. and Campbell, J. S., Shape optimization and sequential linear programming. In Optimum Structural Design, eds R. H. Gallagher and 0. C. Zienkiewicz. Wiley, New York, 1973.