Optimum design using viconopt, a buckling and strength constraint program for prismatic assemblies of anisotropic plates

Compurers & Swucrures Vol. 43. No. 4, PP. 699-708. 1992 Printed in Gnat Britain.

004s7949192 ss.00 + 0.00 0 1992 Pergamon Press Ltd

OPTIMUM DESIGN USING VICONOPT, A BUCKLING AND STRENGTH CONSTRAINT PROGRAM FOR PRISMATIC ASSEMBLIES OF ANISOTROPIC PLATES R.

BUTLERt

and F. W.

WILLIAMS

$3chool of Mechanical Engineering, University of Bath, Claverton Down, Bath BA2 7AY, U.K. IDivision of Structural Engineering, School of Engineering, University of Wales College of Cardiff, Cardiff CF2 IXH, U.K. (Received 9 April 1991)

Abstract-A computer program for obtaining the optimum (least mass) dimensions of the kind of prismatic assemblies of laminated, composite plates which occur in advanced aerospace construction is described. Rigorous buckling analysis (derived from exact member theory) and a tailored design procedure are used to produce designs which satisfy buckling and material strength constraints and configurational requirements. Analysis is two to three orders of magnitude quicker than FEM, keeps track of all the governing modes of failure and is efficiently adapted to give sensitivities and to maintain feasibility. Tailoring encourages convergence in fewer sizing cycles than competing programs and permits start designs which are a long way from feasible and/or optimum. Comparisons with its predecessor, PASO, show that the program is more likely to produce an optimum, will do so more quickly in some cases, and remains accurate for a wider range of problems.

learnt from the early simplified methods. PASCO included various practical design features and an approximate method, used to improve the VIPASA analysis of the overall modes of shear loaded panels. Later studies [8] showed that the results of such approximation were not always reliable. The later analysis program VICON [9] extended VIPASA to deal sufficiently accurately with all modes, but did not have the design capability of PASCO. Between the development of PASCO and the program which forms the subject of this paper, was the development of an approach [lo] which included ad hoc VICON modifications to PASCO’s VIPASA analysis. This approach, although inefficient, overcame the previous inaccuracy of PASCO and served to evaluate the structural efficiencies of various design concepts. VICONOPT (VICON with OPTimization) is a major step towards a design program which will ultimately supersede PASCO, having capabilities which either replace, include or are additional to those of PASCO. VICONOPT combines the analysis capability of VIPASA and VICON with a newly developed and efficient design technique to produce low-mass (i.e. hopefully near-optimum) designs satisfying buckling and material strength constraints and configurational requirements. Buckling analysis uses exact member theory and design uses the mathematical programming optimizer, CONMIN [l 11. VICONOPT is a 35,000-line FORTRAN 77 program available to U.S. users through COSMIC with details of release to other users available from the second author. Readers having access to the version of the program which was released in April 1990 [ 121

I. INTRODUCTION

Buckling criteria govern the design of almost half of the structure of most civil aircraft and more than half of most military aircraft. The buckling strength of composites compares favourably with isotropic materials and recent reductions in production costs have increased their availability to space, military and civil aircraft designer alike. Composite technology presents a manifold array of structural control parameters, making the problem of least mass design highly complicated and well suited to automated design by computers. Such automation has been an area of widening research over the last 20 years. Early studies resulted in computer codes for design of laminated plates [l] and specific plate asembly configurations [2,3] subject to simple loadings (usually axial compression). These studies, containing simplified buckling analysis techniques, served to give insight into the characteristics of efficient design. A more rigorous buckling analysis capability became available in the program VIPASA [4], which permitted a fairly general loading and configuration and considered the full range of local, intermediate and overall buckling responses. Various design methods using less rigorous analysis applied to less general configurations have been presented [5,6]. These may produce designs for which the analysis is inaccurate but are useful preliminary design tools and may present additional mass savings by allowing for some post-buckling. VIPASA buckling analysis was contained within the well established and widely used design program PASCO [7] which also adopted the design techniques 699

R.

700

BUTLER

and F. W.

should note that an earlier paper [13] describes the design features of that version. Some of the design features described in this paper, including strength constraints, have been added to the program prior to its second release so that this paper represents a considerably enhanced and extended form of the earlier paper. Analysis features have been presented elsewhere [14] and Sec. 2 of this paper describes the way in which these features have been adapted for use in design. Section 3 presents the means of combination of the analysis with the design process, whilst Sec. 4 gives example problems which first compare VICONOPT with its predecessor, PASCO, and then illustrate the wider range of problems to which VICONOPT can be applied. Typical computational requirements are also given. 2. ANALYSIS

METHOD

VICONOPT uses VIPASA and VICON analysis to converge on the eigenvalues of assemblies of anisotropic plates. These eigenvalues may be buckling load factors or natural frequencies of vibration, although current design use is restricted to buckling. Laminated plates are defined when the user specifies a series of layers with orthotropic material properties, orientated in any direction, and combines these to form the layup, or wall, of each plate. Plates are offset, rotated and connected to make up prismatic assemblies (e.g. see Fig. la) and are loaded by the longitudinally invariant NL, N, and N, of Fig l(b). The NLs can optionally be calculated by distributing a specified total axial load or an axial load per unit width, as might be obtained from a global finite element analysis of a complete vehicle. In addition, they may be calculated by distributing bending moments about the y-axis and z-axis or a nonlinear bending moment about the y-axis, due to a pressure or eccentricity over a portion of the cross-section, using beam-column theory. Multiple load cases can also be handled. In VIPASA analysis the plate differential equations are solved exactly to give stiffnesses corresponding to the U, u, w and Y displacement amplitudes, see Fig. l(b), of each longitudinal plate edge. Modes of buckling are assumed to vary sinusoidally with half wavelength 1 in the longitudinal direction and the theory includes coupling between in-plane and out-

(4 Fig. 1. Indication of the range of prismatic plates assemblies to which VICONOPT applies. (a) Representative cross-sections. (b) A component plate, showing the basic force system and the plate axis system, with .. .its associated displacement amplituaes.

WILLIAMS

of-plane systems and uses a fully populated [A] matrix. In addition there exists the option to make classical thin plate assumptions or to consider transverse shear deformation. An overall stiffness matrix [K] of order 4N is assembled for the complete structure using the stiffness matrix method and multi-level substructuring is used to reduce N, the number of nodes required to model the full structure. Longitudinal line supports may be prescribed to restrain any of u, u, w and Y. Transverse supports (parallel to the y-axis of Fig. lb) nearly always occur in practice (e.g. wing ribs and fuselage frames) but cannot be prescribed in VIPASA analysis. However, when all plates are either isotropic or orthotropic and carry no shear, the lines of zero displacement of the buckled plate assembly are straight and parallel to the y-axis. This is consistent with simply supported end conditions which are parallel to the y-axis provided the buckling half-wavelength i divides exactly into the length 1 of the plate. Otherwise, i.e. if any plates are anisotropic or more especially, loaded in shear, the solutions obtained approximate the results for such end conditions. Here the stiffness calculations involve complex arithmetic which produces skewed lines of zero displacement resulting from phase differences between displacements along transverse cross-sections. These results can, however, be expected to be realistic at short wavelength (e.g. i < 1/3) but will become excessively conservative as 1 approaches 1. PASCO used an approximate method to model the assembly as infinitely wide in the above circumstances rather than infinitely long, assuming better results would be obtained with boundary conditions satisfied transversely. Hence to satisfy the VIPASA requirement of invariant loading and stiffness in the infinite direction a ‘smeared stiffness’ solution was obtained with the assembly represented by a single orthotropic stiffness. A subsequent study[8] showed that the method could be up to 50% inaccurate. VICON analysis was developed to overcome the above difficulties by using Lagrangian multipliers to couple VIPASA stiffness matrices for an appropriate set of half-wavelengths, i, so that compatibility with an arbitrary supporting structure was achieved. Supports repeat at intervals of 1, hence a plate assembly of finite length f with simply supported ends may be modelled reasonably accurately by representing the simple supports by a line of rigid supports at x = 0. The results assume that the mode repeats over a length L = 21/t for some value 0 < 5 < 1. Each value of 5 generates an infinite series of I, although a small finite number chosen by selecting a suitable value of q in the following equation usually gives acceptable results 1

J.=<

+2m'

(m =O, +l, +2 ,...,

fq).

(1)

VICON and VIPASA use a special parabolic interpolation method [ 151 to solve the transcendental

Optimum design using VICONOPT eigenvalue problem represented by an equation of the type

Kl ID) = 101,

(2)

where {D} is the vector of nodal displacement amplitudes and the coefficient matrix [K,] is identical to [K] for VIPASA analysis, whereas for VICON analysis Ir<,] includes the coupled [K]s for each 1 given by eqn (1) and a matrix of constraints which couples them, while {D} includes Lagrangian multipliers as well as the displacement amplitudes for each 1. The analysis method [14] selects one of a list of alternative determinants derived from [K,] and, subject to certain conditions being met, uses parabolic interpolation through three known points on the determinant vs trial load factor curve to predict the eigenvalue as the zero of the determinant. An algorithm [16] is used to ensure that the points bound the required eigenvalue. The method recovers to a bisection step if more than o . eigenvalue is bounded, if the determinant plot contains a pole (i.e. vertical asymptote) between points or if the parabola through the three points is not sufficiently close to the straight line through their outermost pair. When several determinants satisfy these conditions the one with the least separation between the parabola and the straight line is chosen. The user is required to select values of 1 for VIPASA analysis and values of C and q for VICON analysis such that the full range of local, intermediate and overall modes is covered. VICON analysis is necessarily slower than VIPASA analysis involving a single value of 1, and use of the latter for short wavelength modes (e.g. 1 c f/3) and the former for longer wavelength modes is recommended. A recent study [14] showed that for the set of values of r = 0.0, 0.25, 0.5, 0.75 and 1.O, using q = 2 gave a sufficiently converged VICON result for the long wavelength mode of a shear loaded panel. The difference of about 10% between this result and that given by a finite element analysis arises because of the longitudinal continuity that VICON assumes (continuity which often occurs in practice along the length of aircraft wings and fuselages). The finite element analysis took [14] over 150 times longer than the critical VICON analysis and over 1100 times longer than any one of the non-critical, short wavelength VIPASA analyses that were performed to guard against local buckling.

701

design variables to become dependent variables. Linking maintains geometric consistency and allows time saving reductions in the number of design variables. Since VICON and VIPASA analyses both yield transcendental eigenvalue problems sensitivities cannot be obtained by techniques associated with the more usual linear eigenvalue problems. The partial derivative is therefore approximated to the following difference expression

aFyaxj

z (F; - Fi)/axj ,

for i = 1,. . . ,n,,

and

j = 1,. . . ,nDy,

(3)

where n,, is the number of buckling modes considered and nDy the number of design variables. Fi is the eigenvalue of the unperturbed (nominal) design and Fi the eigenvalue with thejth design variable x,, perturbed (incremented) by axj. VICONOPT first converges to a relative accuracy 100~ on Fi using the parabolic interpolation described above wherever possible. (c, the user-specified relative accuracy used by the program when sensitivities are not required should never be greater than 1 x 10-j and is typically 1 x 10m6.)For each perturbed design the revised plate stiffnesses are calculated before the n, perturbed eigenvalues are found using a quick approximate method. The approximate method assumes mathematical similarity between perturbed and unperturbed designs by displacing the parabola which gave Fi in order to predict Fi, see Fig. 2. The method requires only one iteration at load factor Fi to give the value of the determinant of [K:], i.e. IK:I where [K:] is the coefficient matrix of eqn (2) following perturbation. Appendix I shows how a suitable value of a is chosen and Appendix III illustrates the typical accuracy of the method compared with the full method of the following paragraph. The method recovers to the full method whenever any of the following checks are

Sensitivity analysis

The sensitivities of eigenvalues to small changes in design variables are used to direct the design modification sequence, preventing violation of critical and potentially critical buckling responses. Design variables may include plate breadths, layer thicknesses and layer ply angles. Other breadths, thicknesses and ply angles plus any plate (or substructure) rotations and offsets can be held fixed or can be linked to the

Fig. 2. Perturbed eigenvalue finding by the approximate method. The full line is the parabola for the unperturbed structure. Points 1, 2 and 3-[(F,,,I1y,(,), (Fii2,Iic12) and (F,, . IK..._,_ I, Jl are the three closest iterations u&d to find E. . .,_. the unperturbed eigenvalue and g, is the slope of & parabola as it passes through F,. )Kil, the perturbed

determinant at load factor F,, is found (point 4) and the perturbed eigenvalue F; is then determined by displacing the unperturbed parabola horizontally to make it pass through point 4.

R.

702

BUTLERand

activated: (1) convergence on F, did not use parabolic interpolation; (2) the maximum horizontal difference between the parabola and the straight line passing through the outer two points is greater than 1% of the horizontal distance betwen these two points, i.e. c* > O.O1(Fi,,-F,,,), see Fig. 2; (3) point 4 of Fig. 2 reveals a difference in the number of poles that are exceeded relative to points 1, 2 and 3; and (4), the bound on F, having the same sign determinant as point 4, e.g. point 2 of Fig. 2, has a different number of eigenvalues exceeded (calculated by the algorithm of [ 161)compared with the number exceeded at point 4. Activation of checks (1) and (2) for any mode i forces recovery for all corresponding j = 1, . . , n,, sensitivities whereas checks (3) and (4) force recovery for only the current j. The full method finds F: under all circumstances with certainty by analysing the perturbed problem fully. Convergence is enhanced by supplying close upper and lower bounds, FL, and FL, in the form

F. W.

WILLIAMS

t

nitial Stabilization

I

Sensitivity Analysi

Ft, = F,( 1 + 3~) FL = F,(l - 3a),

(4)

where 3ctF, is assumed to be an upper limit of eigenvalue change caused by a relative design variable perturbation 2. If either of these values does not bound F: a recovery procedure moves it progressively until it does. The choice of a in this case is discussed in Appendix II. 3. DESIGN METHOD

The design method is presented as a sizing strategy (see Fig. 3) consisting of a sequence of steps, each of which is described below. The initial analysis step performs VIPASA buckling analyses followed by VICON buckling analyses at the user supplied values of 1 and l, respectively. Analysis is accelerated with a ‘fast’ option, which omits unnecessary calculation of any eigenvalue which exceeds the lowest one found already. The initial stabilization step which follows factors design variable layer thicknesses to converge on a feasible design. Hence the thickness factor FSa which ensures that all material strength constraints and buckling constraints are ‘just’ satisfied is found and applied to the initial design. For reasons of efficiency F,, the factor for material strength feasibility is found first and applied as a lower bound during buckling stabilization. The latter is based on an earlier method [17] which uses the algorithm [16] to indicate the number of eigenvalues exceeded at any one iteration, corresponding to a trial value of Fsa. Convergence is by bisection and to save iterations, and hence reduce computer time, convergence accuracy has a low default value (1 x 10m3), which the user can adjust. The VIPASA or VICON analysis which gave the lowest (critical) eigenvalue during the initial analysis step is considered first, and this and all remaining

Fig. 3. The VICONOPT sizing strategy.

analyses are checked for stability so that F,, is only updated when an iteration shows that the design load exceeds the buckling load. Increase in plate thickness will always reduce layer stress or strain but an unusual form of design variable linking may cause layer stress or strain to increase. This situation is very unlikely and should it occur design is terminated. Similarly, plate thickness increases always cause monotonic increases in buckling load but changes in the linked dependent variables may cause buckling loads to decrease for some modes. This is equally unlikely but design proceeds without both the stabilization steps and the CONMIN cycling described below should it be detected. Stabilization keeps within any specified design variable bounds and design is terminated if bounds prevent stabilization. The constraint and sensitivity analysis step calculates a set of critical constraints and their sensitivities. This set is drawn from the buckling and material strength requirements specified by the user and is updated whenever a new sizing cycle is entered. The buckling constraints comprise the analyses, i.e. the user-supplied values of I for VIPASA analysis and < for VICON analysis, for which a single iteration at a limiting load factor FB showed that eigenvalues were exceeded. FB has default value 1.15, i.e. a load 15% greater than design load, for analyses involving long wavelengths and 1.05 for shorter wavelength analyses. If material strength is defined using maximum allowable value Sallowof stress or strain, or of some stress/strain function, a full set of material strength constraints is calculated for layers within the

Optimum design using VICONOP’T assembly. Calculated values S of the stress, strain or stress/strain function are found using classical thin plate theory and only the most critical S over the whole plate assembly, per layer specified in data, is stored in the set. If any constraint within this set is greater than or equal to C,, i.e. (S/S,,,,,) - 1 k Cy, where Cy has default value -0.25, it is added to the set of critical constraints. In the CONMIN optimization step CONMIN iteratively applies design variable changes using the method of feasible directions [l l] to minimize subject to

W((X})

G,,,({X})
m=l,...,n,

(5) (6)

and IX,) G PI G c&/l9

703

returns to the design at the start of the sizing cycle and automatically adjusts the limits placed on the CONMIN design variable change (i.e. the CONMIN move limits), re-using the constraints and sensitivities calculated in the time consuming constraint and sensitivity analysis step. This is done to ensure that eqn (8) remains a reasonable approximation to the buckling and material strength constraint values for the design returned from the CONMIN optimization step. The stabilized mass of the assembly will then approach the lowest achievable by the current round of the sizing cycle and design convergence will permit large CONMIN moves when the design is far from optimum and smaller moves as the optimum is approached. The move limits of the move limit calculation step for each sizing cycle k and CONMIN cycle n can be represented as follows:

(7) iXU} = IX) + R,@R)“{X)*

where (X} is the vector of design variables, W({X}) is the objective function, i.e. the plate assembly mass, n, the number of critical constraints and G,({X}) the critical constraint functions. {X,} and {X,} are upper and lower bounds on {X} and comprise the most limiting of the user supplied bounds and the CONMIN move limits described below. The values of critical constraint functions for the (q + 1)th CONMIN iteration are calculated using approximate used by earlier analyses similar to those methods [l, 71 and represented by the following Taylor series G,({XY+‘)

= G,({X]‘) + VG,({XD[{Xj 4+ ’ - (X}ql, (8)

in which VG,({X}) is the matrix of critical constraint sensitivities calculated in the constraint and sensitivity analysis step. The buckling and material strength constraint values produced by such an approximate analysis are incorrect following CONMIN optimization and the stabilization step of Fig. 3 returns infeasible or over-feasible designs to a ‘just’ feasible condition by factoring thicknesses in the way described above, i.e. as in the initial stabilization step. All user-supplied values of 1 and/or r are considered, starting with the critical 1 or < of the previous stabilization. If F,, becomes excessively large or small, indicating that a ‘just’ stable buckling condition cannot be achieved, the CONMIN optimization step is automatically re-entered with move limits increased four-fold. If stabilization is then unsuccessful move limits are reduced to a quarter of their original value. If both these recoveries fail design is terminated. The sequence of move limit calculation step, CONMIN optimization step, stabilization step and CONMIN cycle convergence test constitute one of the CONMIN cycles of Fig. 3. CONMIN cycling

(9) {XL] = {Xl-

&%),{X]*>

in which {Xc} and {X,} represent the upper and lower bounds on {X}, the vector of design variables at their sizing cycle start value. S is P, for the first sizing cycle and the smallest of (P,)k-’ and 36’/2 for the kth sizing cycle where R, and P, are the sizing cycle move limit factors. These have default values that can be adjusted by the user of 0.45 and 0.9, respectively, and 6’ is the value of S which gave the lowest value of stable mass during the previous sizing cycle. The {X}* is a vector of the largest of {X} and {X}i, the vector of design variables at their initial value. (R,), is the CONMIN cycle move limit factor for each CONMIN cycle n, and has value 1.OOfor the first CONMIN cycle of each new sizing cycle [(Rc ),I but is thereafter adjustable. (R,), is always halved for the second of the CONMIN cycles unless the mass of the plate assembly following CONMIN optimization, m,, is greater than that following stabilization, mos,in which case it is doubled. Figure 4 shows how the third CONMIN cycle searches for the value of (R,), that will give a minimum stabilized mass for the sizing cycle by using parabolic interpolation. This only occurs when: (1) the middle (R,), value of the three [(R,),, M,,] points required for interpolation [e.g. (R,), of Fig. 41 gives an moSvalue less than mk, the mass at the start of the sizing cycle; (2) the larger (R,), gives an mu3higher than that of the middle (R,), value, and; (3) the smaller (R,), gives an mos higher than that of the middle (R,), value. When (1) has not been achieved the next (R,), used is half of the smallest tried value of (R,),. When (2) has not been achieved the largest (R,), is doubled. The three points used in interpolation may include the stable mass at the start of the sizing cycle mkr for which (R,), = 0, because this coresponds to a CONMIN optimization step for

R. BUTLER and F. W. WILLIAMS

704

Fig. 4. Calculation of CONMIN cycle move limit factor (R,), using parabolic interpolation. The minimum of the parabola through three suitable points on the m, versus (R,), curve gives (R,), in this case, where m, is the plate assembly mass after the optimization step, m, is the mass after the stabilization step, rnkis the mass at the start of the kth sizing cycle and n is the CONMIN cycle number.

which {X) is not allowed to change. Four CONMIN cycles are the default maximum per sizing cycle but the user may specify more, in which case subsequent move limit calculation steps continue to work with the three best points given by the three lowest stable masses of previous CONMIN cycles satisfying criteria (l)-(3) above. The ‘stop CONMIN cycling’ of Fig. 3 occurs when a stable mass lower than mk has been found and either the CONMIN cycle number is equal to the maximum or the relative mass difference between two successive optimized configurations is within a limiting default value of 1% which can be adjusted by the user. If a stable mass lower than m, has not been found after the maximum number of sizing cycles has been reached further cycles are performed with (R,), + , = (R,),/2. Hence smaller and smaller moves are taken until either a lower stable mass is found, in which case the next sizing cycle is entered, or the move limit R,6(R,), of eqn (9) becomes smaller than a limiting value. In the latter case design is terminated since CONMIN cannot further reduce the stable mass and the final analysis described below is executed for the design at the start of the sizing cycle.

When CONMIN cycling has stopped a new sizing cycle is entered, see Fig. 3. The initial design of this sizing cycle is the design which has mass mk, where mk is the lowest m,, of the previous sizing cycle. New sizing cycles are entered until sizing cycling stops, i.e. until the user specified maximum number of sizing cycles is reached, or until convergence based on the relative mass difference of subsequent sizing cycles is within a limiting value, a user-specified minimum number of sizing cycles having been exceeded. The limiting value is user-adjustable with a default of 2%. The final analysis which follows calculates the constraint values of the final design. It is accelerated by use of the ‘fast’ option described previously, when describing the initial analysis step of Fig. 3, and optionally plots buckling mode shapes and/or gives sensitivities. 4. ILLUSTRATIVERESULTS The following five examples are presented to illustrate the capability of VICONOPT design compared with its predecessor, PASCO, including application to the shear loaded problems for which PASCO can be very inaccurate. CPU times, where given, show the efficiency of the program and refer to runs on a SUN4330. Strength constraints were only used for the variants of the first example described when presenting Fig. 6 and for the fifth example. Comparative PASCO and VICONOPT axially loaded compression panel

design for

PASCO design omits several of the steps given in Fig. 3. Thus it commences with the first constraint and sensitivity analysis step for the first sizing cycle, does not have a CONMIN cycle and proceeds as a series of sizing cycles each containing one constraint and sensitivity analysis step and one CONMIN optimization step. There is no stabilization and move limits are reduced in a prescribed way (which is not always successful) to encourage convergence. A set of critical constraints are retained for each sizing cycle mL VICONOPT ______.PASCO

NX (b)

Fig. 5. J-stiffened panel. (a) One of the eight identical cross-section portions showing design variable reference numbers j as superscripts. (b) Isometric view of complete panel with 0 defining positive ply angle. All dimensions are in mm.

Fig. 6. PASCO and VICONOPT design convergence for Examples 1 and 2.

Optimum design using VICONOPT Table 1. Definition of layups for plates of Fig. 5(a). I denotes layer thickness for initial design in mm, 0 layer ply angle and m material reference number. Material with m = 1 has: E, = 75.8424 GPa; Er = 5.51581 GPa; E,, = 2.06843 GPa; vu = 0.34 and; p = 1390 kg/m3. Material with m = 2 has: E, = 220.632 GPa; E2 = 18.6159 GPa; E,, = 6.41213 GPa; vu =0.21 and; p = 2000 kg/m3. j denote layer thickness design variable reference numbers Skin I

20a 20a 6a 4a

J-upper flange 0

mj

45 -45 0 90

1 1 2 1

4 4 5 6

t

8

m

j

4a 4a 35a 4a

45 -45 0 90

1 1 2 1

7 7 8 9

Symmetric J-web

J-lower flange

t

0

mj

4a

45 -45 0 90

1

4a

8a 4a

1 2 1

7 1: 9

t

9

m

j

4a

45 -45 0 90

I

11 11 12 13

4a 350 4a

1 2

I

Svmmetric Note: a = 0.0508. and an approximate

sensitivity

analysis

method

is

used to reduce analysis time. Example 1 is based on the panel of Fig. 5 simply supported along all four edges and loaded in longitudinal compression by a load per unit width, N, of 1152.98 kN/m (N,X,.= 0). The plates were all of symmetric carbon fibre composite construction with eight layers and Table 1 defines their layups. Substructuring was used to reduce the number of nodes in the final structure to two (one node at each longitudinal support). The design variables were the unthickened skin breadth, i.e. breadth between flanges, the stiffener lower flange breadth and the stiffener web breadth, plus all the layer thicknesses. Figure 5(a) and Table 1 indicate these design variable by giving their reference number i. The plate offsets shown in Fig. 5(a) were linked to the design variables which affected them and the unthickened and thickened skin breadths were linked to maintain a constant stiffener pitch of 152.40mm. Lower bounds of 10.00 and 0.0762 mm were placed on all design variable breadths and layer thicknesses respectively and upper bounds of 66.20 and 5.00mm were placed on the unthickened skin breadth and all layer thicknesses respectively. VICONOPT was made to use the VIPASA analysis inevitably used by PASCO. This was sufficiently accurate (anisotropy produced only slight incompatibility with transverse supports) and allowed direct comparison of the two designs since the problem comprised the balanced, symmetric laminates assumed in PASCO’s VIPASA analysis. The panel length I was 762 mm and analysis was applied for half-wavelengths 1 = I/m where m had values 1,2,... ,30. The L = I analysis was critical at the start of design with eigenvalue 0.865. The design was therefore not acceptable, but was nevertheless reasonably well chosen. PASCO initial move limits were set

705

at their recomended initial value of 20% with reduction of 20% at the start of each subsequent sizing cycle [i.e. in eqn (9), R, = 0.2, 6 is always (Ps)k-‘, where P, = 0.8 and (R,), is always 1.0 because there is no CONMIN cycling in PASCO]. All VICONOPT control parameters had their default settings. Figure 6 shows the panel mass produced by each sizing cycle k, divided by the mass mop, of the final VICONOPT design for design of Example 1 using PASCO and VICONOPT. The first five of the six sizing cycles required by VICONOPT used the maximum of four CONMIN cycles per sizing cycle and the sixth sizing cycle used only two. The mass and configuration of the final design is shown in Table 2. PASCO reached a minimum of different configuration and 3% higher mass after the recommended 10 sizing cycles. The VICONOPT run took just over 6 min of CPU, about the same time taken by PASCO. Three variants of Example 1 were run using VICONOPT and including strength constraints at different values of allowable material strain. Figure 7 shows the final design mass mf for each run compared to the final mass mop, obtained when, as above, strength constraints were not considered. Example 2 compares PASCO and VICONOPT design when a poorly chosen (highly unstable) initial design is used. It is the same problem as Example 1 except all layer thicknesses started at their lower bound of 0.0762 mm producing a critical initial eigenvalue of 0.00437, i.e. an initial buckling load of only about 0.4% of the design load, for 1 = l/6. The initial value of PASCO move limits was increased to 45% whilst equivalent VICONOPT parameters where kept at default settings. Figure 6 and Table 2 show that although the minimum achieved by VICONOPT was of different configuration, it was, as would be hoped, of similar mass to that of Example 1. PASCO, however, was unable to achieve a minimum. Its final design was unstable with critical eigenvalue 0.221 for L = l/4. The VICONOPT run in this case took Table 2. VICONOPT design of Examples l-5. Final values of design variables, x,, in mm and initial mass m,, mass at start of first sizing cycle m, and final mass m,, in kg Example number 3 4

j

1

2

1 2 3 4 5 6 7 8 9 10 11 12 13

66 14 36 0.57 0.08 0.25 0.26 2.23 0.13 0.92 0.20 1.01 0.17

64 15 38 0.13 0.08 0.99 0.37 1.28 0.23 1.19 0.28 0.61 0.32

mi mr m,

10.8 11.9 6.9

1.4 13.2 6.9

66 15 38 0.49 0.08 0.30 0.21 3.33 0.16 1.52 0.08 0.68 0.08 10.8 14.6 7.6

65 10 30 0.68 0.19 0.08 0.21 2.23 0.21 1.25 0.21 4.94 0.24 10.8 14.6 8.4

5 66 14 39 0.46 2.80 0.08 0.55 5.00 0.08 2.38 0.08 3.83 0.08 10.8 45.7 20.2

R. BUTLERand F. W. WILLIAMS

706

Allowable a.,, E2, 2~~ f5 Fig. 7. Mass of final VICONOPT design, m,, obtained with strength constraints at three different values of allowable strain compared with Q,,, the mass of the final VICONOPT design when strength constraints were not considered.

nearly 8 min of CPU, about 63% of the PASCO run time. YICONOPT

design for shear loaded panel

Example 3 was chosen to illustrate the use of VICON analysis in VICONOPT for the efficient and sufficiently accurate design of shear loaded panels. The problem is the same as Example 1 except that longitudinal compression and shear loads of N\ = 1576.14 and N,, = 78.807 kN/m were applied and the stiffener upper flange and skin directly above this flange were modelled as a single asymmetric plate to prevent the redistribution of shear which would otherwise be necessary. (It should be noted that the latter requires increased computational effort during stiffness calculations since the [B] matrix of this plate is non-zero.) Additionally, the final structure was modelled with eight nodes, one at each skin/stiffenerweb intersection, and point supports restraining w displacements were applied at each such node at x = 0 in the VICON analysis, to represent the transverse simple line supports. Longitudinal simple line supports also restrained w displacements, whereas Examples 1 and 2 had restrained both u and w displacements for comparison with the restraint of PAX0 simple supports. A combined VIPASA and VICON analysis was used to efficiently cover all buckling modes likely to govern design with sufficient accuracy. The analysis first considered seven values of VIPASA’s i, in the range l/30 to l/3 followed by five VICON responses corresponding to five evenly spaced values of 5 (1.0, 0.0, 0.75, 0.5 and 0.25 in that order). VICON constrained overall and intermediate buckling by coupling the largest two to five values of i in eqn (1) (i.e. q was set to 2) whilst VIPASA constrained local buckling. The critical eigenvalue in the initial analysis was 0.634 for c = 1.0. Design achieved 30% saving in mass after 47 CPU min involving six sizing cycles, and the critical eigenvalues

were then 1.00 for both VIPASA’s I = 85.7 mm and VICON’s 5 = 1.0. Table 2 shows the final values of design variables. A higher accuracy VICON analysis was performed for this final design by setting q = 5 in eqn (1). This gave the same critical eigenvalue of 1.00 for the l = 1.0 mode, indicating the sufficient accuracy of design using the combined analysis. Isometric and contour plots of the governing mode shapes for the final design are shown in Fig. 8. Example 4 illustrates the effects of using only VICON analysis to cover all modes that may govern design and indicates why this is not a wise course of action. It is the same as Example 3 except that it uses no VIPASA analysis and instead covers local buckling by using a value of q = 7 in eqn (1). Design achieved only a 22% saving time in mass after six sizing cycles taking 2 h and 46 min of CPU, i.e. over three times longer than Example 3, to achieve the minimum of Table 2. Here the mass savings of CONMIN optimization steps during the last three sizing cycles were subsequently lost through stabilization. This indicates cycling between modes which are not retained during the constraint and sensitivity analysis steps of these cycles but are violated following optimization. Such cycling also occurs in methods which use finite element analysis in design where it is similarly very difficult to identify a few representative critical modes from the large number of similar modes that can occur for typical panels. The combined analysis of the previous example presents no such problems because different modes are classified by separate values of J for VIPASA analysis whilst using low values of q for VICON analysis. The design of Examples l-4 did not consider constraints on material strength since no values of the principal strains c,, t2 and yj2 exceeded the chosen alowable values of + 0.004, kO.004 and +O.Ol, respectively, in any layer at any stage of design. Example 5 is the same as Example 3 but is more heavily loaded, causing material strains to govern design. It has an N,r of 7880.70 kN/m and an N,, of 394.035 kN/m and uses the same combined VICON and VIPASA analysis of Example 3. Design started with a critical strength constraint, (S/S,,,,, - l), of 2.61 and a critical eigenvalue of 0.127 for the VICON

Fig. 8. Isometric and contour VICONOPT plots of governing VICON < = 1.0 and VIPASA k = 85.7 mm modes, respectively, for final design of Example 3.

Optimum design using VICONOPT l = 1.0 mode. Initial stabilization required a thickness factor Fss of 4.70 to bring the governing strength constraint to 0.00. Design involving six cycles took 23 CPUmin and reduced the initial stable mass m, by 56% whilst satisfying all buckling and material strength constraints. Table 2 again shows the values of design variables for the final design in this case. 5. CONCLUDING

REMARKS

The way in which the VIPASA and VICON analyses of earlier programs have been adapted and combined in the optimum design program VICONOPT has been presented. Design permits the kind of in-plane loaded assemblies of laminated, composite plates that often occur in advanced aerospace construction and considers buckling and material strength constraints and configurational requirements. Computational demands are kept low and minimum mass designs satisfying constraints are found. It should however be stated that optimization problems of this type are always badly defined so that it is impossible to guarantee that any minimum represents a global optimum. For problems to which the earlier design program PASCO and VICONOPT both apply without approximation, VICONOPT is more likely to converge on a feasible design of minimum mass and will do so more quickly when the initial design is far from feasible. The features which account for these advantages are summarized by the following three paragraphs. Initial stabilization adjusts design variable thicknesses to achieve a ‘just’ feasible configuration, i.e. one that just carries the design load. Subsequent optimization is concerned with only reducing mass and unlike PASCO, convergence is successful when initial designs are poorly chosen. A critical set of buckling constraints are calculated from the eigenvalues of a critical set of VIPASA and VICON analyses. Each eigenvalue is found using a series of parabolic interpolations through the best three of the currently available load factor vs buckling determinant points (each point corresponds to a single iteration). The last such interpolation, which gives the eigenvalue, is reliably re-used in an approximate, sufficiently accurate sensitivity calculation method requiring one iteration per sensitivity. The method automatically recovers to a full method when data is unsuitable for approximation. In the first CONMIN cycle, optimization based on the critical set of constraints and their sensitivities is followed by stabilization which maintains feasibility, correcting the approximations of the optimizer. Optimization is then re-applied in a new CONMIN cycle using the same constraints and sensitivities with move limits tailored with the aim of giving a lower mass following stabilization than the previous CONMIN cycle. CONMIN cycling causes convergence c&i43/4-G

707

in typically six sizing cycles compared with the ten or more of PASCO. CONMIN cycling also adjusts initial values of move limits so that optimization brings together the eigenvalues of buckling analyses which are critical in the final design, thus permitting the time saving deletion of non-critical constraints. VICONOPT also extends the capability of PASCO by including VICON analysis. Constraints and sensitivities calculated for the overall or intermediate buckling of shear loaded, transversely supported plate assemblies using this analysis are sufficiently accurate, whereas the adjusted VIPASA analysis used by PASCO for these modes could be very inaccurate. An efficient and accurate method of combined analysis has been described which covers the full range of buckling responses for such panels. The method uses separate long wavelength VICON analyses (covering overall or intermediate buckling) and short wavelength VIPASA analyses (covering local buckling) to keep track of an adequate range of critical and near-critical modes of failure. Cycling from one set of critical buckling analyses to another which occurs when either VICON analysis is used alone or a general method of finite element analysis is used, is therefore prevented. VICONOPT analysis runs are typically two to three orders of magnitude quicker than finite element analysis and this efficiency is fully appreciated in the design process, permitting fairly interactive design use. Acknowledgements-The authors would like to thank Mr D. Kennedy of the University of Wales College of Cardiff (UWCC) and Dr M. S. Anderson of Old Dominion University for their many helpful suggestions. The work at LJWCC was supported by NASA, British Aerospace and the Science and Engineering Research Council.

REFERENCES

1. L. A. Schmidt

2.

3.

4.

5. 6. 7.

and B. Farshi, Optimum design of laminated fibre composite plates. Inr. J. Numer. Meth. Engng 11, 623-640 (1977). B. Argarwal and R. C. Davis, Minimum-weight designs for hat-stiffened composite panels under uniaxial compression. NASA TN D-7779 (1974). W. J. Stroud and N. Agranoff, Minimum mass design of filamentary composite panels under combined loads: design based on simplified buckling equations. NASA TN D-8257 (1976). W. H. Wittrick and F. W. Williams, Buckling and vibration of anisotropic or isotropic plate assemblies under combined loadings. Int. J. Mech. Sci. 16.209-239 (1974). J. N. Dickson and S. B. Biggers, Design and analysis of a stiffened composite fuselage panel. NASA CR- 159302 (1980). - D Bushnell, PANDAZ-program for minimum weight design of stiffened, composite, locally buckled panels. Comput. Srrucr. 25, 469-605 (1986). W. J. Stroud and M. S. Anderson, PASCO: structural panel analysis and sizing code, capability and analytical foundations. NASA TM-80181 (1981). [Supersedes NASA TM-80181 (1980).]

R.

708 8.

9.

10.

II.

12.

13.

14.

15.

16.

17.

BUTLER

W. J. Stroud, W. H. Greene and M. S. Anderson, Buckling loads of stiffened panels subjected to combined longitudinal compression and shear: results obtained with PASCO. EAL. and STAGS cornouter programs. NASA TP-22 I5 (1984). M. S. Anderson, F. W. Williams and C. J. Wright, Buckling and vibration of any prismatic assembly of shear and compression loaded anisotropic plates with an arbitrary supporting structure. Int. J. Mech. Sci. 25, 585-596 (1983). G. D. Swanson, Z. Gurdal and J. H. Starnes, Structural efficiency study of graphite epoxy aircraft rib structures. 29th AIAAlASMEIASCEIAHS Struct., Struct. Dynam. and Mat. Conf., Long Beach, California, AIAA Paner 88-2218. DD. 85-97 (19881. G.-N. Vanderphrats and F.*Moses, Structural optimization by methods of feasible directions. Compuf. Slruct. 3, 739-755 (1973). F. W. Williams, M. S. Anderson, D. Kennedy, R. Butler and G. Aston, User manual for VICONOPT. NASA CR-181966 (1990). R. Butler and F. W. Williams, Optimum design features of VICONOPT, an exact bucklina oroeram for orismatic assemblies of anisotropic plates: 31sr AiAA/ ASMEIASCEIAHSIASC Struct., Struct. Dynam. and Mat Cot$, Long Beach, California, AIAA Paper 901068, pp. 128991299 (1990). F. W. Williams, D. Kennedy and M. S. Anderson, Analysis features of VICONOPT, an exact buckling and vibration program for prismatic assemblies of anisotropic plates, 31st AlAA/ASME/ASCHE/AHS/ ASC Struct., Struct. Dynam. and Mat. Conf, Long Beach, California, AIAA Paper 90-0970, pp. 920-929 (1990). F. W. Williams and D. Kennedy, Reliable use of determinants to solve non-linear structural eigenvalue problems efficiently, Inr. J. Numer. Meth. Engng 26, 1825-1841 (1988). W. H. Wittrick and F. W. Williams, An algorithm for computing critical buckling loads of elastic structures. Q. J. Mech. Appl. Math. 24, 263-284 (1971). F. W. Williams and C. J. Wright, A design-orientated alternative to finding the initial buckling loads of prismatic plates assemblies. Aero. J. Royal. Aero Sot. 399-402 (1978). APPENDIX I: PERTURBATION APPROXIMATE SENSITIVITY

SIZE FOR METHOD

The relative design variable perturbation c( is chosen on the basis of a simple prediction of the likely change in eigenvalue (F; - F,), that will be caused by a design variable change ax,. The intention is to use a value of a which restricts the perturbed point 4 to the vicinity of the unperturbed points 1, 2 and 3 (see Fig. 2) whilst avoiding an ill-conditioned perturbed determinant (K:( A desirable condition would be that IrC:.( has similar value to the determinant of the third closest point to Fi. In the case of Fig. 2 the following would then be approximately correct -IK,l,lg,--F:-F,.

(Al)

Now no way has been foreseen in which a relative design variable perturbation CLcould cause a relative change of eigenvalue of more than 3a. Hence for the desirable approximate condition of eqn (Al) to be satisfied a value of u causing IR:I to be less than 11;(,1,would be u = -I~,l,/3g,~,.

WILLIAMS

and F. W.

642)

Table Al. Comparative sensitivities by approximate and full methods for VIPASA and VICON analyses of the shear loaded panel of Example 3. The error in the approximate method compared with the full method in this case is always less than 2% VIPASA I = l/6 Approx. Full

j

1

VICON < = 1.0 Approx. Full

2 3 4 5 6 7 8 9 IO 11 12 13

-37.7 125 14.0 841 1100 549 2190 653 739 1110 492 39.6 119

-37.8 125 14.0 839 1100 549 2190 652 739 1110 492 39.6 119

-1.76 25.2 30.4 129 211 77.9 665 92.5 176 340 44.0 172 13.2

- 1.79 25.3 30.5 127 210 77.7 663 91.9 176 339 44.0 172 13.1

a

0.000100

0.000333

0.000 154

0.000333

The maximum of this value and IOOOtis used for the a of the approximate method, where IOOOt is applied when the differences between F,,, and F,,r, and F,.r and F,, are of similar order. In this case the a of eqn (A2) may fail to perturb point 4 out of an ill-conditioned region close to F,.

APPENDIX

II: PERTURBATION SIZE SENSITIVITY METHOD

FOR FULL

The full method of sensitivity analysis finds F, and F: to an accuracy of + IOOeF,and + 100tF: respectively. If sensitivities are required to an accuracy of +p(aF,/ax,) say, so that the maximum absolute error for (F: -F,) does not exceed /3I(Fi - F,)I, the difference between F: and F, must be at least 2OOeF,/jI.This requires a relative perturbation of at least a = 2OOe/3fl(assuming again that 3aF, is an upper bound of eigenvalue change caused by a design variable perturbation of ax,). To be sure of a minimum difference between F, and F: it is necessary to use a higher value of a such that a = lOOyc/3fi

(A3)

and y is given a value based on experience (default of y = 5).

APPENDIX

111: COMPARATIVE ACCURACY SENSITIVITY METHODS

OF

The results of Table Al show the typical accuracy of sensitivities calculated by the approximate method compared to those calculated by the full method and the typical values of a used in each case. Sensitivities were found for the J-stiffened panel of Example 3 at the start of design (i.e. before the initial stabilization step) for the VICON 5 = 1.0 and VIPASA I = 116 analyses, for which the eigenvalues were 0.662 and 1.99, respectively. Values of 6 and j of 1 x lo-’ and 0.05 respectively were used and sensitivities calculated by the approximate method required 26 iterations, whereas those calculated by the full method required 122 iterations.