- Email: [email protected]
A simple non-linear analysis of multi-layered rectangular plates
compwers & strllcfnr.?s Vol. 26. No. 4. pp. 597-M)& 1987 Printed is Great Britain.
0
004s7949/87 1987
s3.00 + 0.00 Ltd.
Pergamon Journals
A SIMPLE NON-LINEAR ANALYSIS OF MULTI-LAYERED RECTANGULAR PLATES J. T. MOTTRAM Applied Mechanics Group, School of Engineering and Applied Science, University of Durham, Science.Laboratories, South Road, Durham DHl 3LE, U.K. (Received
22 December 1986)
Al&ret-A classical, two-dimensional, geometric, non-linear, finite element is reported which predicts displacements, strains and stresses in laminated plates. The rectangular 4-noded displacement element introduced by Brebbia and Connor involves a new description of the non-linear strains, as given by Von Karman, to guarantee that all stiffness terms can be explicitly determined. Analytical solutions for isotropic and laminate test cases are used to establish limitations of the analysis. Numerical results are compared with measured surface strains and transverse displacements. From the comparison it is shown that the element is inaccurate when predicting experimental deformation for reasons given.
1.
INTRODUCTION
Development of multi-layered, fibre-reinforced composites is of interest in structural situations where high specific strength and stiffness are required. The advent of these materials has enabled further scope in design of lightweight structures [ 1,2]. A necessary priority to the successful design of composite structures is a reliable knowledge of strength and stiffness properties. To date, only uniaxial deformation has heen examined in depth[3]. Little data has been collected for more complex forms of deformation. A new biaxial test based on a plate-bending arrangement has heen devised for multi-layered laminates [4], such as used now in aircraft wings and fuselage. These plate structures are constructed from 16 or more uni-directional layers, with lamina orientations at O”, 90”, +45” and -45” to the global x-axis system. Reported here is a simple non-linear displacement, rectangular Cnoded element for generally symmetrical orthotropic plates. In particular the element was formulated to analyse the test, which required at least 36 elements in a quarter-plate model. For this reason special attention was given to reducing computing effort. The program was limited by computing resources (IBM 370/168), with a maximum array storage of 1 MB and a single job run of 1000 set CPU. To minimise computing effort, formulation of the element omitted effects of shear deformation [5], shear stresses [6], material non-linearities [I, and exact position of the neutral axis [8]. Omission of shear deformation had potentially the most serious consequence on performance. Fortunately, it has been shown that this could be neglected providing the support span to thickness ratio, S, is greater than 30 [4,9]. This meant that Kirchhoff’s assumption that normals remain straight and normal applied and only one degree of freedom (DOF) through the thickness was necessary. The geometric non-linearity was
incorporated via Von Karman strains. The basic element was that employed by Brebbia and Connor [lo] to analyse isotropic plates and shells. It will be referred to as ACMBC. The restrictions imposed by these assumptions somewhat limit the modelling ability of the element. However, Mottram and Selby [l 1] have shown that good predictions for the biaxial test are obtained with the linear (bending) element. Transition to non-linear deformation was shown to occur when the maximum central displacement exceeded the plate thickness. They also demonstrated that plate modelling was important when laminates contained lamina at +45”. Their existence provides coupling terms between bending and twisting which makes the application of quarter-plate meshes restrictive. Acceptable accuracy may therefore require a full plate model. This observation placed further emphasis on the need to formulate an efficient element, or else the governing tangential stiffness matrix would become unwieldy. The element is validated by comparing present results of the non-linear analysis for isotropic cases with those analytical solutions of Levy [12]. The comparison of computed values with experimental, obtained from the plate bending test [4], serves to show the difficulty in analysing multi-layered plates. 2. FORMULATIONOF ELEMENT
ACMBC
Finite element techniques have inherent numerical errors (e.g. lack of exact displacement representation and round-off errors), which must be moderated. To reduce further error resulting from excessive element distortion, by maintaining element aspect ratio s 3, when modelling the test, and at the same time providing nodes at positions where strains were recorded [4], a quarter-plate mesh consisted of a minimum of 36 elements. The plates to be analysed were constructed of orthotropic laminae of carbon-fibre-reinforced plastic
591
J. T. MO~~RAM
598
arranged in a generally symmetrical configuration about the mid-plane. Although a plate could be modelled by 36 elements, the number of lamina in half the plate was up to 20. Element stiffness matrices were therefore the sum of several ply contributions. The time to compute the matrices was thus important, see Spilker [13], especially since a non-linear analysis requires several iterations. It was decided that the number of nodes should be a minimum. Five DOF/node were needed to define the non-linear response [three bending (w, dw/dx and dw/dy), and two in-plane (u and v)]. Hence a rectangular element with comer nodes was used. Using this type of element the governing tangential stiffness matrix for a 36-element mesh is 245 x 245, which even though sparse, underlines the resources consideration. TO avoid numerical integration schemes when evaluating stiffness terms the element must have simple displacement functions. In-plane deformation is defined by eight DOF which can be represented by the following polynomial expressions [ 141:
Here Brebbia and Connor’s isotropic element has been furthered to cope with generally orthotropic layers, as introduced by Webber [28], whereas, their analysis terminated with the solution of nodal displacements, strains and stresses are now predicted. There is one major difference between the two formulations. To remove numerical integration schemes with ACMBC, it was decided not to use the substitution of either the full [lo, 221 or truncated form [23,29] of the bending displacement function into matrix [w]: (see Brebbia and Connor [lo, eqn 27b]). Instead the rotation terms in [WI: have been assigned values using a new description. In this work the terms dw /dx, and dw/dy, have been thought of as respective average element rotations (at the present deformation), in the x and y local directions. This leads to the following weighted means over the area of the element -=-
u=a,+a2x+a,y+a,xy
corresponding bending element has 12 DOF (15,161. The displacement function is
The
where dw/dx and dw/dy are derived from the bending function (2) and element DOF, {a},. Equation (3) can be rewritten in matrix notation as:
1 -= 1
dw = Qg dx, [
w = a, + alOx + ally + q2x2 + tl13xy
dw
dye
+ q8y3 + uIgx3y
+ a,xy3.
(2)
Combining the two linear elements provides an element (ACMBC), capable of analysing non-linear problems [lo, 171. Early work on non-linear deformation of structures paved the way for finite element plate-bending analyses. To reduce the complicated three-dimensional situation, certain assumptions had to be made, and some restrictions were placed on displacement fields. This leads to no significant loss in accuracy. Simplifications of this type produced Kirchhoff’s assumption [18] and Von Karman strain expressions [19]. Applying these approximations with classical thin plate theory assumptions [20] has proved to yield excellent analytical results, see Levy [12], and good finite element results [lo, 17,21-241. Element ACMBC was used to solve arbitrary isotropic problems by Brebbia and Connor [lo]. Their consistent formulation assumed that strains and products of rotations were negligible with respect to unity. This restriction could only be removed at great expense [25]. To produce buckling analysis the Newton-Raphson iterative procedure was used [26]. The general features of the non-linear finite element procedure are well known [27] and so will not be presented here in detail.
Qz
{6},[Ab]-’
{6}e[Ab]-’
[
Methods attacking the solution of the simultaneous algebraic equations directly are overwhelmingly the most popular. Many different forms are available, including direct iteration [30], Newton-Raphson [24, 26,3 11, incremental [32], and initial value methods [33]. All these techniques trace a load displacement response through a significant range of non-linear behaviour. Essentially they operate with intervals of loading, the solution for one being the data for the next. Numerical experience has shown that it is impractical to proceed from the initial to the final in a single step due to numerical instability. No one alogrithm has gained universal preference, as each is suited to particular situations. This is the reason why so much attention has been spent on improving this area [24,34]. ACMBC was formulated with the basic Newton-Raphson technique on each iteration. Two time-saving procedures have been taken from other analyses and introduced into ACMBC. Kawia and Yoshimura [22] presented an accurate analysis for isotropic plates using the Langrangian bending function proposed by Greene[35]. The rectangular element [four nodes (three DOF/node)] was cornbined with the constant strain element as defined by eqn (1). Matrices defining non-linear components
599
A simple non-linear analysis of multi-layered rectangular plates were large and coefficients were evaluated by numerical integration schemes. Solution of equations was performed by the novel method of Yoshiki et al. [36]. They discovered that the iteration method was often unstable. If the vector sum of internal and external forces JI {a}, (where n is the n th iteration), was calculated using only {6},_, the solution sometimes failed to converge. To avoid such a difficulty they introduced a technique which Fujino and Ohsaka [37] used to solve large deformation of suspension bridges. When choosing the assumed displacement values to be substituted into the expression of additional nodal forces, for calculation of the correction A{S In, not only {a}, _ ], but also (6)” _ z was considered. The following weighted mean was used:
in which OsB 5 1. Kawia used /I = 0.3 for large deflection problems, and a larger value of fl = 0.5 was found permissible when a system was linear. Bergan and Clough [23] developed a highly efficient formulation for isotropic plates and shells. The bending element was the doubly curved quadrilateral based on the Q-19 expressed of Clough and Fellipa [38]. When combined with the in-plane ZienkiewiczIrons iso-parametric element [21], the non-linear element had 29 DOF. The nine internal degrees were eliminated by static condensation. One feature of their procedure minimised the number of iterations when establishing convergence. In order to determine convergence the following vector was defined 6”=
46, -,-).
AS2
6 I,ref
&f
. .-
AGoor W
where A& is the change of displacement component k during iteration cycle n. Every component was then scaled by a reference displacement. For plate problems, all transverse displacement DOF were scaled by the maximum transverse displacement and similarly for rotations and in-plane terms. A non-dimensional measure of the change in displacement vector during a cycle was obtained using the maximum norm,
The convergence criterion (6 I< E could then be applied, The range for c used by Bergan and Clough was approximately 1O-2-1O-5 depending on the problem. 3. ISOTROPIC
RESULTS
The predictive performance of element ACMBC was investigated with two isotropic test cases. The plate was square with uniform pressure loading
Y A
e
x
SimpLy supported
Cbmercl
d*.“,O
d*.&“.“.O
da
ay
da
Fig. 1. Finite element 4 x 4 quarter-plate mesh.
(consistent load vector), and the bending boundary conditions were either all edges simply supported or clamped. In-plane boundary conditions could take two forms, but only Levy’s [12] restraints were treated. In this case the parallel nodal in-plane displacements to the direction of the edges are free, whilst the corresponding perpendicular DOF are fixed. Wang [39] is attributed with the second type, where all in-plane DOF along the edges are fixed. This latter condition makes the plate stiffer. Analytical results given by Levy (who solved Von Karman’s expressions using a double Fourier series), have been quoted as having a likely error ~2%. For the least number of iterations with ACMBC the value of fl in eqn (5) was always 0.5. To determine convergence using relationships (6) the terms cp, (U and D), trot (dw/dx.and dw/dy), and L, (w), were defined separately. E.+was given the smallest value, then trot, and lastly .+. This was advantageous since it had been found that a small change in Ed, caused a much larger change in the number of iterations, than if either cw or trot were equally changed, without improvement in accuracy. For the test cases examined 4, = 5 x IO-‘, cIOt= 2.5 x lo-* and cw= 10m2. The examples were analysed with a quarter plate mesh consisting of 16 uniform elements (125 DOF), as illustrated in Fig. 1. This mesh has been used as a standard in previous element evaluations [lo, 22-24,40-43]. Results are presented normalised using the parameters for transverse displacement, @ = w/h, and stress, g = aA*/Eh*. These were evaluated at load increments given by the load factor 4 = qA4/Eh4. A is the length of the sides, h plate thickness, q the uniform pressure loading, and E Young’s modulus. Poisson’s ratio v was 0.3. For the test cases, q, and the number of iterations to each convergence, are given in Table 1. Computed values at each load increment are given by symbols in Figs 2 and 3. The solid lines give analytical solutions from Levy. Stresses presented are for the tensile surface. Positions A-D and edge
J. T. MOTTRAM
600
Table 1. Load increments for isotropic test cases Simply supported No. of iterations Q 32.8 98.3 163.8 229.4 294.4
6 4 6 6 6
Clamped supported No. of iterations 4 31.3 93.8 156.2 218.8 281.3
12 10 11 12 12
boundary conditions are shown in Fig. 1. Solution for five increments to 4 = 290 took 50 + s for the simply supported, and 140 + s for the clamped example. ACMBC gave similar modelling characteristics for both isotropic examples. The value of d, (Figs 2a and 3a) was always less than analytical and the difference increased with load.
5
The computed extreme surface tensile stress at the centre, b,,o = dx,p,,D+ c?~,~,~(Figs 2b and 3b), and at the mid-side of edges, d,,, [for the clamped example (Fig. 3b)], were within 1% of analytical. The in-plane stress at the centre, dX,p,,D(Figs 2c and 3c), was always greater than analytical and vice versa for bending, bx,b,D(Figs 2d and 3d). Again as load increased so did the error. The in-plane stresses at the mid-sides, u,,~,,~ and c,,p,,e, for the simply supported case were much lower than analytical (Fig. 2~). The warping stress at comer A, 7xy,b,A,was less than analytical for 4 < 200, above which it was higher (Fig. 2d). To compare with other elements, percentage errors at load factor 4 = 200 have been collected for the test cases in Tables 2 and 3. Percentage error was determined by:
s
Fig. 2. Comparison of non-dimensional values for simply supported, isotropic (v = 0.3) square plate under uniform pressure using element ACMBC. (a) Central transverse displacement; (b) central surface tensile stress; (c) in-plane surface tensile stresses; (d) bending surface tensile stresses.
A simple non-hear
analysis of multi-layered rectangular plates
601
Fig. 3. Comparison of non-dimensional values for clamped, isotropic (v = 0.3) square plate under uniform pressure using efement ACMBC. (a) Central transverse displa~ment; (b) surface tensile stresses; (c) in-plane surface central tensile stress; (d) bending surface central tensile stress.
Analytical - finite element x 100o/ 0 analytical where superscript L means that the error was less than analytical, and H higher. Comparing ACMBC with other element types in Tables 2 and 3 it can be seen that the latest element, along with the other 4-noded elements, is inaccurate. A marked improvement is to be found with higher order elements, but at computational expense. The matrix inversion which is generally performed on every/or alternate iteration requires computing time which approximates to the cube of the total number of DOF. Assuming that all DOF in the 4 x 4 mesh were involved then the ratio of the time required with each type of element would be: Cnoded (ACMBC) 1
@noded (Serendipity) 17.58
9-noded (Lagrangian) : 39.15
Few pub~~tions reported CPU time [13,20], and when reported comparison was not practical due to different hardware. It was difficult to ascertain if any of the elements provided exact solutions. Only Bergan and Clough [23] separated Cs,,, into its components, but results for stresses were poor, see Table 2. It does appear that the predicted distribution of 6x-D, into its components with higher order elements, is not ZX*pl.D and iix,b,Dt the same as found with ACMBC (see Figs 2 and 3), since both 9, and a,,, are generally accurate, ACMBC has been seen to be poor when predicting non-linear displacements for isotropic problems. However, the principle of minimum potentiai energy appears to hold as the total stress is correct. It is suggested that the exact stress exists throughout, but that the distribution between components is incorrect. Results from ACMBC indicate that the bending stress is too low, and the in-plane too high at the centre, and vice versa towards the edges. A schematic
J. T. MOTTRAM
602
Table 2. Percentage errors with non-linear finite elements for the isotropic simply supported test case at 4 = 2M)
DQF Authors
Element (nodes)
Mottram
Analysis
4x4
6x.D
~,,,n
~x,p,.o
Simply polynomial ACMBC (4)
Laminated
125
20.0 L
l.OL
32.5 L
21.0 H
Bergan and Clough [23]
Q- 19; bending Lagrangian in-plane (4)
Isotropic
125
1.1 L
15.3 L
15.3 L
5.0 L
Pica, Wood and Hinton [40]
Heterosis QH (9) Serendipity QS (8) Langrangian QL (9)
Isotropic
425 325 405
0.3 L 0.3 L 0.3 L
0.1 H l.OH 0.9 H
Reddy and Chao [42]
Serendipity (4 linear) Serendipity (8 quartic)
125
8.4 H
1.3 L
325
1.5 L
1.3 L
325 405
l.OL l.OL
4X4X1 325
12.3 H
*‘o
Laminated
Chang and Sawamiphakdi [41]
Serendipity (8) Langrangian (9)
Laminated
Kuppusamy and Reddy 1431
3D iso-parametric (8)
Laminated
76.0 H
L, lower than analytical; H, higher than analytical.
representation of this deformation is shown in Fig. 4. The exact deformation is displayed by the solid line. The above exercise indicates that acceptable accuracy may only result with high order elements. A 4 x 4 mesh will then contain more DOF than when using 4-noded elements. A true comparison should have been based on an equivalent number of DOF. A very small improvement of 1% with ACMBC was noted, when modelling with a 6 x 6 (245 DOF) uniform element quarter-plate mesh. 4. MODIFYING
ACMBC
Brebbia and Connor reported Ku for the clamped test case. At q = 200 their central displacement differed by 23% from ACMBC, see Table 3. Now the major difference between the two finite element formulations was in the description of the rotation terms in [w]:. The expressions for dw /dx, and dw /dye given by (4) are approximate and so it should be reasonable
to take a fraction of the weighted averages. Figure 5 illustrates modelling characteristics for the simply supported example when the terms are equally scaled by 0.75-0.833. At each load increment the arrowed line shows the change as the scale factor decreased. By reducing the terms in [w]: the contribution due to non-linearity in the tangential stiffness matrix was lessened, thus increasing transverse displacements and stresses. Higher loading produced larger change. The value of &, (Fig. 5a), followed the analytical when the scaling factor was about 0.8. The value for ex,u (Fig. Sb), which was originally within I%, is now too high, since both d,,,,, (Fig. 5c) and r?X,b,D (Fig. 5d) have increased. It can also be noted that the magnitude of change was position sensitive (i.e. compare ~x,pl.Dwith ~X.,l,~). Good modelling was unattainable by just changing the two rotation terms, and in the case of (5,,n, cX,,,,,u and &,Pl,c the error was worse. This modification emphasised the original behaviour where d,,,,n was always in excess, and 6x.6,Dwas always below
Table 3. Percentage errors with non-linear finite elements for the clamped support test case at 4 = 200 Analysis
DGF 4x4
CD
=xX.0
0X.C
Simply polynomial ACMBC (4)
Laminated
125
16.0 L
l.OL
5.0 L
Brebbia and Connor [IO]
Simply polynomial ACMBC (4)
Isotropic
125
1.5 H
Thomas and Gallagher [24]
Inconsistent quad. ‘4 triangles’ (4)
Isotropic
?
4.4 H
Pica, Wood and Hinton [40]
Heterosis QH (9) Serendipity QS (8) Langrangian QL (9)
Isotropic
425 325 405
1.8 L 1.8 L 1.7L
2.0H 3.0 H 3.5H
1.4L 4.0 L 1.3L
Reddy and Chao [42]
Serendipity (4 linear) Serendipity (8 quartic)
125
18.5 H
325
4.0 H
Authors
Element (nodes)
Mottram
Laminated
603
A simple non-linear analysis of multi-layered rectangular plates Simply supported or ctamped B or C
Fig. 4. Transverse displacement along section BD isotropic test cases.
for
the analytical. It was therefore proposed that by changing the magnitude of one of the linear stiffness matrices a correction to this could be made. As the bending contribution had to increase relative to in-plane (in the central region), this could be achieved by reducing the magnitude of the bending stiffness
coefficients. Mottram and Selby [ll] have shown that the bending element predicted accurately linear measured response of multi-layered plates. This posed a restriction on scaling if it was not to lead too poor linear modelling. Scaling bending coefficients by 0.9, which was the maximum reduction permissible, produced a change at 4 = 200 (from the original ACMBC), in both the central displacement and central surface stress of < 1%. This change was insufficient to make this modification to ACMBC useful. Figure 6 shows the effect of increasing in-plane stiffness coefficients by a factor from 1.3 to 1.5. In the cases of Go (Fig. 6a), and dX,D(Fig. 6b), the change was not easily recorded. do was now less than that obtained with basic element. Of relevance from this exercise was the fact that cX,.owas now less than analytical. This occurred because 6,,,,n had substan-
(b) 200
-
si
00
100 -
0
csD
f
tiAL__ IO
20
3’0
)-
Fig. 5. Comparison of non-dimensional values for simply supported, isotropic (v = 0.3) square plate under uniform presure when scaling non-linear stiffness coefficients of element ACMBC from 0.75 to 0.833. (a) Central transverse displacement; (h) central surface tensile stress; (c) in-plane surface tensile stresses; (d) bending surface tensile stresses.
604
J. T. MOITRAM
-
Levy
Fig. 6. Comparison of non-dimensional values for simply supported, isotropic (v = 0.3) square plate under uniform pressure when scaling in-plane stiffness coefficients of element ACMBC from 1.3 to 1.5. (a) Central transverse displacement; (b) central surface tensile stress; (c) in-plane surface tensile stresses; (d) bending surface tensile stresses.
tially decreased (cf. Figs 6c and 2c), while bx,b,Dhad
5. LAMINATE RESULTS
slightly increased (cf. Figs 6d and 2d). Good modelling was also not possible by just scaling in-plane stiffnesses. A successful model was suggested through a combination of reducing the rotation terms and increasing the in-plane stiffnesses. By application of a trial and error procedure this was found when dw/dx, and dw/dy, were scaled by 0.75 and in-plane coefficients by 1.3. Figure 7 illustrates the improvement for the simply supported example. Only stresses bX,P,,c(Fig. 7c), and fX,,4A(Fig. 7d), were not within 1% of analytical. The same was found for the clamped case, but for conciseness has not been reproduced here. Comparing the new ACMBC with other elements in Table 2, it can be seen that it now has similar performance to higher order elements, but at a much lower cost.
TO try to find out if the scaling factors were universal for all types of problems, and especially when the plates consisted of generally orthotropic layers, further test cases were sought. Suitable examples were somewhat scarce, since of the available ones, Chang and Sawamiphakdi [41] studied a single honeycomb sandwich problem and both Reddy et al. [42,43] and Noor and Mathers [44] favoured cross-ply and anti-symmetrical problems. Noor and Mathers did study two symmetrical quasi-isotropic laminates but like Reddy were concerned with shear deformation in thick plates. Recently, Noor [45], using a two-step hybrid analytical technique, solved an ideal problem. The plate was square, all edges were clamped (u = u =
A simple non-linear analysis of multi-layered rectangular plates
605
[ I
I
I
.
Levy
-
l
$
modified
200
-
;i
c
IO
2C
20 s
0
DC
DA
200-
4
5
15
IO
Fig. 7. Comparison of non-dimensional values for simply supported, isotropic (v = 0.3) square plate under uniform pressure using the modified element ACMBC. (a) Central transverse displacement; (b) central surface tensile stress; (c) in-plane surface tensile stresses; (d) bending surface tensile stresses.
w = dw/dx = dw/dy = 0), and the loading was uniform pressure. It consisted of 16 orthotropic layers, with the symmetrical lay-up (45”, -45”, 0”, 0”, -45”, 45”, 90”, go”),. Material properties of each lamina were: E 1 l/E22 = 10.0, G 12/E22 = 0.49, v 12 = 0.38. Shown in Fig. 8 are computed values for central displacement, $, using a 6 x 6 full-plate mesh. The non-dimensional load factor was given by 4 = qS4/E22. A value of S = 120 was used by Noor so his results can be those expected with thin plate assumptions. The modified ACMBC element was found to predict nearly identical values to the hybrid technique and a finite element reported by Chia [46]. The very small discrepancies between the three approaches was very encouraging as they demonstrated that the modified ACMBC could analyse laminates, and in particular the biaxial test. Biaxial tests were performed on a range of laminated rectangular plate specimens [4]. In the test the specimen was supported on four steel ball supports, C.A.S. 2414-E
near to each corner. Patch loading was applied centrally through a square steel loading bar with a rubber foot pad. Measurements of central transverse
zcool
15m-
i
l
Hybrid technique (NOOR 45) Finite element (CHIP, 46) Modified ACMBC
IOCXJ -
500-
a
Fig. 8. Comparison of non-dimensionalised central displacement for clamped square 164ayer laminate plate under uniform pressure.
J. T. Mo-rra~~
606
0
0 Point support
I
.
B
6
:-
-
Fig. 9. Test arrangement and 36-element quarter-plate mesh for a 34-layered specimen.
displacement and strains were taken after each load increment until first fibre failure. Illustrated in Fig. 9 is the test arrangement and a 3delement quarterplate model for a 3Clayered specimen. Parameters for the test and modified stiffnesses[ll] are: lay-up (90”, O”, 90”, 45”, O”, -45”, O”, O”, -45”, O”, 45”, O”, o”, 45”, 0”, -45”, Oo), A = B = 260 mm, A, = B, = 200 mm, A, = BP = 20 mm, h = 4.06 mm; lamina material properties: El 1 = 0.149 x 10” N/m*, E22 = 0.289 x 10” N/m*, G12 = 0.577 x 10” N/m*, v 12 = 0.3, t = 0.1194mm. For this specimen S was 49.3, and so thin plate analysis was applicable. Figure 10 presents non-dimensional numerical values with corresponding measurements. Results from ACMBC were the average of the two quarter-
plate models which compensate for warping [l 11. To solve .five increments of load to Q = 18350, 970 set of CPU were used. The load factor is given by Q = qpAi/E22h4, where q,,= P/A,B, is the pressure over patch area. P is the applied load. After completion of the program the IBM 370/168 had been replaced by an AMDAHL 47O/V8. With this computer nine increments to q = 30550 took 740 sec. Results from the basic element have not been given because the analysis went numerically unstable at 4 = 14,250. Up to q = 22,400 the modified ACMBC can be seen to give a comparable picture to that found for isotropic problems with the initial element, see Figs 2 and 3. At this load the measured central transverse displacement, EC, was underestimated by
I
. -
l
0
,
,
Experiment Modified ACMBC Linear ACMBC
.
. . .
.
. .
.
.o .
.
0
.
. .O .. . . 2 (lo-3)
Fig. 10. Expezimental-computed
4
6
6
IO
12
(YC
comparison with modified element ACMBC of non-linear behaviour for a 34-symmetrical layered square plate.
14
.
A simple non-linear analysis of multi-layered rectangular plates
over lOO%, whilst the error for central surface strains &, and & (and thus stresses), was very small. For loading above 4 = 22,400 the numerical strain in the x global direction grew at a much higher rate than measured, and vice versa in the y direction. Over the same load range the computed displacement increased at a much slower rate than that measured. It is only speculation, but from other tests on a range of plates with various lay-ups and test parameters [4], it is suggested that perhaps the non-linear behaviour of the laminates has contributing factors not allowed for by the assumptions. These additional factors could involve effect of support movement, material non-linearity, position of neutral axis, geometric nonlinearity, resin failure under tension, and preferential bending about one of the global axes [4]. The comparison in Fig. 10 indicates that the effect of these became dominant when the maximum displacement exceeded about three plate thicknesses. 6. CONCLUSIONS
A highly efficient, non-linear, laminate, finite element is developed based on classical thin-plate theory and Von Karman strain expressions. The simple displacement 4-noded rectangular element is referred to as ACMBC, after Adini and Clough, and Melosh who formulated the linear element, and Brebbia and Connor who developed the isotropic non-linear element. A great amount of effort was spent formulating a powerful element capable of solving quarter-plate models consisting of a least 36 elements and up to 20 layers in half the thickness, without excessive need for computing resources. To this end, numerical integration schemes for evaluation of all stiffness coefficients were removed by a new description of element rotations. Computed results are presented for isotropic test cases. Good modelling was unobtainable using the basic element. Central transverse displacement was greatly underestimated. Total stresses (or strains) were excellent, but the distribution between bending and in-plane components was very poor. This was corrected by scaling non-linear and in-plane stiffness coefficients. To demonstrate that the scaling factors were universal a laminate test case was employed. Results using the modified ACMBC element were compared with measurement from bending of a rectangular multi-layered plate. The investigation showed that computed solutions gave a response for laminates up to a central displacement of three plate thicknesses similar to that achieved for isotropic problems with the basic element. For central displacements above this, additional factors are believed to result in larger numerical error. REFEBENCJB 1. J. V. Noyes,
Composites in the construction of the Lear Fan 2100 aircraft. Composites 14, 129-139 (1983).
607
2. E. Fitzer, Carbon based composites. J. Chim. Phys. 81, 717-733 (1984). 3. S. W. Tsai and H. T. Hahn, Introduction to Composite Materials. Technomic, Stamford, CT (1980). 4. J. T. Mottram, The failure of carbon fibre reinforced laminated plates under biaxial stresses. Ph.D. thesis, Department of Engineering Science, Durham University (1984). 5. N. J. Pagano, Exact solution for rectangular bidirectional composites and sandwich plates. J. Comp. Muter. 4, 20-34 (1970). 6. N. J. Pagano, Stress fields in composite laminates. ht. J. Solids Struct. 4, 38%$00 (1978). 7. H. T. Hahn, Non-linear behaviour of laminated composites. J. Comp. Mater. 7, 257-271 (1973). 8. R. M. Jones and H. S. Morgan, Bending and extension of cross-ply laminates with different moduli in tension and compression. Comput. Struct. 11, 181-191 (1980). 9. A. K. Noor and S. J. Hartley, Effect of shear deformation and anisotropy of the non-linear response of composite plates. In Development of Composite Materials (Edited by G. Holister), pp. 55-65. Applied Science, Barking, U.K. (1977). 10. C. Brebbia and J. Connor, Geometric nonlinear finite element analysis. J. Engng Mech. Div., Proc. Am. Sot. Civ. Engrs 95 (EM2), 463483 (1969). 11. J. T. Mottram and A. R. Selby, Bending of thin laminated plates. Comp. Struct. 25, 271-280 (1987). 12. S. Levy, Large deflection theory for rectangular plates. Proc. Symp. Appl. Math AMS, 1, 197 (1969). 13. R. L. Spilker and D. M. Jakobs, Hybrid stress reducedMindlin elements for thin multi-layered plates. Int. J. Namer. Methoa!s Engng 23, 555-578 (1986). 14. K. C. Rockey, H. R. Evans, D. W. Griffiths and D. A. Nethercot, The Finite Element Method. A Basic Introduction to Engineers. Crosby Lockwood Staples (1970). 15. R. J. Melosh. Basic derivation of matrices for the direct stiffness method. J. AIAA 1, 1631-1637 (1963). 16. A. Adini and R. W. Clough, Analysis of plate bending by the finite element method and report. Natl. SC; Foundation USA. G7337 (1961). 17. T. M. Roberts and D. G. Ashwell, The use of tinite element mid-increment stiffness matrices in the postbuckling analysis of imperfect structures. Int. J. Sofiak Struct. 7, 805-823 (1971). 18. G. Kirchoff, Uber das Gliechgewicth und die bewegungeiner elastischen scheile. J. Reine Angew. Marh. 40, 51-88 (1850). 19. Y. C. Fung, Foundation of Solid Mechunics. PrenticeHall, Englewood Cliffs, NJ (1965). 20. J. E. Ashton and J. M. Whitney, Theory of laminated plates. Progress in Material Science, Vol. IV. Technomic, Stamford, CT (1970). 21. 0. C. Zienkiewicz, The Finite Element in Structural Continuum Mechanics. McGraw-Hill, London (1976). 22. T. Kawia and N. Yoshimura, Analysis of large deflection of plates by finite element method. Int. J. Numer. Methods Engng 1, 123-133 (1969). 23. P. G. Beraan and R. W. Clot&. Lame deflection analysis ofl plates and shallow she& us&g the tinite element method. Int. J. Numer. Method Engng 5, 534-556 (1973). 24. R. G. Thomas and R. H. Gallagher, A triangular thin shell finite element: nonlinear analvsis. NASA CR-2483 (1975). 25. C. Prato, A mixed finite element for thin shell analysis. Ph.D. thesis, Department of Civil Engineering, Massachusetts Institute of Technology (1968). 26. R. H. Gallagher and J. Padlog, Discrete element approach to structural stability analysis. J. AIAA 1, 1437-1439 (1963). 27. 0. C. Zienkiewicz, The Finite Element Method in Engineering Science. McGraw-Hill, London (1971).
J. T. MOTTRAM
608
28. J. P. H. Webber, Experimental and theoretical bending deflections of laminated plates. J. Struct. Anal. 7, 87-96 (1972). 29. R. H. Gallagher, R. A. Gellatly, J. Padlog and R. H. Mallett, A discrete element procedure for thin shell instability analysis. J. AIAA 5, 138-145 (1967). 30. S. Mau and R. H. Gallagher, A finite element procedure for nonlinear pre-buckling and initial postbuckling analysis. NASA CR-1936 (1972). 31. J. T. Gden, Numerical formulation of nonlinear elasticity problems. Proc. Am. Sot. Ciu. Engrs 93(ST3), 235-255 (1967).
32. L. D. Hofmeister, G. Greenbaum and D. Evenson, Large strain elastic plastic finite element analysis. J. AZAA 9, 1248-1254 (1971). 33. R. D. Richtmeyer and K. N. Morton, D@rence Methoak for Initial Value Problems. Interscience, New York (1967). 34. R. H. Gallagher, Finite element analysis of geometrically nonlinear problems. In Theory and Practice in Finite Elemenr Structural Analysis (Edited by Yamada), pp. 109-125. SEM 74/32064 (1973). 35. B. E. Greene, Stiffness matrix for bending of a rectangular element with initial membrane stresses. Structural Analysis Research Memo. no. 45, The Boeing Co., Seattle (1962). 36. M. Yoshiki, T. Kawia and N. Yoshimura, Matrix methods of analysis of ship structures, 3rd Report. J. Sot. Nau. Arch. Japan 123, 188-195 (1968). 31. T. Fujino and K. Ohsaka, Static structural analysis of
suspension bridges. Mitsubishi Nippon Heavy-Znd Tech. Rev. 3, 17-23 (1966). 38. R. W. Clough and C. A. Fellipa, A refined quadrilateral element for-analysis of plate-bending. Pro;. &d Conf. Matrix Methods in Structural Mechanics. AFFDL-TR68-150, Wright Patterson AFB, OH, ‘pp. 399-440, (1968). 39. C. T. Wang, Non-linear large deflection boundary value problems of rectangular plates. NACA TN 1425 (1948). 40. A. Pica, R. D. Wood and E. Hinton, Finite element analysis of geometrically non-linear plate bending behaviour using a Mindlin formulation. Compur. Struct. 11, 203-215 (1980). 41. T. Y. Chang and K. Sawamiphakdi, Large deformation analysis of laminated shells by the finite element method. Comput. Srruct. 13, 331-340 (1981). 42. J. N. Reddy and W. C. Chao, Large deflection and large amplitude free vibration of laminated composite material plates. Compur. Srruct. 13, 341-347 (1981). 43. T. Kuppusamy and J. N. Reddy, A three-dimensional nonlinear analysis of cross-ply rectangular composite plates. Comput. Struct. 18, 263-272 (1984). 44. A. K. Noor and M. D. Mathers. Nonlinear finite element analysis of laminated composite shells. Computing Methods in Nonlinear Mechanics. Inst. Comp. Mechs, Texas, pp. l-999 (1974). 45. A. K. Noor, Hybrid analytical technique for nonlinear analysis of structures. J. AZAA 23, 938-946 (1985). 46. C. Y. Chia, Nonlinear analysis of plates, McGraw-Hill, New York (1980).
0
004s7949/87 1987
s3.00 + 0.00 Ltd.
Pergamon Journals
A SIMPLE NON-LINEAR ANALYSIS OF MULTI-LAYERED RECTANGULAR PLATES J. T. MOTTRAM Applied Mechanics Group, School of Engineering and Applied Science, University of Durham, Science.Laboratories, South Road, Durham DHl 3LE, U.K. (Received
22 December 1986)
Al&ret-A classical, two-dimensional, geometric, non-linear, finite element is reported which predicts displacements, strains and stresses in laminated plates. The rectangular 4-noded displacement element introduced by Brebbia and Connor involves a new description of the non-linear strains, as given by Von Karman, to guarantee that all stiffness terms can be explicitly determined. Analytical solutions for isotropic and laminate test cases are used to establish limitations of the analysis. Numerical results are compared with measured surface strains and transverse displacements. From the comparison it is shown that the element is inaccurate when predicting experimental deformation for reasons given.
1.
INTRODUCTION
Development of multi-layered, fibre-reinforced composites is of interest in structural situations where high specific strength and stiffness are required. The advent of these materials has enabled further scope in design of lightweight structures [ 1,2]. A necessary priority to the successful design of composite structures is a reliable knowledge of strength and stiffness properties. To date, only uniaxial deformation has heen examined in depth[3]. Little data has been collected for more complex forms of deformation. A new biaxial test based on a plate-bending arrangement has heen devised for multi-layered laminates [4], such as used now in aircraft wings and fuselage. These plate structures are constructed from 16 or more uni-directional layers, with lamina orientations at O”, 90”, +45” and -45” to the global x-axis system. Reported here is a simple non-linear displacement, rectangular Cnoded element for generally symmetrical orthotropic plates. In particular the element was formulated to analyse the test, which required at least 36 elements in a quarter-plate model. For this reason special attention was given to reducing computing effort. The program was limited by computing resources (IBM 370/168), with a maximum array storage of 1 MB and a single job run of 1000 set CPU. To minimise computing effort, formulation of the element omitted effects of shear deformation [5], shear stresses [6], material non-linearities [I, and exact position of the neutral axis [8]. Omission of shear deformation had potentially the most serious consequence on performance. Fortunately, it has been shown that this could be neglected providing the support span to thickness ratio, S, is greater than 30 [4,9]. This meant that Kirchhoff’s assumption that normals remain straight and normal applied and only one degree of freedom (DOF) through the thickness was necessary. The geometric non-linearity was
incorporated via Von Karman strains. The basic element was that employed by Brebbia and Connor [lo] to analyse isotropic plates and shells. It will be referred to as ACMBC. The restrictions imposed by these assumptions somewhat limit the modelling ability of the element. However, Mottram and Selby [l 1] have shown that good predictions for the biaxial test are obtained with the linear (bending) element. Transition to non-linear deformation was shown to occur when the maximum central displacement exceeded the plate thickness. They also demonstrated that plate modelling was important when laminates contained lamina at +45”. Their existence provides coupling terms between bending and twisting which makes the application of quarter-plate meshes restrictive. Acceptable accuracy may therefore require a full plate model. This observation placed further emphasis on the need to formulate an efficient element, or else the governing tangential stiffness matrix would become unwieldy. The element is validated by comparing present results of the non-linear analysis for isotropic cases with those analytical solutions of Levy [12]. The comparison of computed values with experimental, obtained from the plate bending test [4], serves to show the difficulty in analysing multi-layered plates. 2. FORMULATIONOF ELEMENT
ACMBC
Finite element techniques have inherent numerical errors (e.g. lack of exact displacement representation and round-off errors), which must be moderated. To reduce further error resulting from excessive element distortion, by maintaining element aspect ratio s 3, when modelling the test, and at the same time providing nodes at positions where strains were recorded [4], a quarter-plate mesh consisted of a minimum of 36 elements. The plates to be analysed were constructed of orthotropic laminae of carbon-fibre-reinforced plastic
591
J. T. MO~~RAM
598
arranged in a generally symmetrical configuration about the mid-plane. Although a plate could be modelled by 36 elements, the number of lamina in half the plate was up to 20. Element stiffness matrices were therefore the sum of several ply contributions. The time to compute the matrices was thus important, see Spilker [13], especially since a non-linear analysis requires several iterations. It was decided that the number of nodes should be a minimum. Five DOF/node were needed to define the non-linear response [three bending (w, dw/dx and dw/dy), and two in-plane (u and v)]. Hence a rectangular element with comer nodes was used. Using this type of element the governing tangential stiffness matrix for a 36-element mesh is 245 x 245, which even though sparse, underlines the resources consideration. TO avoid numerical integration schemes when evaluating stiffness terms the element must have simple displacement functions. In-plane deformation is defined by eight DOF which can be represented by the following polynomial expressions [ 141:
Here Brebbia and Connor’s isotropic element has been furthered to cope with generally orthotropic layers, as introduced by Webber [28], whereas, their analysis terminated with the solution of nodal displacements, strains and stresses are now predicted. There is one major difference between the two formulations. To remove numerical integration schemes with ACMBC, it was decided not to use the substitution of either the full [lo, 221 or truncated form [23,29] of the bending displacement function into matrix [w]: (see Brebbia and Connor [lo, eqn 27b]). Instead the rotation terms in [WI: have been assigned values using a new description. In this work the terms dw /dx, and dw/dy, have been thought of as respective average element rotations (at the present deformation), in the x and y local directions. This leads to the following weighted means over the area of the element -=-
u=a,+a2x+a,y+a,xy
corresponding bending element has 12 DOF (15,161. The displacement function is
The
where dw/dx and dw/dy are derived from the bending function (2) and element DOF, {a},. Equation (3) can be rewritten in matrix notation as:
1 -= 1
dw = Qg dx, [
w = a, + alOx + ally + q2x2 + tl13xy
dw
dye
+ q8y3 + uIgx3y
+ a,xy3.
(2)
Combining the two linear elements provides an element (ACMBC), capable of analysing non-linear problems [lo, 171. Early work on non-linear deformation of structures paved the way for finite element plate-bending analyses. To reduce the complicated three-dimensional situation, certain assumptions had to be made, and some restrictions were placed on displacement fields. This leads to no significant loss in accuracy. Simplifications of this type produced Kirchhoff’s assumption [18] and Von Karman strain expressions [19]. Applying these approximations with classical thin plate theory assumptions [20] has proved to yield excellent analytical results, see Levy [12], and good finite element results [lo, 17,21-241. Element ACMBC was used to solve arbitrary isotropic problems by Brebbia and Connor [lo]. Their consistent formulation assumed that strains and products of rotations were negligible with respect to unity. This restriction could only be removed at great expense [25]. To produce buckling analysis the Newton-Raphson iterative procedure was used [26]. The general features of the non-linear finite element procedure are well known [27] and so will not be presented here in detail.
Qz
{6},[Ab]-’
{6}e[Ab]-’
[
Methods attacking the solution of the simultaneous algebraic equations directly are overwhelmingly the most popular. Many different forms are available, including direct iteration [30], Newton-Raphson [24, 26,3 11, incremental [32], and initial value methods [33]. All these techniques trace a load displacement response through a significant range of non-linear behaviour. Essentially they operate with intervals of loading, the solution for one being the data for the next. Numerical experience has shown that it is impractical to proceed from the initial to the final in a single step due to numerical instability. No one alogrithm has gained universal preference, as each is suited to particular situations. This is the reason why so much attention has been spent on improving this area [24,34]. ACMBC was formulated with the basic Newton-Raphson technique on each iteration. Two time-saving procedures have been taken from other analyses and introduced into ACMBC. Kawia and Yoshimura [22] presented an accurate analysis for isotropic plates using the Langrangian bending function proposed by Greene[35]. The rectangular element [four nodes (three DOF/node)] was cornbined with the constant strain element as defined by eqn (1). Matrices defining non-linear components
599
A simple non-linear analysis of multi-layered rectangular plates were large and coefficients were evaluated by numerical integration schemes. Solution of equations was performed by the novel method of Yoshiki et al. [36]. They discovered that the iteration method was often unstable. If the vector sum of internal and external forces JI {a}, (where n is the n th iteration), was calculated using only {6},_, the solution sometimes failed to converge. To avoid such a difficulty they introduced a technique which Fujino and Ohsaka [37] used to solve large deformation of suspension bridges. When choosing the assumed displacement values to be substituted into the expression of additional nodal forces, for calculation of the correction A{S In, not only {a}, _ ], but also (6)” _ z was considered. The following weighted mean was used:
in which OsB 5 1. Kawia used /I = 0.3 for large deflection problems, and a larger value of fl = 0.5 was found permissible when a system was linear. Bergan and Clough [23] developed a highly efficient formulation for isotropic plates and shells. The bending element was the doubly curved quadrilateral based on the Q-19 expressed of Clough and Fellipa [38]. When combined with the in-plane ZienkiewiczIrons iso-parametric element [21], the non-linear element had 29 DOF. The nine internal degrees were eliminated by static condensation. One feature of their procedure minimised the number of iterations when establishing convergence. In order to determine convergence the following vector was defined 6”=
46, -,-).
AS2
6 I,ref
&f
. .-
AGoor W
where A& is the change of displacement component k during iteration cycle n. Every component was then scaled by a reference displacement. For plate problems, all transverse displacement DOF were scaled by the maximum transverse displacement and similarly for rotations and in-plane terms. A non-dimensional measure of the change in displacement vector during a cycle was obtained using the maximum norm,
The convergence criterion (6 I< E could then be applied, The range for c used by Bergan and Clough was approximately 1O-2-1O-5 depending on the problem. 3. ISOTROPIC
RESULTS
The predictive performance of element ACMBC was investigated with two isotropic test cases. The plate was square with uniform pressure loading
Y A
e
x
SimpLy supported
Cbmercl
d*.“,O
d*.&“.“.O
da
ay
da
Fig. 1. Finite element 4 x 4 quarter-plate mesh.
(consistent load vector), and the bending boundary conditions were either all edges simply supported or clamped. In-plane boundary conditions could take two forms, but only Levy’s [12] restraints were treated. In this case the parallel nodal in-plane displacements to the direction of the edges are free, whilst the corresponding perpendicular DOF are fixed. Wang [39] is attributed with the second type, where all in-plane DOF along the edges are fixed. This latter condition makes the plate stiffer. Analytical results given by Levy (who solved Von Karman’s expressions using a double Fourier series), have been quoted as having a likely error ~2%. For the least number of iterations with ACMBC the value of fl in eqn (5) was always 0.5. To determine convergence using relationships (6) the terms cp, (U and D), trot (dw/dx.and dw/dy), and L, (w), were defined separately. E.+was given the smallest value, then trot, and lastly .+. This was advantageous since it had been found that a small change in Ed, caused a much larger change in the number of iterations, than if either cw or trot were equally changed, without improvement in accuracy. For the test cases examined 4, = 5 x IO-‘, cIOt= 2.5 x lo-* and cw= 10m2. The examples were analysed with a quarter plate mesh consisting of 16 uniform elements (125 DOF), as illustrated in Fig. 1. This mesh has been used as a standard in previous element evaluations [lo, 22-24,40-43]. Results are presented normalised using the parameters for transverse displacement, @ = w/h, and stress, g = aA*/Eh*. These were evaluated at load increments given by the load factor 4 = qA4/Eh4. A is the length of the sides, h plate thickness, q the uniform pressure loading, and E Young’s modulus. Poisson’s ratio v was 0.3. For the test cases, q, and the number of iterations to each convergence, are given in Table 1. Computed values at each load increment are given by symbols in Figs 2 and 3. The solid lines give analytical solutions from Levy. Stresses presented are for the tensile surface. Positions A-D and edge
J. T. MOTTRAM
600
Table 1. Load increments for isotropic test cases Simply supported No. of iterations Q 32.8 98.3 163.8 229.4 294.4
6 4 6 6 6
Clamped supported No. of iterations 4 31.3 93.8 156.2 218.8 281.3
12 10 11 12 12
boundary conditions are shown in Fig. 1. Solution for five increments to 4 = 290 took 50 + s for the simply supported, and 140 + s for the clamped example. ACMBC gave similar modelling characteristics for both isotropic examples. The value of d, (Figs 2a and 3a) was always less than analytical and the difference increased with load.
5
The computed extreme surface tensile stress at the centre, b,,o = dx,p,,D+ c?~,~,~(Figs 2b and 3b), and at the mid-side of edges, d,,, [for the clamped example (Fig. 3b)], were within 1% of analytical. The in-plane stress at the centre, dX,p,,D(Figs 2c and 3c), was always greater than analytical and vice versa for bending, bx,b,D(Figs 2d and 3d). Again as load increased so did the error. The in-plane stresses at the mid-sides, u,,~,,~ and c,,p,,e, for the simply supported case were much lower than analytical (Fig. 2~). The warping stress at comer A, 7xy,b,A,was less than analytical for 4 < 200, above which it was higher (Fig. 2d). To compare with other elements, percentage errors at load factor 4 = 200 have been collected for the test cases in Tables 2 and 3. Percentage error was determined by:
s
Fig. 2. Comparison of non-dimensional values for simply supported, isotropic (v = 0.3) square plate under uniform pressure using element ACMBC. (a) Central transverse displacement; (b) central surface tensile stress; (c) in-plane surface tensile stresses; (d) bending surface tensile stresses.
A simple non-hear
analysis of multi-layered rectangular plates
601
Fig. 3. Comparison of non-dimensional values for clamped, isotropic (v = 0.3) square plate under uniform pressure using efement ACMBC. (a) Central transverse displa~ment; (b) surface tensile stresses; (c) in-plane surface central tensile stress; (d) bending surface central tensile stress.
Analytical - finite element x 100o/ 0 analytical where superscript L means that the error was less than analytical, and H higher. Comparing ACMBC with other element types in Tables 2 and 3 it can be seen that the latest element, along with the other 4-noded elements, is inaccurate. A marked improvement is to be found with higher order elements, but at computational expense. The matrix inversion which is generally performed on every/or alternate iteration requires computing time which approximates to the cube of the total number of DOF. Assuming that all DOF in the 4 x 4 mesh were involved then the ratio of the time required with each type of element would be: Cnoded (ACMBC) 1
@noded (Serendipity) 17.58
9-noded (Lagrangian) : 39.15
Few pub~~tions reported CPU time [13,20], and when reported comparison was not practical due to different hardware. It was difficult to ascertain if any of the elements provided exact solutions. Only Bergan and Clough [23] separated Cs,,, into its components, but results for stresses were poor, see Table 2. It does appear that the predicted distribution of 6x-D, into its components with higher order elements, is not ZX*pl.D and iix,b,Dt the same as found with ACMBC (see Figs 2 and 3), since both 9, and a,,, are generally accurate, ACMBC has been seen to be poor when predicting non-linear displacements for isotropic problems. However, the principle of minimum potentiai energy appears to hold as the total stress is correct. It is suggested that the exact stress exists throughout, but that the distribution between components is incorrect. Results from ACMBC indicate that the bending stress is too low, and the in-plane too high at the centre, and vice versa towards the edges. A schematic
J. T. MOTTRAM
602
Table 2. Percentage errors with non-linear finite elements for the isotropic simply supported test case at 4 = 2M)
DQF Authors
Element (nodes)
Mottram
Analysis
4x4
6x.D
~,,,n
~x,p,.o
Simply polynomial ACMBC (4)
Laminated
125
20.0 L
l.OL
32.5 L
21.0 H
Bergan and Clough [23]
Q- 19; bending Lagrangian in-plane (4)
Isotropic
125
1.1 L
15.3 L
15.3 L
5.0 L
Pica, Wood and Hinton [40]
Heterosis QH (9) Serendipity QS (8) Langrangian QL (9)
Isotropic
425 325 405
0.3 L 0.3 L 0.3 L
0.1 H l.OH 0.9 H
Reddy and Chao [42]
Serendipity (4 linear) Serendipity (8 quartic)
125
8.4 H
1.3 L
325
1.5 L
1.3 L
325 405
l.OL l.OL
4X4X1 325
12.3 H
*‘o
Laminated
Chang and Sawamiphakdi [41]
Serendipity (8) Langrangian (9)
Laminated
Kuppusamy and Reddy 1431
3D iso-parametric (8)
Laminated
76.0 H
L, lower than analytical; H, higher than analytical.
representation of this deformation is shown in Fig. 4. The exact deformation is displayed by the solid line. The above exercise indicates that acceptable accuracy may only result with high order elements. A 4 x 4 mesh will then contain more DOF than when using 4-noded elements. A true comparison should have been based on an equivalent number of DOF. A very small improvement of 1% with ACMBC was noted, when modelling with a 6 x 6 (245 DOF) uniform element quarter-plate mesh. 4. MODIFYING
ACMBC
Brebbia and Connor reported Ku for the clamped test case. At q = 200 their central displacement differed by 23% from ACMBC, see Table 3. Now the major difference between the two finite element formulations was in the description of the rotation terms in [w]:. The expressions for dw /dx, and dw /dye given by (4) are approximate and so it should be reasonable
to take a fraction of the weighted averages. Figure 5 illustrates modelling characteristics for the simply supported example when the terms are equally scaled by 0.75-0.833. At each load increment the arrowed line shows the change as the scale factor decreased. By reducing the terms in [w]: the contribution due to non-linearity in the tangential stiffness matrix was lessened, thus increasing transverse displacements and stresses. Higher loading produced larger change. The value of &, (Fig. 5a), followed the analytical when the scaling factor was about 0.8. The value for ex,u (Fig. Sb), which was originally within I%, is now too high, since both d,,,,, (Fig. 5c) and r?X,b,D (Fig. 5d) have increased. It can also be noted that the magnitude of change was position sensitive (i.e. compare ~x,pl.Dwith ~X.,l,~). Good modelling was unattainable by just changing the two rotation terms, and in the case of (5,,n, cX,,,,,u and &,Pl,c the error was worse. This modification emphasised the original behaviour where d,,,,n was always in excess, and 6x.6,Dwas always below
Table 3. Percentage errors with non-linear finite elements for the clamped support test case at 4 = 200 Analysis
DGF 4x4
CD
=xX.0
0X.C
Simply polynomial ACMBC (4)
Laminated
125
16.0 L
l.OL
5.0 L
Brebbia and Connor [IO]
Simply polynomial ACMBC (4)
Isotropic
125
1.5 H
Thomas and Gallagher [24]
Inconsistent quad. ‘4 triangles’ (4)
Isotropic
?
4.4 H
Pica, Wood and Hinton [40]
Heterosis QH (9) Serendipity QS (8) Langrangian QL (9)
Isotropic
425 325 405
1.8 L 1.8 L 1.7L
2.0H 3.0 H 3.5H
1.4L 4.0 L 1.3L
Reddy and Chao [42]
Serendipity (4 linear) Serendipity (8 quartic)
125
18.5 H
325
4.0 H
Authors
Element (nodes)
Mottram
Laminated
603
A simple non-linear analysis of multi-layered rectangular plates Simply supported or ctamped B or C
Fig. 4. Transverse displacement along section BD isotropic test cases.
for
the analytical. It was therefore proposed that by changing the magnitude of one of the linear stiffness matrices a correction to this could be made. As the bending contribution had to increase relative to in-plane (in the central region), this could be achieved by reducing the magnitude of the bending stiffness
coefficients. Mottram and Selby [ll] have shown that the bending element predicted accurately linear measured response of multi-layered plates. This posed a restriction on scaling if it was not to lead too poor linear modelling. Scaling bending coefficients by 0.9, which was the maximum reduction permissible, produced a change at 4 = 200 (from the original ACMBC), in both the central displacement and central surface stress of < 1%. This change was insufficient to make this modification to ACMBC useful. Figure 6 shows the effect of increasing in-plane stiffness coefficients by a factor from 1.3 to 1.5. In the cases of Go (Fig. 6a), and dX,D(Fig. 6b), the change was not easily recorded. do was now less than that obtained with basic element. Of relevance from this exercise was the fact that cX,.owas now less than analytical. This occurred because 6,,,,n had substan-
(b) 200
-
si
00
100 -
0
csD
f
tiAL__ IO
20
3’0
)-
Fig. 5. Comparison of non-dimensional values for simply supported, isotropic (v = 0.3) square plate under uniform presure when scaling non-linear stiffness coefficients of element ACMBC from 0.75 to 0.833. (a) Central transverse displacement; (h) central surface tensile stress; (c) in-plane surface tensile stresses; (d) bending surface tensile stresses.
604
J. T. MOITRAM
-
Levy
Fig. 6. Comparison of non-dimensional values for simply supported, isotropic (v = 0.3) square plate under uniform pressure when scaling in-plane stiffness coefficients of element ACMBC from 1.3 to 1.5. (a) Central transverse displacement; (b) central surface tensile stress; (c) in-plane surface tensile stresses; (d) bending surface tensile stresses.
tially decreased (cf. Figs 6c and 2c), while bx,b,Dhad
5. LAMINATE RESULTS
slightly increased (cf. Figs 6d and 2d). Good modelling was also not possible by just scaling in-plane stiffnesses. A successful model was suggested through a combination of reducing the rotation terms and increasing the in-plane stiffnesses. By application of a trial and error procedure this was found when dw/dx, and dw/dy, were scaled by 0.75 and in-plane coefficients by 1.3. Figure 7 illustrates the improvement for the simply supported example. Only stresses bX,P,,c(Fig. 7c), and fX,,4A(Fig. 7d), were not within 1% of analytical. The same was found for the clamped case, but for conciseness has not been reproduced here. Comparing the new ACMBC with other elements in Table 2, it can be seen that it now has similar performance to higher order elements, but at a much lower cost.
TO try to find out if the scaling factors were universal for all types of problems, and especially when the plates consisted of generally orthotropic layers, further test cases were sought. Suitable examples were somewhat scarce, since of the available ones, Chang and Sawamiphakdi [41] studied a single honeycomb sandwich problem and both Reddy et al. [42,43] and Noor and Mathers [44] favoured cross-ply and anti-symmetrical problems. Noor and Mathers did study two symmetrical quasi-isotropic laminates but like Reddy were concerned with shear deformation in thick plates. Recently, Noor [45], using a two-step hybrid analytical technique, solved an ideal problem. The plate was square, all edges were clamped (u = u =
A simple non-linear analysis of multi-layered rectangular plates
605
[ I
I
I
.
Levy
-
l
$
modified
200
-
;i
c
IO
2C
20 s
0
DC
DA
200-
4
5
15
IO
Fig. 7. Comparison of non-dimensional values for simply supported, isotropic (v = 0.3) square plate under uniform pressure using the modified element ACMBC. (a) Central transverse displacement; (b) central surface tensile stress; (c) in-plane surface tensile stresses; (d) bending surface tensile stresses.
w = dw/dx = dw/dy = 0), and the loading was uniform pressure. It consisted of 16 orthotropic layers, with the symmetrical lay-up (45”, -45”, 0”, 0”, -45”, 45”, 90”, go”),. Material properties of each lamina were: E 1 l/E22 = 10.0, G 12/E22 = 0.49, v 12 = 0.38. Shown in Fig. 8 are computed values for central displacement, $, using a 6 x 6 full-plate mesh. The non-dimensional load factor was given by 4 = qS4/E22. A value of S = 120 was used by Noor so his results can be those expected with thin plate assumptions. The modified ACMBC element was found to predict nearly identical values to the hybrid technique and a finite element reported by Chia [46]. The very small discrepancies between the three approaches was very encouraging as they demonstrated that the modified ACMBC could analyse laminates, and in particular the biaxial test. Biaxial tests were performed on a range of laminated rectangular plate specimens [4]. In the test the specimen was supported on four steel ball supports, C.A.S. 2414-E
near to each corner. Patch loading was applied centrally through a square steel loading bar with a rubber foot pad. Measurements of central transverse
zcool
15m-
i
l
Hybrid technique (NOOR 45) Finite element (CHIP, 46) Modified ACMBC
IOCXJ -
500-
a
Fig. 8. Comparison of non-dimensionalised central displacement for clamped square 164ayer laminate plate under uniform pressure.
J. T. Mo-rra~~
606
0
0 Point support
I
.
B
6
:-
-
Fig. 9. Test arrangement and 36-element quarter-plate mesh for a 34-layered specimen.
displacement and strains were taken after each load increment until first fibre failure. Illustrated in Fig. 9 is the test arrangement and a 3delement quarterplate model for a 3Clayered specimen. Parameters for the test and modified stiffnesses[ll] are: lay-up (90”, O”, 90”, 45”, O”, -45”, O”, O”, -45”, O”, 45”, O”, o”, 45”, 0”, -45”, Oo), A = B = 260 mm, A, = B, = 200 mm, A, = BP = 20 mm, h = 4.06 mm; lamina material properties: El 1 = 0.149 x 10” N/m*, E22 = 0.289 x 10” N/m*, G12 = 0.577 x 10” N/m*, v 12 = 0.3, t = 0.1194mm. For this specimen S was 49.3, and so thin plate analysis was applicable. Figure 10 presents non-dimensional numerical values with corresponding measurements. Results from ACMBC were the average of the two quarter-
plate models which compensate for warping [l 11. To solve .five increments of load to Q = 18350, 970 set of CPU were used. The load factor is given by Q = qpAi/E22h4, where q,,= P/A,B, is the pressure over patch area. P is the applied load. After completion of the program the IBM 370/168 had been replaced by an AMDAHL 47O/V8. With this computer nine increments to q = 30550 took 740 sec. Results from the basic element have not been given because the analysis went numerically unstable at 4 = 14,250. Up to q = 22,400 the modified ACMBC can be seen to give a comparable picture to that found for isotropic problems with the initial element, see Figs 2 and 3. At this load the measured central transverse displacement, EC, was underestimated by
I
. -
l
0
,
,
Experiment Modified ACMBC Linear ACMBC
.
. . .
.
. .
.
.o .
.
0
.
. .O .. . . 2 (lo-3)
Fig. 10. Expezimental-computed
4
6
6
IO
12
(YC
comparison with modified element ACMBC of non-linear behaviour for a 34-symmetrical layered square plate.
14
.
A simple non-linear analysis of multi-layered rectangular plates
over lOO%, whilst the error for central surface strains &, and & (and thus stresses), was very small. For loading above 4 = 22,400 the numerical strain in the x global direction grew at a much higher rate than measured, and vice versa in the y direction. Over the same load range the computed displacement increased at a much slower rate than that measured. It is only speculation, but from other tests on a range of plates with various lay-ups and test parameters [4], it is suggested that perhaps the non-linear behaviour of the laminates has contributing factors not allowed for by the assumptions. These additional factors could involve effect of support movement, material non-linearity, position of neutral axis, geometric nonlinearity, resin failure under tension, and preferential bending about one of the global axes [4]. The comparison in Fig. 10 indicates that the effect of these became dominant when the maximum displacement exceeded about three plate thicknesses. 6. CONCLUSIONS
A highly efficient, non-linear, laminate, finite element is developed based on classical thin-plate theory and Von Karman strain expressions. The simple displacement 4-noded rectangular element is referred to as ACMBC, after Adini and Clough, and Melosh who formulated the linear element, and Brebbia and Connor who developed the isotropic non-linear element. A great amount of effort was spent formulating a powerful element capable of solving quarter-plate models consisting of a least 36 elements and up to 20 layers in half the thickness, without excessive need for computing resources. To this end, numerical integration schemes for evaluation of all stiffness coefficients were removed by a new description of element rotations. Computed results are presented for isotropic test cases. Good modelling was unobtainable using the basic element. Central transverse displacement was greatly underestimated. Total stresses (or strains) were excellent, but the distribution between bending and in-plane components was very poor. This was corrected by scaling non-linear and in-plane stiffness coefficients. To demonstrate that the scaling factors were universal a laminate test case was employed. Results using the modified ACMBC element were compared with measurement from bending of a rectangular multi-layered plate. The investigation showed that computed solutions gave a response for laminates up to a central displacement of three plate thicknesses similar to that achieved for isotropic problems with the basic element. For central displacements above this, additional factors are believed to result in larger numerical error. REFEBENCJB 1. J. V. Noyes,
Composites in the construction of the Lear Fan 2100 aircraft. Composites 14, 129-139 (1983).
607
2. E. Fitzer, Carbon based composites. J. Chim. Phys. 81, 717-733 (1984). 3. S. W. Tsai and H. T. Hahn, Introduction to Composite Materials. Technomic, Stamford, CT (1980). 4. J. T. Mottram, The failure of carbon fibre reinforced laminated plates under biaxial stresses. Ph.D. thesis, Department of Engineering Science, Durham University (1984). 5. N. J. Pagano, Exact solution for rectangular bidirectional composites and sandwich plates. J. Comp. Muter. 4, 20-34 (1970). 6. N. J. Pagano, Stress fields in composite laminates. ht. J. Solids Struct. 4, 38%$00 (1978). 7. H. T. Hahn, Non-linear behaviour of laminated composites. J. Comp. Mater. 7, 257-271 (1973). 8. R. M. Jones and H. S. Morgan, Bending and extension of cross-ply laminates with different moduli in tension and compression. Comput. Struct. 11, 181-191 (1980). 9. A. K. Noor and S. J. Hartley, Effect of shear deformation and anisotropy of the non-linear response of composite plates. In Development of Composite Materials (Edited by G. Holister), pp. 55-65. Applied Science, Barking, U.K. (1977). 10. C. Brebbia and J. Connor, Geometric nonlinear finite element analysis. J. Engng Mech. Div., Proc. Am. Sot. Civ. Engrs 95 (EM2), 463483 (1969). 11. J. T. Mottram and A. R. Selby, Bending of thin laminated plates. Comp. Struct. 25, 271-280 (1987). 12. S. Levy, Large deflection theory for rectangular plates. Proc. Symp. Appl. Math AMS, 1, 197 (1969). 13. R. L. Spilker and D. M. Jakobs, Hybrid stress reducedMindlin elements for thin multi-layered plates. Int. J. Namer. Methoa!s Engng 23, 555-578 (1986). 14. K. C. Rockey, H. R. Evans, D. W. Griffiths and D. A. Nethercot, The Finite Element Method. A Basic Introduction to Engineers. Crosby Lockwood Staples (1970). 15. R. J. Melosh. Basic derivation of matrices for the direct stiffness method. J. AIAA 1, 1631-1637 (1963). 16. A. Adini and R. W. Clough, Analysis of plate bending by the finite element method and report. Natl. SC; Foundation USA. G7337 (1961). 17. T. M. Roberts and D. G. Ashwell, The use of tinite element mid-increment stiffness matrices in the postbuckling analysis of imperfect structures. Int. J. Sofiak Struct. 7, 805-823 (1971). 18. G. Kirchoff, Uber das Gliechgewicth und die bewegungeiner elastischen scheile. J. Reine Angew. Marh. 40, 51-88 (1850). 19. Y. C. Fung, Foundation of Solid Mechunics. PrenticeHall, Englewood Cliffs, NJ (1965). 20. J. E. Ashton and J. M. Whitney, Theory of laminated plates. Progress in Material Science, Vol. IV. Technomic, Stamford, CT (1970). 21. 0. C. Zienkiewicz, The Finite Element in Structural Continuum Mechanics. McGraw-Hill, London (1976). 22. T. Kawia and N. Yoshimura, Analysis of large deflection of plates by finite element method. Int. J. Numer. Methods Engng 1, 123-133 (1969). 23. P. G. Beraan and R. W. Clot&. Lame deflection analysis ofl plates and shallow she& us&g the tinite element method. Int. J. Numer. Method Engng 5, 534-556 (1973). 24. R. G. Thomas and R. H. Gallagher, A triangular thin shell finite element: nonlinear analvsis. NASA CR-2483 (1975). 25. C. Prato, A mixed finite element for thin shell analysis. Ph.D. thesis, Department of Civil Engineering, Massachusetts Institute of Technology (1968). 26. R. H. Gallagher and J. Padlog, Discrete element approach to structural stability analysis. J. AIAA 1, 1437-1439 (1963). 27. 0. C. Zienkiewicz, The Finite Element Method in Engineering Science. McGraw-Hill, London (1971).
J. T. MOTTRAM
608
28. J. P. H. Webber, Experimental and theoretical bending deflections of laminated plates. J. Struct. Anal. 7, 87-96 (1972). 29. R. H. Gallagher, R. A. Gellatly, J. Padlog and R. H. Mallett, A discrete element procedure for thin shell instability analysis. J. AIAA 5, 138-145 (1967). 30. S. Mau and R. H. Gallagher, A finite element procedure for nonlinear pre-buckling and initial postbuckling analysis. NASA CR-1936 (1972). 31. J. T. Gden, Numerical formulation of nonlinear elasticity problems. Proc. Am. Sot. Ciu. Engrs 93(ST3), 235-255 (1967).
32. L. D. Hofmeister, G. Greenbaum and D. Evenson, Large strain elastic plastic finite element analysis. J. AZAA 9, 1248-1254 (1971). 33. R. D. Richtmeyer and K. N. Morton, D@rence Methoak for Initial Value Problems. Interscience, New York (1967). 34. R. H. Gallagher, Finite element analysis of geometrically nonlinear problems. In Theory and Practice in Finite Elemenr Structural Analysis (Edited by Yamada), pp. 109-125. SEM 74/32064 (1973). 35. B. E. Greene, Stiffness matrix for bending of a rectangular element with initial membrane stresses. Structural Analysis Research Memo. no. 45, The Boeing Co., Seattle (1962). 36. M. Yoshiki, T. Kawia and N. Yoshimura, Matrix methods of analysis of ship structures, 3rd Report. J. Sot. Nau. Arch. Japan 123, 188-195 (1968). 31. T. Fujino and K. Ohsaka, Static structural analysis of
suspension bridges. Mitsubishi Nippon Heavy-Znd Tech. Rev. 3, 17-23 (1966). 38. R. W. Clough and C. A. Fellipa, A refined quadrilateral element for-analysis of plate-bending. Pro;. &d Conf. Matrix Methods in Structural Mechanics. AFFDL-TR68-150, Wright Patterson AFB, OH, ‘pp. 399-440, (1968). 39. C. T. Wang, Non-linear large deflection boundary value problems of rectangular plates. NACA TN 1425 (1948). 40. A. Pica, R. D. Wood and E. Hinton, Finite element analysis of geometrically non-linear plate bending behaviour using a Mindlin formulation. Compur. Struct. 11, 203-215 (1980). 41. T. Y. Chang and K. Sawamiphakdi, Large deformation analysis of laminated shells by the finite element method. Comput. Srruct. 13, 331-340 (1981). 42. J. N. Reddy and W. C. Chao, Large deflection and large amplitude free vibration of laminated composite material plates. Compur. Srruct. 13, 341-347 (1981). 43. T. Kuppusamy and J. N. Reddy, A three-dimensional nonlinear analysis of cross-ply rectangular composite plates. Comput. Struct. 18, 263-272 (1984). 44. A. K. Noor and M. D. Mathers. Nonlinear finite element analysis of laminated composite shells. Computing Methods in Nonlinear Mechanics. Inst. Comp. Mechs, Texas, pp. l-999 (1974). 45. A. K. Noor, Hybrid analytical technique for nonlinear analysis of structures. J. AZAA 23, 938-946 (1985). 46. C. Y. Chia, Nonlinear analysis of plates, McGraw-Hill, New York (1980).