A finite element programme for the analysis of laminated fibre reinforced shells of revolution under axisymmetrically applied loading

UNLIMITED

C o l l e g e of A e r o n a u t i c s Memo 76O6 June 1976

TECHNISCHE HOGESCHOOL DELFT LUCHTVAART- EN RUIMTEVAARTTECHNIEK

BIBLIOTHEEK

Kluyverweg 1 - DELFT

7 DEC. 1877 A FINITE ELEMENT PROGRAMME FOR T H E

ANALYSIS OP LAMINATED PIBRE REINFORCED SHELLS OF REVOLUTION UNDER AXISYMMETRICALLY

APPLIED LOADING

R . C . PEWS R. TETLOW

CRANFIELD INSTITUTE O P TECHNOLOGY CRANPIELD

BEDFORD

Prepared for Procurement Executive Ministry o f Defence

Contract Number A T / 2 0 2 8 / 0 9 V R P E UNLIMITED

All recipients of this report are advised that it must not be copied, in part or in whole, or be given any further distribution outside the Ministry of Defence, without the written approval of the Assistant Director Research, Guided Weapons, Procurement Executive, Ministry of Defence.

College of Aeronautics Memo 7606 June 1976

CRANPIELD INSTITUTE OF TECHNOLOGY

A FINITE ELEMENT PROGRAMME FOR THE ANALYSIS OF LAMINATED FIBRE REINFORCED SHELLS OF REVOLUTION UNDER AXISYMMETRICALLY

APPLIED LOADING

by

R.C. PEWS MSc, CEng, MIMechE. R. TETLOW MSc, CEng, MRAeS.

SUMMARY

A two noded axisymmetric shell element is used as the basis of a finite element programme for the analysis of laminated anisotropic shells of revolution.

Elements form conical frustra, to accommodate tapered or domed shells and both uniform and abrupt variations in shell wall thickness from element to element are catered for.

Output from the programme includes nodal deflections, element stress resultants and stress couples and if desired, layer by layer stress evaluation through the shell thickness, together with a strength assessment.

Work carried out as part of M.O.D, (P.E.) Agreement No .AT/2028/09 V R P E .

CONTENTS ?a^ NOTATION 1. INTRODUCTION 1 1.1 Previous Work 1 1.2 Current Task 1 1.3 Initial Approach 2 1.4 Choice of Element 2 2. STRUCTURAL DEFINITION 2 2.1 Co-ordinate System 2 2.2 Element Construction . 3 2.3 Loading Environment 3 2.4 'step' Discontinuities 3 3. STRESS ANALYSIS 4 3.1 Nodal Displacements 4 3.2 Stress Resultants and Stress Couples '4

3.3 Interlaminar Shear Force 5

3.4 Sign Convention 6 3.5 Layer Stresses and Failure Criteria 6

4. PROGRAMME INPUT AND OUTPUT DATA 7 4.1 Programme Input Information 7

4.2 Programme Output 13 4.2.1 Title and Salient Constants 13

4.2.2 Materials Data . 1 3

4.2.3 Geometry Data 14 4.2.4 Element Data 14 4.2.5 Constraint Data 14 4.2.6 Input Load Vector 14 4.2.7 Nodal Displacements 14 4.2.8 Stress Resultants and Stress Couples 15

4.2.9 Layer Stresses 15 4.3 Data Suppression 15 5. DISCUSSION • . 15

5.1 Limitations 15 5.2 Supplementary Storage Requirements 16

5.3 Possible Future Developments l6

APPENDICES

A ELEMENT STIFFNESS MATRIX B LAMINATED PLATE THEORY

C INPUT DATA GENERATION

D WORKED EXAMPLES AND PROGRAMME LIST

FIGUIIES

2. Axisymmetric Shell Element 2. Layer Definition

3. Transverse Shear Force Evaluation 4. Generalised Hookes Law

5. Stress Resultants and Stress Couples 6. Interlaminar Shear Stress

7. Geodesic Isotensoid End Closure 8. Data Generation Example

9. Element Grid for Pressure Vessel 10. Element Grid for Test Cylinder 11. Isotropic Cylindei' Results 12. Pressure Vessel Results

NOTATION

B Element Strain - displacement matrix D Element elasticity matrix

E Youngs Modulus

F Layer strength, or reserve factor G Shear modulus

K Stiffness matrix L • Length of an element M Stress couple

N Stress resultant

Q Transverse shear force R Internal radius

T. Total shell wall thickness X Longitudinal tensile strength Y Transverse tensile strength Z Shear strength

h Distance of a layer from the reference surface r Hoop radius

Radial co-ordinate s Local co-ordinate

Nodal displacements

P Rotational displacement

5 Column matrix of element nodal displacements e Direct strain

Ö Fibre orientation relative to a meridian X Change of curvature

a Direct stress T Shear stress

^ Cone angle of element

SUBSCRIPTS

I Parameter referred to the direction of fibres

t Parameter referred to a direction at right angle to fibres

It Parameter referred to the l,t plane s Parameter referred to the (s) direction

Ö Parameter referred to a circumferential direction k Parameter referred to the kth layer

1

-1. INTRODUCTION 1.1 Previous Work

Prior to the establishment of the current contract, a programme of work was completed under M.O.D.(P.E.)

Agreement No.AT/2028/059GC, relating to fundamental methods of design with fibre reinforced materials. In the initial stages of this work special emphasis was given to the

application of these 'modern materials' to anti tank guided weapons and a detailed study of a launcher tube component was carried out. (Ref.l). In the execution of this work

laminated plate theory was established and computer

programmes written on cylindrical buckling phenomena (Ref.2), strength evaluation (Ref.3) and the evaluation of elastic

constants for multi layered plates (Ref.3).

In the later stages 'of this work programme' effort was concentrated on the problems of pressure vessel design with fibre reinforced materials and computer programmes were developed for the evaluation of stress resultants and stress couples in common areas of structural discontinuity. (Refs 4 and 5 ) . Attention was also given to the calculation of individual layer stress for these structures and a 'yield type' failure criteria was used to assess strength reserves. (Ref.6).

In conclusion the methods developed in Refs ( 1 - 6 ) were used to evaluate design charts for a number of the more commonly used materials and a parametric study was completed on a typical rocket motor case component. (Ref.7)•

1.2 Current Task

There is a growing interest at M.O.D.(P.E.) in fibre reinforced materials and it has recently been suggested that a series of fibre composite rocket motor cases be produced to a current missile specification, for materials and

manufacturing verification.

Computerised initial design procedures facilitating the variation of fibre orientations to give the optimum

compromise between strength and stiffness requirements,

already exist at Cranfield and can be applied directly to the components in question. The current contract AT/2028/094/RPE was therefore initiated by R.P.E. Westcott to make use of the existing Cranfield work in the verification of their

components and also to extend the programme to encompass the finite element method of stress analysis.

2

-1.3 Initial Approach

In recent years a great deal of research has been directed towards the finite element method and many suites of programmes such as NASTRAN, PAFEC, etc, have been developed. Such packages necessitate extreme generality, so that they

can be applied to a vast number of engineering problems. For specific applications this in built generality often makes the use of these techniques rather cumbersome and it was therefore decided to write a small programme to fit the specific

requirement of R.P.E. Westcott. 1.4 Choice of Element

Most structures of interest to R.P.E. Westcott are axisymmetric and in the interests of economy therefore, the obvious choice is an element that can take account of these axisymmetric properties. Already operational at Westcott is a suite of programmes acquired from the Rohm and Haas Company of Alabama U.S.A. (Ref.13) and it was agreed that the

Cranfield programme should if possible conform with this system. The Rohm and Haas programme uses an axisymmetric triangular element with two degrees of freedom per node and this element must therefore be a contender for development in the current work. This element is however most useful for thick walled structures and though it can be used for shells is considered to be inferior to an element specifically

tailored to shell requirements. The terms of this contract are hinged around a filament wound rocket motor case, which

is primarily a shell type structure and it was therefore

decided to concentrate Initial effort on an ajcisymmetric shell element. The type of element chosen has three degrees of

freedom per node (two displacements and a rotation see Figure 1) and forms a conical frustrum,thereby facilitating the analysis of tapered, or domed shells of revolution, as well as parallel cylinders. Use of the element for isotropic materials is

demonstrated in Ref.8 and a full description of the element stiffness matrix derivation, together with the elasticity matrix modification to accommodate directional materials is given in Appendix A.

2. STRUCTURAL DEFINITION 2.1 Co-ordinate System

The axisymmetric structure is defined in terms of cylindrical co-ordinates (z,r), where z,r indicates a circle of radius (r) located (z) units from the origin, measured along the axis of rotation (see Figure l ) . Individual elements form conical frustra and are therefore two noded, requiring the definition of two co-ordinate circles (z,r)j^ and (z,ryi. A typical element is shown in Figure 1 and it can be seen that the structure as a whole, can be defined by a string of such elements linked end to end.

3

-2.2 Element Construction

The thickness (t) of an element is assumed to comprise a number of layers of directional material. The location of each layer in the overall thickness is described relative to an arbitrary reference surface, as indicated in Figure 2.

The programme treats layer data in 'block' form and splits each input layer into six layers of equal thickness, so that the cross section shown in Figure 2 consisting of four blocks, is idealised to a total of twenty four individual layers. Two fundamental block types are possible and for descriptive purposes will be designated Blocks type A and B. Block type A (blocks 1 and 4 Figure 2) is a cross plied block and is

assumed by the programme to be constructed from uni-dlrectional layers, with filaments orientated alternately (from layer to layer) at ±0° to the direction of the axis of rotation (axial direction). Block types B (blocks 2 and 3 Figure 2) have the fibres of all six of their constituent uni-dlrectional layers, orientated at +6° to the axial direction.

Materials used throughout the cross section are defined by four elastic constants, E,, E^, \x-,^ and G,^ and need not

therefore be directional. If required a total of five different materials can be used in the construction and five blocks

similar to blocks I-4 Figure 4 can be accommodated through the material thickness. It is therefore possible to describe a total of thirty layers at thirty locations through the material thickness and stressing information can be output at each of these stations.

2.3 Loading Environment

The displacement field for the element chosen refers specifically to axisymmetric deformations and therefore all structural loading must be axisymmetric. The most obvious

form of axisymmetric loading is pressure loading and a facility exists within the programme to automatically generate the

pressure loading vector, through linear condensation of the pressure forces to the nodal circles.

Axisymmetric concentrated loads are also catered for and when required, are added automatically into their correct locations in the input loading vector. When concentrated loads are included in the loading data it is important to realise that it is the whole circumferential value of the input load that the programme requires.

2.4 'Step' Discontinuities

Where the total shell thickness is varying from element to element, unless specifically requested, it is assumed during element assembly that adjacent elements are linked via their neutral axes. Some structures requiring analysis may however contain step discontinuities in wall thickness and in such

cases a bending moment appropriate to the shift in neutral axis location is added to the overall structural load vector.

- 2^

-3. STRESS ANALYSIS

3.1 Nodal Displacement

The general equation on which the displacement method

of analysis is based Is as

followsj-{_F}

= [K] {5} ...(1)

Where {^P"} is a vector representing applied loading

[ K ] is the structural stiffness matrix and {_ 6 ^ a vector of

nodal displacements. The unknowns in equation (1; are the

nodal displacements and their derivation requires the solution

of 3N simultaneous equations, where N is the total number of

nodes used to describe the structure,

For a continuous shell of revolution defined by a

series of elements of the type discussed in Appendix A the

stiffness matrix [K] is highly banded (max semi - band width

of six) and the programme solves for the nodal displacements

using a banded solution technique as suggested by Zienkiewicz

(Ref,8).

3.2 Stress Resultants and Stress Couples

With the nodal displacements in terms of global

co-ordinates available, element stress resultants and stress

couples can be evaluated from the following

^^s

= [D] [B] [T]03

.^e-^ ...(2)

Where [D] is the element elasticity matrix, [ B]

the element strain - displacement matrix, L T 4 a transformation

matrix and \_&y the 6 x 1 oolumn matrix of nodal displacements

for the element, (see Appendix A for more detailed definition

of these matrices.)

The matrix C B ] in equation (2) is dependent on the

element local co-ordinates and therefore the stress resultants

and stress couples can be calculated at any point in the element,

During preliminary runs of the programme it was observed that

best results were achieved when co-ordinates appropriate to

the element centre point were used in [ B ] and therefore all

stressing information is referred to the centre of respective ^

5

-3.3 Interlaminar Shear Force

The axisymmetric interlaminar shear force is

numerically equal to the rate óf change of meridional bending moment, ie,

dMg

% =

di-The meridional bending moment Mg is available as an explicit function in (s) from equation (2), but this involves the element shape function, which can be discontinuous in terms of curvature (d^o/ds^) at the element boundaries. It was therefore considered to be unreliable to evaluate Qg from equation (2) directly and in the programme the following

procedure is used,

The centre element values of Mg as derived from equation (2) are stored for each element together with corresponding peripheral co-ordinates (s). Consider the

co-ordinate points referring to elements i to 1 + 3 (Figure 3)

e point i i s

" i + 1 "

" 1 + 2 "

1 + 3 "

defined

11 f t t i by M ( l ) , " M(2), " M(3), " M(4), s ( l ) s ( 2 )

s ( 3 )

s(4)

where M(l) = (Mg)^ - (Mg)^ ^^(2) = (Ms)i + 1 - (Mg)^ etc s(l) = s(i) - s(i) s(2) = s(i + 1) - s(i) etc

The co-ordinates (M,s) are now used to identify the three coefficients of the least squares curve fitting

parabola passing through the four data points, 2 ie M = a + bs + cs *

•••

^ s = f

= ^ ^ 2CS

at element 1 s = 0 .*. (Qg) = b

6

-For (Qs)j^+i elements i+1 to i+4 are used and the cycle repeated until the element (N-3) Is reached, where (N) is the total number of elements used in the structure,

then M(l) = (Mg)j^,^ - (^1^)^^.^ f^(2) = (Mg)N-2 - ( M S ) N - 3 etc. s(l) = (S)N_3 - s(N-3) s(2) = (s)jj_2 - S ( N - 3 ) etc. and (^S)N-3 ^ ^ ( Q S ) N - 2 = ^ + 2cs(2) (^S)N-I = ^ ^ 2cs(3) ( Q S ) N = b + 2cs(4) 3.4 Sign Convention

Stress resultants N and N Q are positive when tensile and stress couples M_ and M^ are positive when they imply

s o

tensile stresses on the inside surface of the shell. 3.5 Layer Stresses and Failure Criteria

As has been mentioned previously, for an axisymmetric deformation there can be no circumferential displacement

component and therefore no residual 'in*^lane' shear stress over the shell thickness as a whole. For laminated fibre reinforced shells however, in-plane shear forces (N g) will exist from layer to layer, but will sum to zero through the total shell thickness. Resulting in-plane or interlaminar shear stress can significantly effect the overall structural strength for certain ply orientations and must be included in any failure criterion applied to the layers,

It can be shown that individual layer stresses are given by the following matrix equation

"e

^ L e J ^ j N j .([f]''-z''[c]'^)Jo

7

-where the 3 x 3 matrices [ e]

defined i n Appendix B,

(k)

[f]

^^^ and [c]^^) are

Layer stresses evaluated from 3 can be rotated to the plane of the layer concerned, (l,t) and applied to the Hill (Ref.9) failure criteria, which

statesi-^l '^t It

...(4) Where (F) is the strength reserve factor for the layer concerned. For layers to be considered as safe under the applied stress system T,, TX^, T,^ this value of (F) must be in excess of unity. Should any one of the layers

fail, redistribution of stresses will occur and the laminate may still be capable of balancing the applied stress system. The capability of reacting other load cases must now be

considered as extremely suspect however and it is recommended that where one, or more of the constituent layers in a

laminate are shown to be in a failed condition, that the laminate as a whole should be declared unsafe.

Interlaminar shear stresses (-sr) result from the axisymmetric transverse shear force (Q_) (see 3.3) and are similar to the parabolic shear stress distribution associated with a rectangular sectional beam subjected to shear loading. For a laminated material these shear stresses are given by the following 1 .1,1+: 'sr k=l (See Appendix B) *^12 ds

(K,

W -

.^-^

)

d \ ds

...(5)

These stresses tend to cause delamlnation and the matrix material (normally resin) must be capable of

withstanding them.

According to an input marker a table of layer stresses is output by the programme for all elements, none of the

elements, or a selective number of elements.

4. 4.1

PROGRAMME INPUT AND OUTPUT DATA Programme Input Information

The following represents a complete list of the data required by the programme. Real variables must be punched with a decimal point and integer variables without a decimal point.

8

-The data is input in card images and in order to Indicate the number and type of variables required pei' card, the actual programmed variable names are used in the

description. If not otherwise indicated punching should start in card column one and variables should be separated by leaving a space.

Card fa^ TITLE

TITLE is an alpha-numeric title, which can occupy the first 72 columns of the first card.

Card (b) NP NE NB NLD NDF NMAT NST II (eight integer variables)

NP is the total number of nodal points NE is the total number of elements

NB is the number of constrained boundary nodes

NDF is the number of degrees of freedom per node = 3 NMAT is the total number of materials used in the

structure

NST is the number of sudden changes or steps in the wall thickness

II is an output suppression marker

If II = 1 the programme output will include the title, the values of NP, NE, NB and NMAT, a list of materials data, a list of element data and a list of constraint data.

If II = 0 all the above is suppressed and the output data will commence with the list of nodal loads.

cards (c) EL(I) ET(l) GLT(l) UMLT(l) ^ ^ ^^ ^^^

SLTT(I) SLC(I) STTC(I) SS(l)

(NMAT pairs of cards with four real variables per card)

E L ( I ) is the longitudinal modulus of elasticity of

material (l) (E,, ie in the direction of the fibres) E T ( I ) is the transverse modulus of elasticity of material

(I) (E., ie at right angles to the fibres)

9

-UMLT(I)

SLTT(I) SLC(I)

is Poisson's Ratio n,. for material (I) er.

ie

s 'i:

_E,

't "l

is the longitudinal tensile strength of material (I) is the longitudinal compressive strength of

material (l)

STTC(l) is the transverse tensile strength of material (l) S(l) is the shear strength of material (l)

All the above parameters refer to a uni-axial fibre composite in its longitudinal - transverse plane.

Cards (d) I C 0 R D ( I , M ) M = 1,2

(One integer and two real variables per card, to be input according to the following format image)

Card Columns

8_ 9 10 11 12 j.3 14 15 16 17; E ± -C0RD(I,1)

C0RD(I,2)

Integer variables are allocated a field width of four and the space must be occupied from i'ight to left, ie if I = 12, 12 must occupy column 3 and 4,

I is the node numbei"

C0RD(I,l) is the (z) co-ordinate of node I

C0RD(I,2) is the (r) co-ordinate of node I

There will be a total of (NP) cards (d) specifying the location of the (NP) nodes used in the structure

Cards (e) I MARKD(l) NOP(l,M) M = 1,2

(Four integer variables per card, to be input according to the following format image)

10

-Card Columns

1 .2 3 4 5,6,7 8 9 10, 11 ,12 13 14 15 16 17 18 19 20

MARKD(I) N0P(I,1) N0P(I,2)

is the element number

MARKD(l) is a marker for the omission of unnecessary elasticity matrix evaluation. If the element (l) is of the

same 'lay up' and the same thickness as the element (l-l) then MAFIKD(I) should be input as unity;

otherwise input as zero.

(MARKD(l) is always input as zero)

NOP(1,1) is the first node of number of element I NOP(I,2) is the second node number of element I

There will be a total of (NE) cards (e) to specify the (NE) elements used in the structure.

Cards (f) NBC(l) NFlX(l) I = 1,NB

(Two integer variables per card)

NBC(I) is the node number of the Ith constrained node NFIX(I) is the degree of constraint of the node NBC(l), •

This is a three digit number containing ones or zeros. One indicates fixation of the degree of freedom and zero indicates that the degree of

freedom is free. The first digit of the number in NFIX refers to the (u) degree of freedom, the

second to the (w) degree of freedom and the third to the (P)degree of freedom, ie if node I is fully fixed NFIX(I) = 111, whereas if it is only fixed in the radial direction NFlX(l) = 010 etc.

There will be a total of (NB) cards (f) to specify the (NB) constrained boundary nodes.

Cards (g) NS(l) I = 1,NST

(maximum of eight integer variables per card) The following five data cards (h - l) refer to the 'lay up' of the structure and a new batch of cards is

required whenever ply orientations or layer thicknesses change from element to element. In scanning through the element data (cards e) the programme expects to find a batch of cards (h - l) for every element having a 'MARKD' value of zero (see cards (e)).

11

-Card (h) NS

(single integer variable)

NS is the number of layer blocks used through the wall thickness (NS would be punched as 4 for the structure of Figure 2)

Card (1) MARK(I) I = 1,NS

(NS integer variables)

MARK(l) is a marker distinguishing the type of layer block. If MARK(l) is punched as one, then the Ith layer block has all its plies orientated at a single angle (e) to the axial direction. If MARK(l) is punched as zero, then the Ith layer block is angle plied, ie alternate plus are orientated at plus and minus (ö) to the shell axial direction. For the structure of Figure 2 MARK would be punched 0 1 1 0 .

Card (.1) MMT(I) I = 1,NS

(NS integer variables)

MMT(l) specifies the material of the Ith layer block. If MMT(l) is punched as one, then the material of the Ith layer is that defined by the first pair of cards c. If MMT(l) is punched as two, then the material of layer I is as defined by the second pair of cards c etc,

Card (k) ALLR(l) I = 1,NS (NS real variables)

A L L R ( I ) identifies the angle of orientation of the Ith layer block. For the structure of Figure 2 ALLR would be punched 60,0 90.0 90.0 -60.0.

Card (1) H(I) I = 1,NS+1

(NS+1 real variables per card)

H(l) specifies the location of layer (l) in the overall shell thickness. The elements of the array (H) are the numerical equivalents of the dimensions (h) in Figure 2. Dimensions referring to layers below the reference surface must be punched with negative sign, ie For the structure of Figure 2

(H) would be input as the numerical equivalents of -h-^ hg h^ h^ h^ (hg = 0.0),

12

-Cards (h to l) must be repeated as many times as there are zeros in the 'MARKD' array input on cards e, Cards (m) ENDL(l) I = 1,NST

(NST real variables)

E N D L ( I ) is the longitudinal end load in the shell at the location of the Ith sudden increase in wall

thickness.

(Pull circumferential values must be used, with a maximum of five variables per card)

• Cards (m) must be omitted if NST has been input as zero.

Card (n) NLN

(Single integer variable)

NLN is the total number of nodes where concentrated loads are applied.

Cards (o) NQ(l) R ( K ) K = 1,3| I = 1,NLN

(One integer and three real variables per card) NQ(l) is the node number of the Ith node where a

concentrated load is applied.

R ( 1 ) is the value of the concentrated load applied to the node NQ(l) in the axial direction (z).

R ( 2 ) is the value of the concentrated load applied to the node NQ(l) in the radial direction (r).

R ( 3 ) is the value of the concentrated bending moment . (M^) applied to the NQ(l).

N.B. It is the total circumferential values of the concentrated loads that must be fed into the programme and loads are positive in the directions of positive

displacements (see Figure l ) . If NLN is input as zero then cards (o) must be omitted,

Card (p) P

(Single real variable)

P is the internal pressure applied to the vessel. Card (q) MOUT

13

-MOUT is a marker for selective layer stress evaluation. If MOUT is punched as zero layer stresses are output for all elements.

If MOUT is punched as one no layer stresses are evaluated, or output.

If MOUT is punched as two layer stresses are output for selective elements.

If MOUT = 0 or 1 the following two cards (r and s) must be omitted.

Card (r) NLAY

(single Integer variable)

NLAY is the total number of elements requiring layer stress evaluation.

Card (s) NLOUT(I) 1 = 1 , NLAY (NLAY integer variables)

NLOUT(l) is an array of NLAY specific element numbers where layer stressing Information is required.

If information is required for more than one load case ie if NLD> 1, then cards m to s must be repeated NLD times,

4.2 Programme Output

4.2.1 Title and Salient Constants

Under the above heading the following data can be output

a) Title:- Read in as data, can be 72 characters long. b) Number of nodal points

c) Number of elements ^ d) Number of constrained boundary nodes

e) Number of materials 4.2.2 Materials Data

The above appears as a heading in the output data and a table containing the following data for all the

materials used in the structure is generated. a) Longitudinal Young's Modulus, E, b) Transverse Young's Modulus, E^

14

-c) Modulus of rigidity, G,. d) Poisson's Ratio, |x, .

e) Longitudinal tensile strength f) Longitudinal compressive strength g) Transverse tensile strength

h) Shear strength in the I - t plane

All data refers to the uni-dlrectional composite. 4.2.3 Geometry Data

Under this heading a list of nodes together with respective z and r co-ordinates is given.

4.2.4 Element Data

A list of elements with identifying node numbers. 4.2.5 Constraint Data

A list of the constrained boundary nodes together with the degree of fixation. The fixation parameter is a

three digit integer, comprising ones and zeros, denoting fixed, or free, degrees of freedom respectively. The first column location refers to the (u) degree of freedom, the second to the (w) degree of freedom and the third to the (p) degree of freedom, so that 111 denotes full fixation, 010 denotes fixation of the w degree of freedom only etc. 4.2.6 Input Load Vector

A full list of nodal loads and moments

P(U) is the load applied to the node in the (u) direction. P(V) is the load applied to the node in the (v) direction. M(S) is the applied meridional bending moment at the

node.

4.2.7 Nodal Displacements

A table of nodal displacements resulting from the applied loading of 4.2.6.

(U) is the nodal displacement in the (u) direction. (V) is the nodal displacement in the (v) direction. T H E T A is the meridional rotation angle at the node 6.

15

-4.2.8 Stress Resultants and Stress Couples

At each element the following stress resultants and stress couples are output

a) Longitudinal membrane force per unit circumferential length N(S).

b) Circumferential membrane force per unit longitudinal length N(THETA).

c) Meridional bending moment per unit circumferential length M(S).

d) Circumferential bending moment per unit longitudinal length M(THETA).

e) Transverse shear force Q(s) per unit circumferential length.

4.2.9 Layer Stresses

Where layer stressing Information is requested the following data is output.

a) Layer location:- The location of each Individual layer in the material thickness relative to the reference surface (layer 1 is the innermost layer). b) Layer stresses:- Longitudinal, transverse and shear

stresses relative to the specific (l,t) axis system of the layer concerned,

c) Strength factor:- The 'Hill' strength factor for each layer,

d) Interlaminar shear stresses:- The interlaminar

shear stress distribution resulting from the element transverse shear force (Qg),

4,3 Data Suppression

Input markers allow all data prior to para 4,2.6 and the entire table of layer stresses to be suppressed if so desired by the user.

5. DISCUSSION 5.1 Limitations

The programme listed in Appendix C has dimension statements set to accept structures of up to 100 elements. This size has proved adequate for the problems analysed to date, but may not be sufficient for some future users. It is a relatively simple matter to increase the element

16

-capacity of the programme, merely involving the modification of the relevant array sizes in the 'Common' storage areas. This however will increase the meiximum core size required to operate the programme and because of the increased number of simultaneous equations, must increase the overall computer run time. For the 100 element programme listed in Appendix C the maximum core size is approximately 25K I.CL, storage units,

The maximum semi-band width of the structural stiffness matrix is set at six and the current programme can therefore only be used for continuous shells of

revolution, where element node numbers are consecutive, The element stiffness matrix is derived for axisymmetric deformations only and reliable results can

therefore only be expected for structures conforming to this limitation. This does not necessarily limit the programme to 'balanced' laminates, as coupling between membrane and bending forces is catered for in the element forraulationi it does however mean that where angle plied blocks are used there must be as many layers at +Ö as there are at -0 so that no overall shell twisting deformations take place.

The maximum number of layer blocks catered for in the element elasticity matrix formulation is five, giving a maximum of thirty individual layers through the shell wall thickness, or thirty stations per element, where stressing information can be output. Upto five different materials (fibre reinforced, isotropic or mixed) can be used throughout the structure.

5.2 Supplementory Storage Requirements

The element elasticity matrix and several other arrays and matrices used in its derivation, are required after the solution of the nodal deflections for the

evaluation of element and more particularly layer stressing Information. To store this data in core could result in an excessively large programme and it was considered inefficient to recalculate the information in the stressing routine.

Consequently the relevant data is written to disc and the I.CL, subroutines USEFILE, PUTARRAY, GETARRAY and FREEFILE are used for this purpose. If the specific system of the user does not Incorporate these routines then programme modification will be required to access the particular disc or magnetic tape facility available.

5.3 Possible Future Developments

The current programme is intended for application to axisymmetric shells. Some structures of interest to M.O.D, have significantly large wall thicknesses however and shell assumptions may therefore no longer apply. Typical examples of such structure are nozzles and general motor end plates,

17

-where wall thicknesses can exceed one inch on diameters of six inches, or less. There is an obvious need therefore, for an element capable of predicting stress variations in a radial sense through the wall thickness. This will require the development of a new programme and work on this objective is already in progress at Cranfield,

Future developments on the existing programme depend entirely upon M,O.D, directive, but the following may be

possible,

a) Extension of the element stiffness matrix to accept thermal loading.

b) Enlargement of the stored structural stiffness matrix to cater for assemblies of axisymmetric

shells (ie where element node numbering is no longer consecutive).

c) Automatic data generation. The types of structure considered in Appendix C lend themselves very

readily to this.

Development of a new element to cater for non

symmetric deformations and,or non symmetric loading.

18

-6.

REFERENCES 1, TAYLOR, P.T. 2. TAYLOR, P.T. 3. TAYLOR, P,T, 4. FEWS, R.C. TETLOW, R. FEWS, R,C. TETLOW, R. FEWS, R.C. TETLOW, R,

7.

FEWS, R.C. TETLOW, R.

MP-ATGW Launcher Tube Study,

Cranfield Institute of Technology ATGW/COA/2 (Confidential).

Computer Programme for the Determination of Buckling

Coefficients of Fibre Reinforced Cylindrical Shells.

Cranfield Memo 102,

Cranfield Institute of Technology. Computer Programmes for the

Determination of Elastic

Properties and Strength of Thin Laminated Fibre Reinforced

Composite Plates. Cranfield Memo 92,

Cranfield Institute of Technology. Computer Programme for the

Determination of Discontinuity Forces and Couples at the Junction of a Laminated Fibre Reinforced Cylinder and Domed End Closure. College of Aeronautics Memo 7501, Cranfield Institute of Technology, October 1974.

Computer Programme for the Determination of Discontinuity Forces and Moments At a

Reinforcing Ring on a Circular Anisotropic Pressure Cylinder.

College of Aeronautics Memo 7502, Cranfield Institute of Technology, November 1974.

Computer Programme for the Stress and Strength Analysis of

Axisymmetrically Loaded Laminated Anisotropic Cylinders.

College of Aeronautics Memo 7505? Cranfield Institute of Technology, September 1975.

Design With Fibre Reinforced Materials.

College of Aeronautics Memo 7506, Cranfield Institute of Technology, March 1975.

19

-8, ZIENKIEWICZ The Finite Element Method in

Engineering Science, McGraw-Hill Book Company London. 9. HILL, R. 10. ZICKEL, J. 11. TIMOSHENKO. S,P, WINOWSKY-KRIEGER, S, 12, PEWS, R,C, TETLOW, R, 13. BRISBANE, J,J.

The Mathematical Theory of Plasticity.

Clarendon Press, Oxford, 1950. Isotensoid Pressure Vessels, ARS Journal, Vol,32, 1963, PP 950-951.

Theory of Plates and Shells, McGraw-Hill Book Company.

End Closure Shapes for Filament Wound Pressure Vessels,

Cranfield Institute of Technology, College of Aeronautics Memo 7605. Heat Conduction and Stress

Analysis of Solid Propellant Rocket Motor Nozzles,

Technical Report S-198, Rohm and Haas Company,

20

-Typical Element

Detail A

Detail A

21 -^

©_

®

•60 90 90 ±60 Reference Surface

22

-(MsVr^i^i

(MsVSj

Elements

23

-^ ^ .

^n

* ^ 1

FIG 4 GENERALISED HOOKBB LAW

24

25

-<>i

Q

*^s*.*S.s«**

26

-Meridions

a08r

007

006

-c o»

L-o

E

ê

ft> u c 10 ** _IA 5 005 004 0-03 0 0 2 001 -69

i

20 Elements O-001 M Pitch" ro co

Distance From Origin (M)-J I l _ 1 10 Elements i_ O 002 M Pitch 5 Elements __ a 0 0 4 M Pitch 0.01 0 0 2 0 0 3 0 0 4 0 0 5 FIG 9 ELEMENT GRID FOR PRESSURE VESSEL

29 -7 2 E = 10 l b / i n

r

= 0-3 5 0 " 1 lb/in 1 1.0'

M:

_^ 20 At 0 0 5 " Pitch 0 0 1 " 2.0"

{

I

2 0 ^ ^

30

-0003

_(a)

0000 FE. Programme

Analytical

Distance From Free End (in) ^ 3 04 0.5 0.6 0.7 0-8 09 1-0 0001 L-c 0 0 6 r ( b ) c £ .^ c 0; o 2 m c •o

1

c o "D ',_ _fti 2 0 0 5 0 0 4 003 002 001

0D1 L Distance From Free End

31 -0-33 d 0 32 «70 31 • * • o A0-30 S «-^ cO-29 Oi E 0» «0.28 a U) ^0-27 (0 "2 50 26 r -_ 0 -— -— 0 \ 0 \ 0 N. 0 0 0 0 025 -oooo F E. Programme Analytical Programme (Ref A)

t'8°R

216.0

2 U . 0

z

*-12

0 c 0* E

| 1 0 0

Bendin g O D 75 6 0 c o ^ 4 0 0/ 2 0 -Datum o 0 0 0 •o cr

-©-Distance From Datum (M) J I I

0005 001 0015 002 0 025 0 0 3

Distance From Daturn (M)-0 (M)-0 (M)-0 5

o o o 0025 003

Al

-APPENDIX A

ELEMENT STIFFNESS MATRIX

1. GENERAL

The element used in this analysis is a two noded axisymmetric shell element, which can be cylindrical or form a conical frustrum as shown in Figure 1. Provision is made in the formulation for three degrees of freedom per node including displacements in the axial and radial directions, (u and w) together with a meridional rotation (P). For a general shell deformation six degrees of freedom per node are required and the three degrees of freedom used in this element tailor it specifically for axisymmetric deformations. The element must therefore not be used in loading environments,

or structural configurations that result in a non axisymmetric displacement field.

2. 2.1

STRAIN DISPLACEMENT MATRIX C B ] Displacement Function

The displacement of a node (i) is defined by the following three components

tei

=

...(1)

A two noded element (i,j) therefore possesses six degrees of freedom, ie

{6/ =

...(2) Let the variation of u and w within an element be defined by u = a, + Qp s 2 3 w = a^ + a,, s + ttr- s + a/; s 3 ^^ 5 D ie in matrix form 1 s 0 0 0 0 0 0 1 s s^ s^

A2 -or 2.2

= [a] {a}

Boundary Conditions At node 1 s = 0, u = u., w w, At node j s = L, u = u., w = w. J J

...(3)

...(4)

Substituting relationships (4) into (3) and solving for the a's

gives:-ra^'^

a, a

₃

a, a,

V"6

1

0

0

0

0

0

0

1

0

.^>

0

0

0

'A

_L

2

0

VL

0

0

0

0

0

0

0

0

0

0

0

0

^/L -V:

Pn u, w,

p,

or

{_a'}

= tb]{:57

2.3 Transformation Matrix

...(5)

For the transformation of displacements from an orthogonal axis system (u,w) to an orthogonal axis system (ü,w) through a rotation 6, use is made of the following transformation matrix m - n 0 n m 0 0 0 1

Vl

r

u

where m = cos 0 n = s i n e le

{_^\ =

[T] CB}

_...(6)

A3

-2.4

Element Strain Components

With the angle / of the conical frustrum held

constant the four element strain components are given by

to =

V =

2.5

CO

du/ds

(w cos e + u sin /)/r|

d^w/d^s

- sin ^/v "^^ds

...(7)

Final Strain Displacement Matrix Formulation

From (3) and (5)

w

=

[a] [b]

{_^-f

le

= L

N]

{€f

where L N ] =

...(8)

1-s', 0 , 0 , s', 0 ,

0 , l-3s'^+2s'^, L(s'-2s'^+s'^), 0 , 3s'^-2s'^, (-s'^+s'^)L

t _ s

From (7) and (8)

C O = [B]

{j^f

= [CB],, C B ] J

{_€f

...(9)

Where the matrix [ B ] is defined by differentiating

L N ] according to the relationships stipulated at (7)

ie:-- V L

0 0

(l-s')sinx$/r ( l - 3 s ' ^ + 2 s ' ^ ) c o s / / r L ( s ' - 2 s ' ^ + s •^)cos/ii/r

[ B ] . = 0 ( - 6 + 1 2 S ' ) / L 2 ( - 4 + 6 S ' ) / L

0 (6s'-6s'2)sin)zJ/(rL) ( - l + 4 s ' - 3 s ' 2 ) sin;zJ/r

. . . ( 1 0 a )

A4

-CB]j =

l/L O O

s ' s i n ^ / r ( 3 s ' ^ - 2 s ' ^ ) c o s / / r L(-s'^+s'^)cos;zJ/r

O O

(6-12s')/L'

(-2+6s')/L

(-6s'+6s'^)sin^/(rL) (2s'-3s'^)sin/zJ/r

...(lOb)

From (6)

LU

[T] O

_ 0 [T]_

ie in terms of global co-ordinates

CO =

3.

[B]. [ T ] , [B]j [T]

ELASTICITY MATRIX D

tB}'

•^\e ...(11)

Restating equation (12) from Appendix B

_ The matrix C B 3 of Appendix B has been renamed L C -1 to avoid confusion with the strain displacement matrix [ B ] used in this Appendix.

In Appendix B equation (12) refers to an arbitrary reference surface. For the Elasticity Matrix [ D ] however, all parameters must be referred to respective shell neutral axes and the following modification is required.

CNig= [AD^ie- tEl^xi

e

...(13)

Where L E ] = [c] - L A1 [c]

[A]""*-Also from Appendix B (Equation l6)

tM\,

=

W ^ N U -

w m

se •se

For bending about the shell neutral axes this equation reduces to

A5

-For the element in question the elasticity matrix [ D ] is defined as [ D ] = [ D ] T h e r e f o r e from ( 1 3 ) a n d ( I 4 ) ' 1 1 1 2 ^ 1 1 A 21 '11 A E 22 "^21 1 2 1 1 E 21 E 22 "^21 1 2 '22 '12 '22

The matrix [ D ] takes account of possible coupling effects between membrane and bending forces, but in view of the nature of the element displacement function cannot cater for shear coupling.

4. 4.1

following

STIFFNESS MATRIX General Expression

The element stiffness matrix is defined as the

C F T

= [K]^ iU

. . . ( 1 5 )

Where {_Fj is a column matrix of nodal forces,

which are statically equivalent to the boundary stresses and distributed loads on the element. By imposing an arbitrary nodal displacement and equating external work done, with the total internal work obtained by integrating over the volume of the element it can be shown

that:-r

[ K ] ^ = 4 . 2 [B]-^ [ D ] [ B ] d ( v o l ) I n t e g r a t i o n ..(16)

For the element in question the Integrations are carried out over the element area (A) since axisymmetry effectively reduces the formulation to two dimensions.

A6 -i e dA = 27rrds = 27rrLds' X LKT = f [ B ] " ^ [ D ] [ B ] 27rrLds' O Upon s u b s t i t u t i o n , a 3 x 3 element s u b m a t r i x LK] . of t h e m a t r i x [K] ^ i s g i v e n b y : -^ T [ K ] , j = C T ? [ B ] ^ - ^ [ D ] [ B 1 ^ r d s ' > [ T ] 2TrL . . . ( 1 7 ) ( s e e e q u a t i o n s 10 and 11) Let [ B ] ^ ^ [ D ] [ B ] r = f ( s ' ) t h e n t r a n s f o r m i n g t h e l i m i t s 1

ƒ g(s')ds'

L_-l [ T 3 2TrL where g ( s ' ) = f < j | ( s ' + l ) >

1 L J

L e t I g ( s ' ) d s ' = I - 1 t h e n u s i n g f i v e p o i n t G a u s s i a n Q u a d r a t u r e 0.5689 g(0) + 0.4786 g ( 0 . 5 3 8 5 ) + 0.4786 g ( - 0 . 5 3 8 5 ) + 0.2369 s ( 0 . 9 0 6 2 ) + 0.2369 g ( - 0 . 9 0 6 2 ) r

The repetition of this process for the four 3 x 3 [ K ] . . matrices gives upon assembly the final 6 x 6 element stiffness matrix [ K ] ® .

BI

-1. 1.1

APPENDIX B

LAMINATED PLATE THEORY

TYPE OF MATERIAL

Generalised Hooks Law

With reference to Figure 4* the generalised Hooke's Law relating stresses to strains in an elastic body, can be written in matrix form as follows.

a '11 '22 '33 ^23 13 '12 '11 '12 '13 '14 '15 '12 '13 14 15 16 '22 ^23 24 ^25 26 '23 ^ 3 ^34 ^35 ^36 =24 ^34 ^44 ^45 ^46 '25 ^35 ^45 ^55 ^ 6 '16 ^26 ^36 ^46 ^56 ^66 The matrix C ij ...(1) is termed the 'material stiffness matrix' and is always symmetric, thereby containing 21 independent elastic constants.

1.2 Orthotropic Symmetry

If the structure of a material is such as to produce three mutually perpendicular planes of elastic symmetry,

then it can be shown that the stiffness matrix reduces to the following. ^ 1 1 C , 2 ^ 1 3 0 0 0 ^ 1 2 C22 ^ 2 3 0 0 0 ^ 1 3 ^ 2 3

^33

0 0 0 0 0 0

=w

0 0 0 0 0 0

c

0 55 66 ...(2)

ie the number of independent elastic constants is reduced to 9. Such a material is said to possess orthotropic sjnnmetry and is normally termed an orthotropic material.

B2

-1.3 Plane Stress

If it is assumed that a material is completely anisotropic, but very thin so that it can be considered as being two dimensional in a state of generalised plane stress; then equation (1) becomes

'11 '22 > -'12 '11 '12 '12 '22 '16 '26 '16 ^26 ^66 ^11 22 ^12 ...(3) If it is now further assumed that the material is orthotropic and the co-ordinate axes coincide with the natural axes, equation (3) reduces

to:-1^11

'22 / — '12 '11 '12 '12 '22 0 33 1 1 22 1 2 ...(4) The number of independent elastic constants is now reduced to four and it is this type of material that is used as the basic 'building block', for the composite

structure discussed in other sections of this report.

2. 2.1

MATERIAL STIFFNESS MATRIX Stress-Strain Relationships

For a uni-dlrectional composite under plane stress conditions with the fibres aligned in a direction (I)

V.'^it, ^11 h2

0

1 P PP 0 0 i '33 It ie

i

Oit = ^S\ tOit

...(5)

Where [S]^. is the local compliance matrix, which is defined simply

as: B3 as:

-Ls]i,=

' \ -^^VB.

-Hl / E .

V.

2.2 O o Axis Rotation O O

VG,.

Rotation from the (l-t) axis system, through an angle Ö, to a general axis system (x-y) is achieved by the action of the following transformation matrices

'y

m

n

n

m

mn -mn 2 2 -2mn 2mn m -n

t^U^^it

' o, m n2 mn

n

m

-2mn 2mn 2 2 -mn m -n

...(6)

...(7)

le Where m = cos e n = sin Ö ' denotes transpose.

2.3 Lamina Compliance Matrix

The compliance matrix for a unidirectional layer

rotated through an angle Q from its local co-ordinate system

is as follows.

From equations (5) and (6)

B4

-From e q u a t i o n s (7) and (8)

l e

{ - ^ i y = ^ ^ \ t O x y . . . ( 9 a ) Where [ S ] ^ = [ T ] [S]j^ [ T ] ' i s the lamina compliance

m a t r i x i n terms of g l o b a l c o - o r d i n a t e s .

2.5 S t i f f n e s s Matrix For A Multi Layered P l a t e From e q u a t i o n (9a) (k) (k) {^Oxy = t^^ t O x y . . . ( 9 b ) Where [ c] = ICsJ^ in a composite plate. (k) r ^(K:)-\"^

[ c] = \CsJp J and (k) denotes the kth layer

Consider the laminated construction shown in Figure 2. If it is assumed that strains remain constant through the plate thickness such

that:-then the total plate stresses can be written as follows. (P)

tOxy = ^^nOxy ...(10)

Where [ A ] is the effective plate stiffness matrix and is given by:-n (k) _ [CJ k^l

[A] = i

^

[C]

(h^^, -

h^)

2.6 Effective Laminate Elastic Constants From (10) it can be seen that

(P)

i'\y

= ^^HOxy ...(11)

B5

-Comparing (11) with (5) it can be seen that

P V

^x " ^11 '

xy 33

t^:

xy

= - 12/,

11

EI

= Va

22

_'yx

^2Ve

'22

3.

3.1

STRESS RESULTANTS AND STRESS COUPLES

Sign Convention

Replacing the plate co-ordinates x and y by the (s,e)

co-ordinate system as used in Appendix A to describe shells

of revolution; the stress resultants and stress couples are

positive as shown in Figure 5.

3.2

Basic Relationships

In the derivation of the following equations use is

made of an arbitrary reference surface located at some station

through the plate thickness (see Figure 2). Relative to

this reference surface the stress resultants and stress

couples can be listed as follows,

A

_{11 "^12 ^13}

'21

A

22 "^23

^31 ^32 ^33

B

₁₁

B

_{12 B}

B

B,

B

21 "22

B

B

13

23

B,

31 "^32 "33

le

tN}se=^AKOse- C B ^ ^ A e

B

₁₁

B

₂₁

B

₃₁

% 2 % 3

B22 Bg^

^32 ^33

11

12 "13

...(12)

D,

D

D,

D

21 "'22

D

D

₂₃

D.

31 32 ^33

ie

Wse= tB]0}se' t°KXi

_e

. . . ( 1 3 ) Where

C A ] =

n _{( k )}

B6

-1 f (k) p p

k=l

n

(K)

[D] = i ^ [O] (h^,^h,,5)

k=l ...(14)

3.3 Generalised Parameters

From equations (12) and (13)

The suffix (o) signifies reference surface strains.

Substituting (15) into (13) gives

Where

[b] = [B] [A]-l

1 ...(17)

[d] = L D ] - [B] [A]-l [B]

The matrix [d] represents the shell bending

stiffnesses about its neutral axes and is useful in the

definition of the 'Elasticity Matrix' for the finite element

discussed in Appendix A.

The matrix [b] gives the location of the shell

nautral axes and is also used in Appendix A for the

elasticity matrix derivation.

4. LAYER STRESSES

4.1 Strain Displacement Relationships

Assuming plane sections remain plane and that strains

through the shell thickness are small,the strain - displacement

relationships can be written as follows

bu

;N2

B7

-èu^

hv^

> > 2 ,

«, o , o ^„ o w

^se - ïïë"

•

*

" ^ i ~ " "^^ ^sïïë

Where u and v^ are the reference surface displacerneiits.

_{o o}

These equations can be written in matrix notation as follows

COse = tOo " 2 t^3se ...(18)

4.2 Strains in the kth Lamina

Substituting (15) into (18) making use of

relationships (17) yields

(k)

=

[G3

j^e I

+([b3 - z [53)

Ns

Nse

Where L G 3 = [A3~"'" and [53 is the 'Kronecker

Delta'.

4.3 Stresses in the kth Lamina

From (9b) and (I8)

,. ^ (k) ^ (k) (k) ( k ) ^ ..

Substituting for \ ^ \ Q ^^om (15)

. . . ( 2 0 )

Where [ e ] = [ C] [Gj

(k) (k)

[ f 3 = Id

[ b ] '

4.4 Interlaminar Shear Stress

In order to secure continuity of slope and radial

deflection (w) in regions of structural discontinuity

axisymmetric shells often require the action of a transverse

shear force (Q ) (see Figure 5 ) . (This force is assumed

in the programme to be positive as shown in the figure).

BB

-Where such forces exist, they will develop shear stresses ( T ) similar to the stresses developed in a simple beam under shear loading. With reference to the free body diagram shown in Figure 6 it can be seen from equilibrium of forces that

1,1+1

Tg^ dsde = 1 <;^n:^ ds)> dr

...(21) Substituting for the stresses from (20) it can be shown that

1,1+1 sr

k=l

^ L L'' ^ ^^' ' ''

fe^J

...(22)

where g-,; ' , hA ' and p,| ' are elements from the 3 x 3 matrices [g] [h] and [ p3 which are defined as follows (k) ^ (k)

LS3

=

( V l

-

\ ) [e]

,(k) (k)

Ch3 = (h^^^ - h^) [f]

(k) n P P (k)

[P3

- \

(Vl - V )

i:c3

A more detailed derivation of this equation is given in Ref.6.

5. FAILURE CRITERION

To give an indication of the strength reserves for structures analysed by the programme a failure criteria suggested by Hill Ref.9 is used to compute layer by layer strength reserve factors.

For the plane stress assumption the theorem can be stated simply as follows.

2 2 2

1

_^\ _

<^i=^t + ft_- +

''^^

F^ X^ X^ Y^ Z^ ...(23) In order to calculate the reserve factor (P) from

(23) knowledge of the layer stresses relative to their

local co-ordinate axes (l,t) is required. This information can be generated from equation (20) through the action of the transformation matrix [ T ] ' (see 2.2).

B9 -i e ( k ) ( k )

t O i t = t^^' tOse

A d e t a i l e d d e r i v a t i o n of the ' H i l l ' formula i s

shown i n R e f . 9 .

C l

-APPENDIX C

INPUT DATA GENERATION

1 . GErJERAL

The volume of data required by the finite element programme is not excessive as only one element is used through the material thickness. The input data list must however be precise and manual preparation can be both time consuming and prone to error. Items requiring definition include the overall vessel geometry, fibre orientations and stacking sequence, together with the overall shell thickness distribution. Much of this information is

available from analytical relationships, which are discussed in this Appendix and used in a small, data preparation,

computer programme,

2. 2.1

END CLOSURE SHAPES Basic Equations

The most efficient fibre reinforced composite

end closures are integrally wound with the parallel portion of the pressure vessel in a manner as indicated in Figure 7. Employing netting analysis procedure the contours can be established to give zero stress at right angles to the fibres from the following relationship (see Ref.10).

r^(z^)"

2 - tan a

z"" [l + (z*)'x (z*)'] " "" ...(1)

xvhere z* = V R , r* = V R (see Figure 7)

2.2 Geodesic Isotensoid

Before equation (l) can be solved the helix angle

(a) must be expressed in terms of the r^,7.^ co-ordinate

system. This requires prior knowledge of the winding pattern, which for most efficient head contours should be geodesic, ie

s m a

yields

...(2)

Substituting (2) into (l) and integrating both sides

(r^)^ dr^ _^

{i-(r'')2}{Cr'')^-aJ{(r«)2-aJ

1 2 + K

-i.(3)

C2 -v/here a, ag - - 2

•^m

L (^ - "-o ) J

,

-2 4r

r(-'-o^)j

+ 1

Equation (3) is an elliptic integral of the third kind and in general requires numerical solution.

2.3 Balanced 'In-Plane' Head Shapes

In practice head shapes are seldom wound to

geodesies, as this involves considerable complication in the winding machine. A more favoured procedure is to apply the fibres in a longitudinal plane tangent to the polar bosses. For this wind pattern no explicit solution to equation (l) exists as the helix angle (a) is a function of z* and r*, which can only be established for a known head contour. The desired head shape solution can be achieved hov/ever, by guessing an initial contour (ie geodesic contour) evaluating the helix angle relationship on this contour and substituting into equation (l) to give a modified solution. The helix angle relationship can now be re-evaluated and the procedure repeated until an acceptable head contour is achieved.

A full discussion of both Balanced 'In-Plane' and Geodesic head contours is given in Ref.12.

2.4 Shell Wall Thickness

Por both geodesic and balanced in plane head contours, t'lere is a progressive increase in shell wall thickness from equator to boss

ie:-t =

t cos a

o o X

r cos a

S). DATA GENERATION PROGRAMME

3.1 General

The programme listed at (4) is a modification of the head shape definition programme developed in Ref.12.

The overall pressure shell is split into a number of segments containing elements of equal axial length and the programme defines geometry, element and stacking data, in formats consistent with the input requirement of the finite element programme. Both geodesic and balanced, in-plane, head

C3

-to cylinder diameter ratiow Helix angles and thicknesses are evaluated at the central point of individual elements and data generation can commence at any location on the head contour. Different diameter polar bosses at either' end of the pressure vessel can be accommodated for the balanced, in-plane, winding method.

3.:

Description of Programme Input Data

In the description of the input data list it is Convenient to use the actual variable names appearing in the programme. Real variables must be punched with,, and integer variables without a decimal point and punching should commence on each card from card column one, separating each variable by at least one space. Card 1 LSYM

Single Integer variable

LSYM is a marker for structural symmetry. If the pressure ve.ssel is symmetrical about a central axis, ie

both end closures are the same, the finite element programme only requires definition of half the structure and LSYM must be punched as '2'. For dissimilar end closures LSYI4 must be punched as 'l'. Card 2 Rl R2 R AL TO

Five real variables

Rl is the polar boss radius at end 1. R2 is the polar boss radius at end 2.

R is the internal radius of the cylindrical portion of the pressure vessel.

TO is the equatorial wall thickness. Card 3 ZDAT RS

Tv;o real variables

ZDAT is the longitudinal distance from the origin of co-ordinates of the first node requii'ing data definition.

RS is the radial distance from the origin of co-ordinates of the first node requiring data definition.

Card 4 NS NSTA

C4

-NS is the number of discrete steps, or variations in element size required in the end closure nearest the origin of co-ordinates.

NSTA is the node number to be assigned to the first node. Cards 5 YF(l) NOINC(I) I = 1 to NS

One real and one integer variable per card., a total of (NS) cards being required, to identify the (NS) element groups in the first end closure.

YF(l) is the axial distance from datum 1 (see Figure 8) of the last node in the Ith element step range. NOINC(l) is the number of equal axial length elements in the

Ith element step range.

Referring to the worked example (3.3) there are two element step ranges in the first end closure region and hence two cards 5^ which are input as follows

10.0 5 0.0 10

Therefore the region between the start node (node 5) and an axial location 10mm from datum 1 (measured towards origin and positive in sign, see Figure 8) is split into

five equal axial steps and hence is defined by five elements. The remaining portion of the end closure (from 10mm to the datum location) is divided into 10 equal axial steps

resulting in a total of fifteen elements over the end

closure region. The element mesh is obviously finer in the second element step region than in the first.

Card 6 NCI

Single integer variable.

NCI is the number of element step ranges in the cylindrical portion of the pressure vessel.

(similar to 'NS' Card 4 ) . Card 7 P(l) Nl(l) 1 = 1 to NCI

One real and one integer variable per card; a total of NCI cards,

These cards are very similar to cards 5, but refer to the cylindrical portion of the pressure vessel Instead of the fii'st dome. Axial distances P(l) identify the Ith step range and are still measured relative to datum 1 (see Figure 8 ) . Nl(l) defines the number of elements in the Ith step range.

C5

-If LSYM (Card 1) is punched as '2' then this is all the data required by the programme. For LSYM values of '1'

however, additional data must be input to describe the second end closure and the following cards are required.

Card 8 NS

Single integer variable

NS is the number of element step ranges required. Cards 9 YF(l) NOINC(l) I = 1 to NS

One real and one integer variable per card. These are similar to Cards 5 and 7 and a total of NS cards will be required. For the last step range, ie when I = NS, YF(I) must be replaced by ( R F ) , which is the radial

co-ordinate of the last node in the structure. (There will always be a limitation in radial co-ordinate, due to the existance of a stationary point in the head contour

relationship).

3.3 Worked Example

3.3.1 Problem Description

It is required to generate data for the finite element programme describing a 90mm diameter polar wound

pressure vessel, with a wall thickness of 1mm in its parallel region. End closures are to be planar wound to polar bosses of 20mm and 40mm diameter and the vessel is to have a

parallel length of 45mm. Data must be generated from a radial location of 14mm at the 20mm boss end, to a 23mm radial location at the 40mm boss end and provision must be made for four elements to be added at the 20mm boss.

3.3.2 Programme Input Data 1 10.0 20.0 45.0 45.0 1.0 0.0 5.0 14.0 2 5 10.0 5 0.0 10

3

10.0 10 35.0 10 45.0 10 2 10.0 10 28.0 5 15.23.

C6 -3-3-3 Programme Output 5 0.50000E 01 6 0.84467E Öï 7 8 9 10 11 1 2 -13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 4 9 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 0.11893E 02 0.15340E 02 0.18787E 82 0.22233E 02 0.23233E 02 OétA233E 02 0.25233E 02 0.26233Ê 02 0.27233E 02 0.28233E 02 0.29233E 02 0.30233E 02 0.31233E 02 0.32233E 02 0.33233E 02 0,34233E 02 0.35233E 02 0.36233E 02 0.37233E 02 0.38233E 02 0O9233E 02 0.40233E 02 0.41233E 02 0.42233E 02 0.44733E 02 0.47233E 02 0.49733E 02 0.52233E 02 0.54733E 02 0.57233E 02 0.59733E 02 0.62233E 02 0.64733E 02 0.67233E 02 0.68233E 02 0,69233E 02 0.70233E 02 0,71233E 02 0.72233E 02 Ö.73233E 02 0.Ï4233E 02 to.fe233E 02 0 * ^ 2 3 3 E 02 0.77233E 02 0.78233E 02 0.79233E 02 0.80233E 02 0.81233E 02 0.82233E 02 0.83233E 02 0.84233E 02 0.85233E 02 0.86233E 02 0.87233E 02 0.9Ó254E 02 0,93274E 02 0.96295E 02 0.99316E 02 0.10234E 03 MO c.e.one.rXi DATA tJoati 1,1 z hNCi4 \ 0.14000E 02 0V28589B Ó W ^ -0*34318E 0?' 0*38158E 0 ^ O44O88OE 02 0.42797E 02 0#43231E 02 0»43609Ê 02 O,43940E 02;^ 0,44223è 02^ 0»44461E 02 0.44655E 0 2 0*44805E 02 0,44912E 02 0.44977E 02 0,45000E 02 0,45000E 02 0#45000E 02 0,45000E 02 0.45000E 02 0,45000E 02 0,45000E 02 0.450OOE 02 Ó.45000C 02 0.45000fe 02 0.45000É 02 0.45000E 02 0,450006 02 0«45000E 02 0.45000E 02 0.45000E 02 O.450OOE 02 0,45000E 02 0.45000E 02 0.450001 02 0.45000C 02 0.45000E 02 0.45000E 02 O.450OOE 02 0*45O00E 02 ^0.45000E 02 0.45000E 02 0.45000E 02 0,45000E 02 0.45000E 02 0.45000E 02 0 . 4 4 9 7 7 c 02 0.44912E 02 0.44805E 02 0.44656E 02 0.44464E 02 0.44229E 02 0.43949E 02 0.43624E 02 0.43252E 02 0.42827E 02 0.41230E 02 ©•39102E 02 0.36316E 02 0.32688E 02 0.28000E- 02 foa