Optimal composite structures by deflection-variable programming

COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING 12 (1977) 1X5-179 0 NORTH-HOLLAND PUBLISHING COMPANY

OPTIMAL COMPOSITE STRUCTURES BY DEFLECTION-VARIABLE

PROGRAMMING

J. J. McKEOWN The Numerical Optimiwtion Centre, The Hatfield Polytechnic, 19 St. Albans Road, Hatfield, Hertfordshire, U.K. Manuscript received 10 November 1976 Revised manuscript received 7 February 1977 This paper deals with the problem of optimising multilaminar fibre-reinforced continua. A formulation is proposed in which primary optimisation variables are the nodal deflections of the finite element model. This leads to a decomposition of the problem into an inner and an outer subproblem; algorithms for solving both are described and numerical results given. It is shown that this approach is convenient for solving the nonlinear mixed integer program involved and allows the optimal number of layers in each element, the layer thicknesses and the fibre angles are to be determined simultaneously for very general inequality constraints on stresses and deflections. The case of structures of maximum stiffness is particularly considered.

Introduction This paper is concerned with the problem of designing minimum-cost structures whose elements are under states of plane stress. The elements of these structures are to be constructed from multiple layers of a linearly elastic orthotropic material; at any point on the structure each layer is distinguished from the others at that point by having a different angle (relative to some datum) for its axis of maximum strength and stiffness (its “fibre angle”). The number of layers may vary from point to point. The optimisation problem, therefore, is to determine at every point the number of layers, their fibre angles and their thicknesses such that some linear function of the thicknesses is minimised subject to constraints on stresses and deflections. To solve this problem numerically, it is convenient to idealise the structure as a set of N finite elements connected at nodes; the geometry of such elements and the form of the strain variation prescribed within them are arbitrary in the sense that the choice of these quantities is not a concern of the investigation to be described. With this idealisation the problem becomes that of determining the optimal number of layers within each element together with their fibre angles and thicknesses. In fig. 1, which ihustrates the layout of the ith element (assumed triangular), these quantities are referred to respectively as Lo,

ei,,ei,,... eLi

and

t;,’ . . . gj.’

The problem may be placed in perspective by considering the corresponding one for isotropic materials. There the variables L, and @f disappear, leaving only a single variable t’ for each finite element. Even with this simplification, however, the problem remains far from trivial. Indeed, it has been the focus of a considerable amount of effort, particularly in the field of aerospace structures, dating back to the early 1960s [ 11. The difficulty arises from the nonlinear relationship

156

J. J. McKeown,

Optinol composite structures b,v deflection-variable

programming

Fig. 1. The ith element.

between the design variable t’ and the deflections and stresses; the constraints involving these quantities are thus nonlinear and usually nonconvex. Gellatly and Berke [2] chronicle the direction of research in this area, which has resulted in some dissatisfaction with the mathematical programming approaches originally used. While these methods have a sound theoretical basis, and have given good results on problems which were not too large, they are expensive to use for the larger problems (several hundred design variables) typical of aerospace structural design. Because of this, less rigorous but more intuitive methods such as the “optimality criterion” approaches have been developed (e.g. [3]). It is not the intention of this paper to discuss the relative merits of the different approaches to the solution of the isotropic problem, but the very fact that mathematical programming methods have been discarded by some workers serves to illustrate the nontrivial nature of the problem. Returning now to the case of multilaminar structures, it is clear that the introduction of two more sets of variables must further complicate the problem as compared with those constructed from isotropic material. Additional complications arise from the fact that the Li are integer variables, and the stresses and deflections vary with the t3; in a very nonlinear way. Because of the integer programming aspect there seems to be no obvious way of extending the orthodox methods of mathematical programming, as applied to isotropic structures, to the multilaminar; and the discussion above would suggest that in any case no great success could be predicted. In fact, attempts to solve the problem have been mainly based either on ad hoc methods, using only a very small number of variables [ 4,5], or on optimality criterion techniques [ 61. In the latter case both the integer variables and the nonlinear character of the 0; have been avoided by restricting the allowable fibre angles to a small discrete set, typically O”, 90” and +45”. In every finite

J, J. McKeown, Optimal corrpositestructures by deflection-variableprogramming

157

element, layers are assumed to occur at every fibre angle in the set so that only the thicknesses of these layers (which may of course be zero) are the variables. The drawback of this formulation is of course the degree of approximation involved in thus restricting the set of allowable fibre angles. It is the aim of this paper to present a formulation of the multilaminar design problem which avoids the need for restrictions on the choice of fibre angle and which allows mathematical programming algorithms to be used. At the heart of this formulation is the proposition that the design variables Li, r3;and t; are not necessarily the best choice as optimisation variables. This distinction appears to have been seldom made by structural optimisers, perhaps because of the long tradition of structural design extending from the period before the availability of digital computers. Reinschmidt et al. [ 7 1 have in fact considered the effect of using simple transformations of the design variables in the isotropic case; in particular, they investigated the effect, in a simple problem, of using the “compliances” (corresponding to reciprocals of the t ‘> as optimisation variables to reduce the nonlinearity of the constraints. They were able to conclude that the problem which they considered could be made more tractable by a suitable choice of transformation, and, in particular, that it was beneficial to choose transformations which made the constraints more linear even at the expense of making the objective function nonlinear. The optimisation variables considered by Reinschmidt and his coworkers were simple functions of the design variables. In the present work a more radical departure from orthodox methods is proposed, namely that the optimisation variables should be chosen to be one of the sets of behaviour variables - in fact, the displacements at the nodes of the finite element mesh. There exists in general no (two-way) one-to-one correspondence between design variables and the displacements, so that this choice of variable imposes a quite different form on the problem; in effect, it is decomposed into two subproblems: One of these, the “inner” subproblem, involves determining the structure of minimum cost having a given set of nodal displacements under the given applied loads; the “outer” subproblem consists in finding the minimum of such fixed-displacement minima over the feasible set of deflections. Such an approach may initially seem a little abstract in that there is little intuitive connection between the optimisation variables and those quantities which can be directly controlled by a designer. However, it will be shown that it can lead to algorithms which solve the mixed-integer problem without difficulty, at the same time providing some useful insights into the nature of the optimal design. The account which follows is necessarily condensed and aims to do no more than to present an outline of the approach, followed by a brief description of some algorithms suitable for solving the inner and outer subproblems. The solution of two sample problems will be described.

1. A deflection-space formulation In this section the general problem described in the introduction will be restricted to exclude both direct constraints on stress and multiple load cases; we shall, however, retain constraints on the deflections. These restrictions will be relaxed in later sections. Let the following vectors of design variables be introduced: t = {ti,

t;, . . . tf Ll'

tf,

. ..

tQ

7

J. J. McKeown,

158

L=

Optimal composite structures by deflection-variable

programming

{L,,...L,} )

where N is the number of finite elements, which is of course fixed; the other quantities have been defined in the introduction. Any design D is completely defined once values have been assigned to t, 8 and L. The overall problem to be solved can now be expressed as follows: Li EyV=E

AiC i= 1

1 *

t;, j=l

subject to

6 = K-‘P, t; > 0

(1.1)

)

8; E Hi ,

L,integer

,

where Hi is the set of allowable angles for the layers in element [O, 71). The quantity 6 is the deflection vector defined as 6 =

(6,) 6,, ... S,}

i - usually

the half-open

interval

)

where M is the (fixed) number of nodal degrees of freedom of the structure. The deflections are due to a fixed load vector P, also of dimension M, and Ai is the (fixed) plane area of the ith finite element. The scalar functions fk(&) may be of any form. The deflections are related to the loads through a stiffness matrix K which in turn depends upon the design variables: 6=K-‘(t,o,L)P,

(1.2)

where K is linear in t but highly nonlinear in 0 so that its inverse is nonlinear in both of these sets of variables. Before considering the formulation of problem (1.1) in terms of deflections, it is convenient to prove the following optimality condition: LEMMA 1.Let D* be a design which solves problem (1 .I) and let 6* be its vector of deflections under load P. Then, of all allowable structures having this set of deflections under the load P, D* is that which has minimum volume. An allowable design is one which satisfies the direct constraints on the design variables in problem (1.1).

J. J. McKeown, Optimal composite structures hy deflection-variableprogramming

159

Proof: Since 6* corresponds to an optimal solution, it also satisfies the constraints fk(6*) < 0. Assume now that the lemma is false; then there exists at least one other allowable design, say D”, which also has deflection 6’ but is of lower volume. But this design would also satisfy all constraints (being allowable, it satisfies the direct constraints on t, 8 and L). It follows that D” is not optimal, which is a contradiction. Hence the lemma is true. Consider now the way in which problem (1.1) must be manipulated in order to express it in terms of 6 only. The deflection constraints are of course already in a suitable form, but the objective function poses a difficulty. In the form given, this function may be regarded as a means of associating uniquely a cost W with a design D(t, 8, L). Ideally, therefore, one would wish to associate a unique design D with every deflection 6 under the load P.Although (in the absence of instability) each D has a uniquely defined 6, the converse is hardly ever true. However, lemma 1 provides a means of associating a unique volume with every deflection. It is only necessary to use the fact that the optimum design D’ is uniquely related to its deflection 6’ in that, of all allowable designs having that deflection under P, the design D* is of minimum volume. If this necessary condition is enforced at every value of 6 under consideration, an objective function W(6)can be defined as follows: Let 6 be any M-vector; then W(6)is defined as the volume of that structure which has a deflection equal to 6 under the load P,which is allowable in the sense of lemma 1, and which is of lower volume than any other structure satisfying these conditions. If no such structure exists, then W(S)is of infinite value. Problem (1 .l) becomes: Min W(6) subject to k= 1,2,...4.

WI

It is clear that, in order to evaluate W(6), it is necessary to solve a subsidiary minimisation problem, and thereby to generate a design (or designs) D(S) which satisfies the necessary condition of lemma 1. Evidently this problem (the inner subproblem) can be expressed as follows: Li

ypLW’6’=k4J t;, ,, i=l

j=l

subject to K(t,e,L)S tpo,

=P, OE

W) Hi,

Li integer,

(where S is of course given). The equality constraints in (P2) are of course the equilibriumcompatibility conditions for the structure. Problem (1 .l> has now been decomposed into two subproblems. The outer problem is (Pl), and to solve it a search is required in the M-dimensional space of the deflections. During this search the subproblem (P2) must be solved each time the objective function W(8)is evaluated;

160

J. J. McKeown,

Optimal composite structures by deflection-variable

programming

(P2) is in terms of the original design variables, and is a mixed integer problem.

At this stage, therefore, it may not be obvious that the problem has been made more tractable by the reformulation. Clearly, (Pl ) does not directly involve integer variables, and it satisfies Reinschmidt’s general recommendation in that the constraints have been simplified at the expense of the objective function. In fact, in many practical problems the functions f’k(6) will be very simple indeed _ often only upper and lower bounds on S. However, W(6) is now an unknown quantity as far as its functional behaviour is concerned. In addition, the subproblem (P2), which will have to be solved repeatedly in the course of solving (Pl), is on the face of it a mixed integer problem not very much easier to solve than the original problem itself. Note that the equality constraints are nonlinear in 8. It will be the aim of the next section of the paper to show that in fact problem (P2) can be solved without difficulty by an algorithm which is a special generalisation of the simplex algorithm of linear programming. Furthermore, the function W(S) will be shown to posess useful properties from the point of view of solving (Pl), particularly with regard to the very important class of maximum stiffness structures.

2. The inner subproblem

(P2)

The inner subproblem will be considered first because the feasibility of the whole approach depends upon the ease of its solution. We begin by writing the constraints explicitly in terms of the design variables, bearing in mind that the displacement vector 6 is given. Let kipibe the stiffness matrix of the ith layer in the ith finite element, referred to some global Cartesian axes and scaled to unit thickness. Then it is shown in [9] that this matrix can be expressed as follows in the case of tibre-reinforced composites: ki*i = ka + k; cos 4ej + k\ sin 4%; + k\ cos 20; + ki sin 281 T

Cl)

where k’, (r = 0, 1, ... 4) are matrices which depend only on the fixed geometry of the ith element, and (3; is the angle between the fibre axis in the (i, i)th layer and the positive .u-axis. Using (2.1) we can write the total stiffness matrix of the structure as: N

K = c

Li c

{k;

+ k;

izl j=l

The equality

constraints

~0s

40; + k’, sin 40; + tii

cos

20;

+

kd sin 26;) tj .

(2.2)

in (P2) can be written: (2.3)

where p;(f?f) is the load vector (per unit of thickness) necessary to cause the (i, j)th layer to deflect by an amount 6 if it were isolated from the remainder of the structure. Using (2.1) it can be written explicitly as: Pj =pcl

i,i +p;i

cos

40; +p’;i

sin 40; + p’;j

sin 2e; +pij

cos 2e;

,

(2.4)

J. J. McKeown, Optimal composite structures by deflection-variableprogramming

161

where p:j = k;b (r = 0, . . . 4) can be evaluated immediately from the known k:j and 6. The inner subproblem can now be written in full as follows: =i

Min 5 A c ti, ’ j=l t,e,L i=l subject to

c c tfpi($) i=1 N

=i

j=l

t; > 0,

=P )

8; E Hi,

(2.5) L i integer.

Problem (2.5) has certain features in common with the classical linear programming problem; indeed, it can be viewed as such but with constraint coefficients which are nonconvex separable functions of variables 8’ whose optimal values (3;are required to be found, in addition to those of the linear variables t,!. In [ 81 the author has described a special algorithm for solving this problem, a detailed description of which is not possible within the limits of the present paper. Briefly, it proceeds as follows: It can easily be shown that a design which solves problem (2.5) cannot have more than M layers in total throughout the structure; the problem therefore is to determine how these layers are distributed and what their thicknesses and fibre angles should be. The algortihm, which will be referred to as the “functional linear programming” (FLP) algorithm, automatically generates an initial design havingM layers and satisfying the constraints in problem (2.5). Such a design is referred to as “basic feasible”, a term borrowed from classical linear programming practice. The procedure then continues iteratively, computing on each iteration the parameters of a new layer to be introduced in some element so as to decrease the total cost, and determining which layer already present should be discarded so as to maintain basic feasibility. It will be recognised that this logic resembles that of classical linear programming; indeed, the main difference between the two algorithms is that in the case of FLP the reduced gradients, whose values determine the choice of a new basic variable, are not scalars but functions of the 8’. Therefore, on every iteration a number of one-dimensional minimisations are required. This poses no difficulty because the functions to be minimised are simple enough not to require an iterative procedure for their solution. It will be clear from this synopsis of the algorithm that the optimal numbers of layers in each element (the Li) are generated automatically during the process of adding and discarding layers. Note that on any iteration every element is considered as a candidate to receive a new layer or to lose one which it already has; thus a layer may be added to element i and another dropped from element j, where i may or may not equal j. It follows no integer variables as appear in the problem. because the fibre angles appear as of the reduced functions, their optimal can be exactly subject to the fact that such values must in the sets Hi, can be either continuous, the algorithm generates infinite sequence which be terminated when convergence criterion has been In it contrasts with classical linear programming, which a finite Ref. 91 contains description of experience which been gained the

162

J. J. McKeown,

Optimal composite structures by deflection-variable

pro,q-amming

FLP algorithm, implemented as a FORTRAN program. Problems solved have been fairly small because of limitations imposed by computer size; typically, M and N would be between 30 and 45. No difficulties have been encountered in any of the cases. Perhaps the most striking result has been the enormous range of designs which can be found, all having identical deflections under a given load. The point is clearly demonstrated by the first of the examples given at the end of the paper. Here, the range was wide enough to include structures whose costs (in this case, volumes) differed from one another in the ratio 10: 1. Since the strains within the elements depend only upon nodal displacements, it follows that all such designs have the same strains. However, the stresses in any layer depend on both the strains and the fibre angle, so that the stresses usually differ widely throughout the range. It was stated at the beginning of this section that the feasibility of the whole deflection-variable approach depends mainly upon the ease with which (P2) can be solved. The FLP algorithm provides a means of overcoming the difficulties inherent in the problem, particularly its integer programming aspect, and is both reliable and reasonably efficient. Because this algorithm is in itself a considerable departure from other methods of constrained minimisation, it seems to offer considerable scope for further development both theoretically and numerically. For the moment, however, it can be claimed that the numerical experience so far gained (and described in part in [ 81 and [ 91 and in section 6 of this paper) provides a sufficiently firm base to justify one proceding to solve the outer problem (Pl ). Before leaving this section, a feature of the solution of the subproblem (P2), and hence of solutions to the overall problem, will be briefly noted. It has been mentioned above that the maximum number of layers in total throughout the structure is equal to M. This result is proved, e.g. [ 81 and [ 91; it follows very simply from the quasilinear programming form of problem (2.5). It is also proved in these references that the maximum number of layers possible in any substructure, such as an individual finite element, is equal to the number M' of independent modes of deformation of that substructure. For a three-node triangle the maximum number of layers is therefore three; for a six-node triangle, nine. For additional alternative load cases these limits are increased by M' - 1 for the first such additional load, M' - 2 for the second and so on. One of the unexpected bonuses of the deflection-variable formulation was this insight which it provided into the nature of the optimum structure.

3. The outer subproblem (Pl ) The solvability of the outer problem (Pl) depends both on the properties of the objective function W(6) and on the particular form of the constraintsf‘k(6). In practice, deflection constraints are seldom complicated, and in fact a very important class of structures will be shown to be characterised by only one linear equality constraint. The key question therefore involves the properties of the function W(6). It has been stated in section 2 that W(6) can be fairly easily evaluated using the FLP algorithm. Clearly, (P2) will not have a solution for every vector in M-space; for example the whole halfspace defined by the unequality P'6< 0 consists of a set of deflections for which any structures would have nonpositive strain energy; such a structure would violate the principle of conservation of energy. By definition the value of W(6) in that half-space is infinite. However, even in the positive-energy half-space there may exist regions in which no structure is defined; all such regions

will be described as “physically infeasible” to distinguish them from the merely infeasible regions for which the arbitrary constraintsfk(6) G 0 are not satisfied. Clearly, if (P2) can be solved for some deflection 6’, then a solution also exists for any deflection CY 6’, cr > 0; in fact the designs are the same for all positive values of 01,except that the thickness vector t and therefore the value of W is scaled by 1la. It follows that the boundaries of the physically infeasible regions are straight line generators passing through the origin in b-space. The existence of such boundaries is a potential source of difficulty with the method, and the behaviour of the function near them requires more investigation. To date, however, they have caused no computations dif~culties. It is shown in [ 91 that W(6)is continuous and differentiable everywhere except in and on the boundaries of this physically infeasible region. The differentiability of W(6)is a very useful feature from the point of view of solving (Pl). In fact the derivatives are easily computed in spite of the fact that W(S) is not of explicitly known analytic form. To see how this is done, consider the dual problem associated with (P2). This can be written [ 8,91: MaxP’X, 6,h

subject to pf(8’) A.< A, )

i=l,2

, ...

N,

where X is the M-vector of dual variables (Lagrange multipliers) straints in (2.5). Either by direct differentiation of (2.5), or by as the derivative of the objective function in (2.5) with respect of the constraint equations, we obtain the following expression

(3.1)

associated with the equality connoting that Xcan be interpreted to changes in the right-hand sides for the derivatives of W(6): (3.2)

where K is the stiffness matrix of the design which is the optimal solution to (3.1) for the given vector 6. Since the evaluation of W(6)involves the solution of (P2) and hence the determination of 1, it follows that eq. (3.2) allows us to compute the derivatives of W(6)with almost no additional effort. The usefulness of this property can be judged from the fact that to evaluate the derivatives approximately using central differences would require 2M + 1 solutions of (P2), where M, the number of nodal degrees of freedom, might well be large. The investi~tion of the properties of the function W(A) requires much more space than it can be given in a condensed account such as this; there is a more detailed description in [ 91. However, some observations will be made. Consider first the properties of 1, the vector of dual variables of (P2). This can clearly be regarded as a virtual deflection, being of the same dimension as B . Eq. (3.2) would then suggest that the vector -VW(&) be interpreted as an associated load vector. If this is done, and the virtual work theorem is applied to the two load-deflection sets (P,‘a) and (-VW(S), A), we obtain the following equation:

164

J. J. McKeown,

Optinal composite structures by deflection-variable

programming

P’1= --(VW@))’ 6. But it can be shown that the quantity ing partial differential equation: W(6) =

-6’VW(6) .

Ptl is equal to W(S), so that W(S) is seen to obey the follow-

(3.3)

Again, if k is interpreted as a virtual deflection then the dual (P2), (3. l), can be seen as the problem of maximising the total virtual energy under the requirement that the individual strain-energy densities p;(B’)l/Ai be less than unity. (Note that ~~(0’) is a unit thickness force vector, so that the strain-energy density interpretation is dimensionally correct). Taking this a little further, and recalling the definition of pi(@) embodied in equations (2.3), it can be seen that the typical constraint of dual (P2) takes on the following interesting bilinear form in terms of 6 and k :

Gfk,@)A/Ai < 1 . This form suggests the question: what is the significance of the special case where S is proportional to 1? In fact, this relationship is a property of maximum stiffness structures, which will be discussed in the next section. The only way to determine whether (Pl) is tractable or not is to choose some problem which involves a particular choice of the constraints fk(6> and to attempt to construct an algorithm which will solve it. In this way it should be possible to judge the effect of regions of physical infeasibility, an effect which at this stage of development cannot be predicted. To this end, it was decided to consider the class of maximum stiffness structures already mentioned - these will be dealt with in the next section. It must be emphasised, however, that the deflection space formulation is by no means limited to this class of problems; it is capable in principle of being applied to any structural problem with stress and/or deflection constraints of a very general form. This point will be further considered in a later section.

4. Maximum

stiffness structures

In this section a class of problems defined by a particular deflection constraint will be considered. A structure of maximum stiffness is defined to be one whose strain energy is a minimum for a given volume, that is one whose strain energy density is a minimum. It will now be shown that such a design is the same as that whose volume is minimal for a given strain energy. Let 6 1 be the deflection under a load P of a design solving the following problem: M$i W(S) , subject to

(4

p’6 = E, (constant) (D is the set of all allowable

designs).

J. J. McKeown, Optimal composite structures by deflection-variableprogramming

165

Assume that the associated structure D, is not of maximum stiffness, i.e. is not a solution to the following problem: Min P’6 D

subject to volume = W(S,) = constant.

09

Let E, be the strain energy of a design D, which does solve (B). We assume therefore that E, < E,. Let 6, be the deflection of D, under P. Thus P’ii, < P’6,

.

Define cr=E,IE,>l.

Scale D, by 1/a to produce a design whose deflection is u6,, whose strain energy i$(E,Ik,)E, , and whose volume is

=

E,

Hence, we have a design of lower volume than D, but of equal energy. But this violates the initial hypothesis that D, solves problem (A); therefore (A) and (B) are equivalent. It follows that the maximum-stiffness problem can be expressed as follows:

yin W(6) subject to P’6 = E,

(4.1)

where W(6) is of course the function defined in section 2, and E is a fixed energy constraint. Thus the problem has only one constraint, and that linear. The Lagrangian function associated with problem (4.1) is clearly L@, P) = W(6) - p(P'6 - E) ,

(4.2)

where p is a Lagrangian multiplier. For a solution to (4.1) we require that the following condition be satisfied: VL(S,p)=VW(6)+~=0.

Using (3.2) together with the equations relating P and 6, we obtain

(4.3)

166

J. J. McKeown,

Optimal composite structures by deflection-variable

programming

(4.4)

K(pu---A)=o. Here, the matrix K as well as the vectors 6 and L are those associated fying (4.1). Unless K is singular, eq. (4.4) implies /.lti=i.

with the structure(s)

satis-

(4.5a) (4.5) by P’ gives

Premultiplying

p = PQ/P’ii ,

(4.5b)

Eq. (4.5a) shows that if a solution to (4.1) exists, then for such a design the deflections 6 are proportional to the dual variables 1. This property of maximum stiffness structures has already been mentioned in the previous section. If a maximum stiffness design has values of the design variables denoted by O*, t* and L* then it can be shown that the inequalities in the dual problem (3.1) are satisfied as equalities for 13;corresponding to O*. Thus S’k’@j)l

=/li

i = 1, __.N .

)

(4.6

Using (4.5) and (4.6) we obtain --1 PfA iw(ef)S Ai Pf6

= 1.

Now P’6 = E is the total strain energy of the structure, timal structure. Thus, we can write: i

(stk’($)ij)

=

Ei/W, = E/W = constant

and P’l

= W(S) is the volume of the op-

,

I

where Wiis the volume of the ith layer. In other words, the strain energy per unit volume within all layers in every element is constant and equal to the strain energy per unit volume of the structure as a whole. This relationship is derived by Venkayya et al. [6] by another approach, and in fact forms the basis of their optimality criterion algorithm. It is a general property of maximumstiffness structures and is not limited to those which are multilaminar. But it will be shown below that the constant strain energy condition necessarily cannot always be satisfied, because of the existence of physically infeasible regions in &space. This problem is shown up very clearly by the deflection space formulation. Relationship (4.5a) suggests the following iterative procedure for solving (4.1): hk+l = ijk +sk(kk

- pkGk) )

(4.7)

where sk is a step-length whose determination will be discussed below. Clearly, when (4.5a) is satisfied, eq. (4.7) will imply no further change in 6. That iteration (4.7) implies downhill direc-

.I.J. McKeown, Optimal composite structures by dejlection-variableprogramming

167

tions on the function W(6) is easily shown as follows: Downhill directions will be defined as vectors d satisfying the condition d’ W(S) < 0 .

(4.8)

Using eq. (3.2) to eliminate the derivative and taking the direction defined by eq. (4.7), we then must prove the following inequality for the kth iteration (the super-fix k is dropped for simplicity):

and using (4.6) this would imply

Using P = Ki3 we obtain

i.e. we must prove that (6’KA)’ G (s’Ki5)(I.‘KA).

(4.9)

Since, by the well known property of stiffness matrices, K is at least positive semidefinite, it can be written in the form K=QtQ. Define the following variables:

Substituting these expressions in (4.9) we obtain

This is the Schwartz inequality, which holds for any vectors x,y in Euclidian vector space; thus the assertion is proved. The equality holds in either of two cases: (i)l ‘~6, i.e. at an optimum, (ii) Kl = o. In case (ii) the stiffness matrix is singular, that is it has at least one zero eigenvalue, which has an eigenvector coinciding with the vector of dual variables The simplicity of iteration (4.7), as well as its downhill property, recommends its use as the

168

.I.J. McKeown,

Optimal composite structures by deflection-variable

programming

basis of an algorithm for finding maximum-stiffness structures. The following is such as algorithm. Set k-= 0. For an arbitrary design, generate its deflections hk and its strain energyP’6k = E. where Kk is the stiffness Step 1 Compute W(sk), Ik, pk, dk = Ik - pkGk and VW(hk) = -Kklk, matrix of the optimal design having deflections ak. Step 2 If ldkl < e1 or ldktV Wk(sk) I < e2 (where e, and e2 are small positive numbers), then: Step 2a Set S* = hk, W* = W(lik), t"= tk,O*= f3k and L” = Lk, stop. Otherwise: Step 3 Using for example a parabolic unidirectional search technique as in [9], find sk in eq. (4.7) such that W(P+i) < W@). Step 4 If Sk < e3 (small positive number), go to step 2a. Otherwise set k = k+l and return to step 1. step 0

k=O

V

Generate

initial

feasible using

Sk and

physically

deflection

2

arbitrary

design

_6k+‘=

Gk - Sk 1

0

Do

Yes k = k+l Fig. 2. Algorithm

for finding

maximum-stiffness

structures.

J. J. McKeown, Optimal composite structures by deflection-variableprogramming

169

Fig. 3. Diagram for algorithm in fig. 2 (contours are those of W(6) projected on the constant energy hyperplane).

Figs. 2 and 3 illustrate this algorithm. For clarity the plane of diagram 3 is chosen to coincide with the hyperplane of constant energy, on which lie all points generated by the above algorithm. For this reason the coordinates axes in the diagram are not the deflections but a linear combination of them. It is important to realise that this algorithm can converge in a number of ways. Firstly, 16k+l - 8 k I can become very small in either of the two ways listed above, that is either at an optimum or at a point where K is singular and 1 is an eigenvector of zero eigenvalue. There is, however, another very significant way in which it can stop, namely when the search directions terminate on the boundary of a physically infeasible region. These regions constitute, in effect, additional constraints on problem (Pl) over and above the arbitrary constraints f’(s) < 0. They are unfortunately not expressible analytically (or do not seem to be), so that they cannot be easily included in the problem statement. However, their existence implies the possibility that the “solution” (the point where A = 1.16)does not lie in a physically feasible region: does not, in fact, exist. In such a case the algorithm can be expected to approach the boundary of the physically feasible region in which it has started and finally to halt when the criterion 1~~I < e3 is satisfied. This means that a point is generated such that either W(S) increases sharply after an initial decrease along the search direction or that a significant move along this direction takes the search into the physically infeasible region in which W(6) is defined to be infinite. There seems to be some evidence to suggest that the function may in fact tend to increase steeply as the boundary is approached from the physically feasible side, thus forming a natural barrier - but this point remains to be investigated. At any rate, such a “non-optimal” convergence may indeed be the best that any algorithm could be expected to do, and that it may be a perfectly satisfactory solution will be demonstrated by one of the numerical examples given in section 6. The comment has already been made that the only way to see whether any algorithm will work is to try it on problems. This is particularly true at the present stage in the development of the deflection-space formulation because of the imponderables involved - mainly the potential effect of physical infeasibility, For this reason the numerical results to be given in section 6 are particularly significant. However, before proceeding to discuss these results a section will be devoted to the removal of a simplification which was introduced in section 1, namely the exclusion from consideration of stress constraints.

170

J.J. McKeown,

Optimal composite structures bv deflection-variable

programming

5. Stress constraints In the statement of the overall problem given in the introduction, stress constraints were included; this would be the most typical case, with the constraint imposing some upper limit either on the individual stress components or on a quadratic function of them, as in the von Mises-Hill criterion. Buckling constraints, where the permitted stress levels depend on the layer thicknesses, will not be considered. In an isotropic structure the strains, and therefore the stresses, can be expressed as linear functions of 6 ; in such a structure it would. be unnecessary to include any stress constraints as such in the formal statement of the problem since they could be regarded simply as deflection constraints. In the FRP (= Fibre-Reinforced Plastic) case the situation is a little less simple. If the (i, i)th layer is subjected to a deflection 6, then the strains in that layer are everywhere defined by that deflection, no matter how many layers exist in element i. However, the stresses induced by such strains depend on the fibre angle in the layer. The stresses can therefore be regarded as depending in part on the fleflections and in part on the fibre angles. Put another way, if the deflections are held fixed as they are for the subproblem (P2), then the stresses depend on the f,ibre angles alone. The problem (2.5) can therefore be written in the following modified form when stress constraints are present: N

Min c

t,O,Li;l

Li

AiJg

t;,

subject to Li

N

c c t;p’(e;) i=]

=P

(5.1)

)

j=]

sp;, 6) given,

< 0

)

$2

i

0,

=

1, . . . N)

6; E Hi,

j = 1) . . . Li )

Li integer.

The stress constraints do not directly involve the thicknesses, and for this reason they do not affect the form of the optimisation algorithm. Their influence is confined to the choice of a new column to enter the basis on each iteration; where previously this involved a series of unconstrained one-dimensional minimisations, these now become constrained. There is no other difference, and, in particular, the upper bounds on the total number of layers and on the number of layers in each element remain unaffected. The modified subproblem (P2) expressed as (5.1) must be solved if the stress constraints are to be applied separately to each finite element. However, it must be remembered that the finite element model is usually itself an abstraction and that the element-by-element stress constraints are an interpretation, in terms of this model, of the stress constraint applied to the sheet as a whole: usually the stress limits would be the same in all elements. In some instances therefore it may be best to treat the overall maximum stress in the model as the quantity to be constrained and to scale the thicknesses of the whole sheet in order to satisfy this constraint rather than to control the stress separately in each element or layer. The appropriateness of this course of action will clearly depend on the deflection constraints applied to the outer problem; these must be such that the stresses in the optimal structure are reasonably evenly distributed. This can be expected

J.J. McKeown, Optimal composite structures by deflection-variableprogramming

171

to be the case with maximum-stiffness structures - this method of stress-scaling is used for example by Khot et al. [6]. It has also been adopted for the purposes of the numerical experiments to be described in the next section.

6. Numerical results In this section some examples of problems solved by the deflection-variable approach will be discussed. Before considering representative problems, however, it may be useful as a means of illustrating the general approach to consider a very simple structural design problem. In [ 81 the following was used to illustrate the ideas behind the FLP algorithm. Referring to fig. 4, a load (P, , P,,) is applied at a point A which is connected to given points along a straight rigid boundary by means of straight, pin-joined uniform elastic bars. Find the cross-sectional areas of the bars such that the loaded point deflects by a given amount (6,, 6,,) and the volume is minimal.

Pv , by

t

Fig. 4. A structure illustrating the FLP algorithm.

This will be recognised as a (P2) subproblem analogous to that for an FRP structure. If the number of support points is finite, then the problem is analogous to (2.5) when N = 1, H’ is a finite set, and M = 2. The vector ~‘(0’) is then the 2-vector of loads at A corresponding to a bar of unit cross-section making an angle 8’ with the horizontal (say) and having a displacement (6,, 6,,) at A. Since this problem is fully discussed in [ 81, both for discrete and continuous H’ , it will not be described further here. However, it can be extended to illustrate also the (Pl) subproblem; more specifically, we can easily investigate the nature of the function W(S) for this simple structure. Table 1. Positions of bars Bar x-coordinate (in)

1

2

3

4

5

6

7

8

-1.5

-1.3

-1.1

-1.0

-0.8

-0.6

-0.4

-0.2

9

10

11

12

13

14

15

16

0.2

0.4

0.6

0.8

1.0

1.1

1.3

1.5

Consider the particular case where P, = 1000/E, P,, = 0, and the permissible set of bars, sixteen in number, have the values for their base coordinates given by table 1. The set H is here discrete, and the corresponding (P2) subproblem is a straightforward linear programming one. Since M = 2,

172

J.J. McKeown,

Optimal composite structures by deflection-variable

programming

the maximum number of bars in any optimal fixed-deflection structure deflection problem was solved for a number of values of the deflection, in fig. 5. Each selected point is marked in the two-dimensional space of each point is the volume W(S) of the optimal fixed-deflection structure; range N,, N, of bars involved in the optimal arrangement. The diagram the load vector chosen.

is also two. The fixedand the results are set out the deflections - above below each point is the is of course specific to

Fig. 5. Results for the fixed-deflection structure (number above each point is W(6), volume of optimal structure; number-pair below each point is Nt , Nz, range of bars involved in optimal layout).

Fig. 5 illustrates the following general properties of the function W(6). (i) Along straight line generators through the origin, such as AB, AC, AD the optimal layout remains constant. (ii) Along such generators the volume decreases as the reciprocal of distance from the origin. (iii) There exists a limiting value of 6,/S,, namely 1.5 shown by the line AE, above which no structure exists. Since the diagram is clearly only the upper half of a symmetrical one, the physically feasible region is confined between two lines making an angle 2 tan-’ 1.5 with the positive b,-axis. The reason for this effect is simple: considering the upper half of the diagram, when 6,/6, exceeds 1.5, all bars in the given allowable set must go into tension. But to equilibrate the horizontal load, at least one must be in compression; hence the deflections become inconsistent with the load. The situation in the lower half is similar, with compression substituted for tension and vice versa. (iv) The value of W(6) increases as the physically infeasible region is approached. This property may be specific to the type of structure being considered.

J.J. McKeown, Optimal composite structures by deflection-variableprogramming

173

The maximum-stiffness structure is easily found in this case; we need only find the minimum value of W(6) subject to the condition P,6, = constant (since P,, is zero). Thus, since the value of the constant does not affect the layout, the maximum stiffness layout is found by seeking the minimum of W(S) along any vertical line in the physically feasible region. This is clearly the line 6, = 0 with a layout using bars 4 and 13 ; it is an orthogonal arrangement, and is in fact the Michell structure. Although the diagram does not record it, along this line the thicknesses of the bars are equal; all other layouts are unsymmetrical. The discussion above by no means exhausts the insights which can be gained from this simple example; a more extended discussion of it is included in [ 91. However, it illustrates the basic features of the deflection-space approach to structural optimisation. When multilaminar sheets are considered, the problem becomes more complex, both because of its increased dimension and because the problem of physical infeasibility is less easy to analyse. However, the main features of W(6) remain unchanged, while the maximum-stiffness problem becomes one of minimising W(S) not along a line but on a hyperplane. Two such problems will be discussed next: Problem 1. A cantilever problem (fig. 6), Problem 2. A cutout problem (fig. 7). Material properties are given in table 2.

Fig. 6. A cantilever problem.

Fig. 7. A cutout problem.

Table 2. Material properties Stiffness

Strength

.!I11= 30.0 X lo3 ksi

ok, (tensile)

E22 = 3.5 X lo3 ksi

q,

“12 = 0.210

0~~ (tensile)

v2r = 0.0189

qu

G

= 0.65 X 10’ ksi

= 198.0 ksi

(compressive) = 390.0 ksi = 11.4 ksi

(compressive) = 44.6 ksi

*LT,,

= 21.1 ksi

174

J.J. McKeown,

Optimal composite structures bJ1 dejlection-variable

programming

The first of these was chosen so as to approximate in outline to a Michell structure, and the finite element boundaries followed roughly the lines of principal stress in such a structure: these were merely devices to enable the layout of the final design to be compared with the Michell structure. The second problem is one solved in a restricted form by Venkayya et al. [ 61. The finite element idealisation used here is, however, very much coarser than that used by Khot et al. [6] because of restrictions caused by computer size. The results obtained on these and other examples by the deflection-space algorithms are described in detail in 191; a summary now follows: Problem 1. This example was solved starting from two distinct points in deflection space, obtained as follows: Starting point I. The deflection vector corresponding to a uniform plate one inch thick with each finite element having a single layer with fibre axis at zero angle. Starting point 2. As above but with fibre angles of -0.8 rad in elements above the centre line, and +0.8 rad below; also the thickness of the plate was not unity but was chosen so that both starting points had equal strain energies and thus ensured that final designs were directly comparable. The main features of the optimisation are summarised in table 3. Table 3. Data for problem

Starting Starting _

point point

1 2

1

Initial volume (in3)

Initial optimum volume (fixed deflection) (in3)

Final optimum Cin3)

37.6940 11.7062

3.5 395 4.4000

3.1357 3.1355

volume

16 - h/PI 161

Iterations

0.002 0.001

5 6

The second column of this table is included to show how much volume reduction was brought about simply by solving the subproblem (P2) once; it can be seen that this reduction could in fact be much greater than that due to the subsequent adjustment of the deflection. However, this comparison of volumes is made on the basis of equal strain energies; when volumes are compared on the basis of equal value of maximum stress, it has been found that the search in deflection space can result in significantly reduced values of this quantity. Table 3 illustrates very clearly the following points: (i) The algorithm converged to significantly the same values of optimum volume from both starting points. In fact, a study of the optimum designs found confirmed that the process had in both cases converged to virtually the same point in design space. fiil The dimensionless ratio I 6 ~ X/p l/l 6 I reached very small values. This quantity is a measure of the angle between the true and virtual deflection vectors; the values found here correspond to an angle of about 0.2”. Since a zero value implies a completely uniform distribution of strain energy density, it follows that the designs generated came very close to achieving this property. There can be little doubt, considering these points, that something very close to the true optimal design has been found in this case; it is illustrated in fig. 8. The result is satisfying for a number of reasons. Firstly, it indicates that an optimum solution can be found in a reasonable number of iterations, and that the (presumed) existence of physically infeasible regions need not cause the slightest difficulty. Secondly, it demonstrates that the algorithm converged to the same solution

175

J.J. McKeown, Optimal composite structures by deflection-variableprogramming

(a)

Fig. 8. Problem 1. a) Fibre angles (in each element, layers are numbered to match list of thicknesses in fig. 8b); b) Layer thicknesses (inches).

from quite distinct starting points - again without difficulty being caused by physically infeasible regions. Problem 2. The sheet-with-hole problem was run from one starting point, a deflection vector corresponding to a uniform plate of unit thickness with zero fibre angles in all elements - that is, all fibres aligned with the applied edge stresses. Table 4 summa&es the main results. Table 4. Data for problem 2 Initial volume (in3)

Initial optimum volume (fixed deflection) (in3)

Final optimum volume (in3)

IS - AlpI ISI

Iterations

17.227

14.267

12.683

0.57

11

These results are quite different from the previous one in two main ways. Firstly, the initial and final reductions in volume are not so spectacular; this is not surprising when it is considered that without the cutout these reductions would have been zero. More important is the fact that the solution is not optimal in the sense of condition 4.5a, i.e. I 6 - L//J I/ I61 % 0. In fact, the sequence of designs terminated close to the boundary of a physically feasible region. It has already been pointed out in section 4 that the algorithm could terminate in this way, and indeed that the point thus reached might well be the best obtainable. The design produced is shown in figs 9a and 9b. The thicknesses have been factored to give a maximum stress value of 1.O; the corresponding volume is 2.34 in3 for the quater-plate. This compares with the design of Venkayya et al. [61 which had a volume of 2.76 in3 and thus could be

176

J.J. McKeown,

o.

Optimal composite structures by deflection-variable

programming

lob

0.013

O.lSl

/

-.-a

\?SJ,

3

0.

I.00 .

w

1

w+

0.0003

0

/ 0.184

’

(b)

O.Olb

Fig. 9. Problem 2. a) Fibre angles (in each element, b) Layer thicknesses (inches).

layers are numbered

to match

list of thicknesses

in fig. 9b).

regarded as a satisfactory outcome; however, it must be stressed that the finite element mesh used in [ 61 was much finer than that used for the present work (138 elements compared with 33), so that an exact comparison cannot be made. These results are representative of those that have been solved using the deflection-variable approach up to the present time. Together with some other problems, they are described in much more detail in [ 91. The summary above does, it is believed, demonstrate the inherent feasibility of

J.J. McKeown, Optimal composite structures by deflection-variable

programming

177

the approach. No problem that has been attempted was unsuccessful so far.

7. Scope of the method Up to the present time, problems solved by the method of deflection-variables have necessarily been restricted in several ways. In the examples described in section 1, the most obvious restrictions have been as follows: 1. Only maximum-stiffness structures have been designed The decision to restrict initial tests to problems in this class was taken because of their simplicity when expressed in deflection-variable terms, combined with their importance as a class of practical problems. However, as has already been said, the method can be applied to a much wider class of structures. For example, if upper and lower limits are applied to the deflections, we have a problem with simple linear constraints and a nonlinear objective function W(S). This could be solved by a variety of available mathematical programming algorithms - for example, the method of projected gradients. The only special requirement would be that care must be taken not to enter physically infeasible regions during linear searches. If the constraints on the deflections are nonlinear, then a more powerful method would be required. Whatever the case, subproblem (P2) remains the same and deals with the integer part of the problem. 2. Multiple alternative loads have not been considered Again, this restriction has been imposed on grounds of simplicity; it is not inherent in the approach. Each additional load case increases the size of the deflection space since every load has a corresponding deflection. However, the space does not increase in dimension by the amount M because such deflections are, by the virtual work theorem, not independent. The alternative load-deflection pairs (pi, St) must obey the equations: Pi’Gi= SjP,

)

i = 1, 2, .. .

Q,

where of course Q is the number of alternative load cases. Thus, if 6, is chosen as a point in deflection subspace corresponding to P,, only M-l components of 6, can be arbitrarily chosen etc., so that the size of the composite space is M + (M-l) + .. . + (M-q+l) = MQ - Q(Q-1)/2. The maximum size of the space is therefore M(M+1)/2, when Q = M. Although the methods described in this paper are clearly still at an early stage of development, it is worth considering the contribution which they might, with very little extra effort, be able to make to the practical design of tibre-reinforced laminates. Their first advantage is the freedom which they allow in the choice of tibre angles. One can imagine the functional linear programming algorithm initially being used in conjunction with one of the optimality criterion methods in the following way. As has already been mentioned, such methods, of necessity, restrict the allowable fibre angles to a small set. However, if 6’ is the deflection of a design found by an optimalitycriterion method, the FLP algorithm could be used to find the minimum-volume design corresponding to this deflection without restriction on the fibre angles or with fewer restrictions. The result may not be truly optimal since 6* will not be an optimal deflection, but it is likely that it will be an improvement over the first attempt. In this way the FLP algorithm could be further tested and developed; at some stage algorithms for the solution of the outer subproblem could be introduced and tested in parallel with the established methods.

178

J.J. McKeown.

Optimal composite structures by deflectiorwariable

programming

The solution to the problems described in section 6 had one characteristic which would often be undesirable in practice; namely, finite elements which were “empty” in the sense of having no layers assigned to them. Normally some lower limit on the total thickness of material in any element would be a part of the problem statement. Such a constraint can be very easily applied as part of the functional linear programming algorithm; empty elements are not an inherent feature of the method. However, the fact that the algorithm can, if allowed, point to areas of the structure which perhaps ought to be empty is a useful feature because it suggests an interactive way of using the method for shape optimisation. Suppose, for example, that the cantilever structure of problem 2 had been defined not by a curve but by a rectangular envelope containing it. An initial application of the deflection-variable algorithms would be expected to produce empty elements around the upper and lower free corners. The FLP algorithm automatically prevents any node from becoming isolated, however, so that some vestigial elements would remain in these regions to allow all nodes to remain connected. A next step would be to redefine the finite element mesh to eliminate nodes which clearly served no purpose, thus allowing larger areas to be emptied on the next run. In this way the optimum outline, as well as the internal design, of the structure would emerge. This had been done in part for the finite element mesh of figs. 8a and 8b, where early experiments had shown that the root area between the two (final) points of support ought to be empty and was therefore not assigned an element. According to some authors, e.g. Willshire [ 1 11 the form in which fibre-reinforced composites will be initially used on a large scale in aeronautical structures will be as reinforcing on isotropic primary structures. Mixed problems of this type would present no difficulty for the methods described in this report; the problem (P2) would simply take on a form with some constraint columns pi fixed and some variable. The functional linear programming algorithm is capable of solving such problems without modification. Finally, any discussion of the scope of the deflection-variable method would not be complete without a reference to its generality; it is not just applicable to multilaminar structures. Problems involving trusses, space frames, rigid frames and isotropic sheets can also be approached in this way. The inner subproblem would differ in its form from one class of structures to another; in some cases it would be of classical linear programming form, and in others FLP. This would seem to be an interesting area for further research.

Acknowledgement The work described in this paper was done as part of an external Ph.D. project in association with the department of Aeronautics, Imperial College of Science & Technology, London. I would like to thank my Supervisor, Mr. Frank Matthews, for his help and advice. The numerical results were obtained using the DEC System 10 at Hatfield Polytechnic.

References [ 1 ] R.A. Gellatly, R.H. Gallagher and W.A. Luberacki, Development of a procedure for automated synthesis of minimum weight structures, AFFDL-TR-64-141 (1964). [2] R.A. Gellatly and L. Berke, Optimal structural design, AFFDL-TR-70-165 (1971). [3] V.B. Venkayya, N.S. Khot and V.S. Reddy, Energy distribution in an optimal structural design, AFFDL-TR-68-156 (1968).

J.J. McKeown, Optimal composite structures by deflection-variableprogramming

179

[4] R.N. Hadcock, Boron-epoxy aircraft structures, in: G. Lubin (ed.) Handbook of fibreglass and advanced plastics composites (Van Nostrand-Reinhold, 1969). [S] C.W. Rogers, Structural design with composites, in: R.T. and H.S. Schwartz (eds.) Fundamental aspects of fibre reinforced plastic composites (Interscience, 1968). [6] N.S. Khot, V.B. Venkayya, C.D. Johnson and V.A. Tischler, Application of optimality criterion to fibre-reinforced composites, Unnumbered AFFDL report (1975). [7] K.F. Reinschmidt, C.A. Cornell and J.F. Brotchie, Interative design and structural optimization, Proc. ASCE, J. Strut. Div. 92 (1966) 281-318. (81 J. J. McKeown, A quasi linear programming algorithm for optimising fibre-reinforced structures of fixed stiffness, Int. J. Comp. Meths. App. Mech. Eng. 6 (1975) 123-154. [9] J.J. McKeown, A deflection-variable technique for the optimisation of fibre-reinforced composite structures, Ph.D. Thesis, Dept. Aeronautics, Imperial College of Science and Technology. To be submitted 1977. [lo] A.C. Ham and A.J. Willshire, Advanced materials and their use in civil aircraft structures, Paper presented at Roy. Aero. Sot. Spring Convention (May 1976).