- Email: [email protected]
Simulation of the stress computation in shells
SIMULATION OF THE STRESS COMPUTATION IN SHELLS M.%LAMAt Jet Propulsion Laboratory, California Institute of Technology. 4800Oak Grove Drive, Pasadena, CA 91103,U.S.A. and S. UTKUS Duke University, Durham, North Carolina. U.S.A.
(Receioed28 Nouaber 1977;‘receicedfor publication 20 Aptif 1978) Abstract-A self-teaching computer program is described, whereby the stresses in thin shells can be computed with good accuracy using the best tit approach. The program is designed for USC in interactive game mode to aiiow the structural engineer to learn about (If the major sources of ~~c~~s and associated errors in the compumtion of stresses in thin she& (2) ~ssib~e ways to reduce the errors, and (3) trade-off between com~tationa~ cost and accuracy. Included are derivation of the compumtion~ approach, program description. and several examples illustrating the programusage.
1. INTRoDucrIO~ Accurate stress computation in the analysis of thin shells by displacement methods requires special attention, even when the displacements have been accurately computed. The stresses are more sensitive than the displacements to numerical i~cc~acies, modeling ~scon~uit~s and approximations of the shell behavior and geometry. Although several authors have dealt with ways to improve upon the computed stresses[ 1-71, it is unfortunate that this particular phase of the discrete shell analysis is least appreciated by practitioners who have limited acquaintance with the subject. One reason for this is that some of the sources of inaccuracies are subtle and can be easiiy overlooked. For example, it is customary to refer the stresses in shells to local o~hogon~ coordinate systems, oriented such that the local third axis coincides with the normal to the middle surface at the node of interest. Slight modeling inaccuracies in the orientation of these coordinates cause gross errors in the computed stresses. Another reason is that solutions found in the literature are often dedicated to dealing with a specific finite element or a specific finite diierence scheme. Although perfectly valid, such so~u~ons do not offer much help to someone seeking a basic unders~ding of the subject. Motivated by these observations, this paper presents a self-teaching computer program by which the stresses in thin shells can be simulated with good accuracy. The program is designed for use in an interactive game mode that allows the user to gain a first hand experience in; (1) the major sources of difficulties and associated errors in the stresses, (2) possible ways to reduce the errors, and 8) the trade off between the compu~tion~ accuracy vs computational cost. The program can also be used in a production mode to compute the stresses at the nodes of a given discrete shell model whose deformations are known. The deformations may have been previously obtained for a finite Member of the Technical Staff. SProfessor of Civil Engineering.
element program such as ELAS, NASTRAN, or SAP, or from a finite difference scheme, or even analytically. With the deformations computed in another program, one may have the stresses computed in this program for checking purposes or to obtain more accurate stresses at certain critical areas in the model of interest. Here, the user is given control over values of several parameters that determine the degree of accuracy of the computed stresses. Z MEl-HODOF %XUTlON A. Preliminaries
To minimize the program size and cost, the following ~s~ptions and i~i~t~ns are imposed; (I) only isotropic elastic thin shells with uniform thickness, at least in the neighborhood of interest, are assumed, (2) no initial strains (due to temperature or otherwise) are allowed, and (3) one shell example is provided for checking the accuracy of the computed quantities. Other assumptions and limitations’ will be discussed as they arise in subsequent sections. Because this paper is concerned only with the computational aspects of the stresses in shells, it is assumed that the shell geometry and deformations ~displa~ments and rotations) are already known with suflicient accuracy at a node set. A node set, a, (i = 0,1, . . . N2) is defined here as the collection of nodal points lying on the shell middle surface. These include the current node, no, at which the stresses are desired, and a set of neighboring nodes, nj, (i=1,2 ,... N2), located in the immediate vicinity of the current node. A right handed Cartesian overall coordinate system xyz (Fig. 1) is used to describe components of the position vector gi, displuc~~t vector &;,and rotation vector &, for any node ni in the node set; Cj+
=
[Xi
_Yi
Zi]
d:=[drx
di,
&‘=[&
Bjy 6:].
diz] (i=O,l,...N2)
(1)
With respect to node no, one may form the relative
580
M. SALAMAand S. UTKU
Both the linear and parabolic fit are special cases of eqn (4). In the linear fit, h,, and hz are the only unknown parameter (E2 = 2). and hS = . . . hs = 0. In the parabolic fit, h,-h5 are unknown (E2 = 5) and h6 = h, = h,, = 0.
Assuming that a representation of the type in eqn (4) is valid, at least in the neighborhood of the node set. eqn
(4) may be evaluated at each neighboring node to obtain N?-condition equations : V’h=f
(5)
in which f’=[f,fi...fN*l v = (VI v2
. .
. VNZ]
vi’ = [a: yf (f:)’ . . . (i’,)‘], - i = 1,. . . N?. Fig. I.
vector of position. c, (Fig. 1). displacements &. and rotations 8:
I%~=[&
4,
(1,J (i=l,2,...N2)
&‘=[&
e;,
e;,]
(2)
where .&=
xi - ~0, . . , &
=
dir - do,, . . . , 8, = &x - @ix, . . ..
The vectors in eqn (2) may also be expressed in a right handed Cartesian local coordinate system, x’y’:‘; whose origin is located at node no, its z’-axis is coincident with the local normal n to the shell’at node no; and its x’-axis is (1) either coincident with the principal direction having the algebraically smaller curvature, or (2) in the plane containing n and the vector joining the current node no to its first neighbor n ,. Thus, eqn (2) may be alternatively written in terms of the local coordinates x’v’z’; PI= [f: &=[d:,
$
Z’] 6:,
e’: = [e’;, $.
In evaluating the state of stress at node no from known information at members of the node set, one is justified in placing different weight (emphasis) on information associated with different nodes. Other than the usual equal weight distribution, a weight assignment scheme is provided in which the weight at node ni is inversely proportional to its distance from’ node no. This scheme has the advantage of preserving the local essential characteristics of the computed stresses and their true variation. Thus, the weighted condition equations are modified by the diagonal weight matrix W so that.
dj:]
(i = I, 2, . . . N2)
(3)
W-VTh= Wf
(6)
W = diag ( Usi), and si = (ZiTZi)“*.
(7)
in which
When N2 = E2. unique solution of h can be found from h = G-‘Wf
(8)
where G = WV? However, when N2> E2, the fit’parameter vector h is found by least squares from
&,I.
h=BGTWf
According to the usual thin shell theory ,& = 0. B. Best fit approach Instead of the discrete description of eqn (3). smooth functional representation may be achieved by fitting either a linear, parabolic, or quadratic continuous function, f, to the Bosition El, each of the displacement components &, d&, a,, and each of the rotation components A, t%,.The number of unknown parameters, E2, in each function determines the minimum number of neighboring nodes N2 that must be used in the computation. For example, in a quadratic fit (E2 = 8);
(9)
where B = (GTG)-‘. Equation (9) is valid whether f represents the shell surface geometry or the shell deformations. In any case. the solution accuracy for the vector of fit parameters h has special significance because it directly affects the accuracy of the computed stresses. Shell surface fitting. The.quadratic form of ur in eqn (4) is appropriately interpreted here by taking f, = I: in eqn (6). For later reference, the vector of fit parameters for this case will be termed ho, and the associated function fo. Notice that elements of ho are related to the direction-cosines of the local normal n at node no by: II= = [Y**Yi*- l]/[(r:,*)* + (fly.)2+ 11”’
vTh=f where VT = [a’, j’, (n’)Z.,‘,‘, (y’)2.f’i’, j’f’. (#I and the vector of fit parameters hr = (h, hz . h,).
(4)
(10)
where & = hoI and 2;. = ho2. Also, the normal and cross curvatures at node no, and f’ and .v’directions are: r;,,, = hoa;
f;l.C. = h,;
?,_:= ho5
(11)
Simulationof the stress computationin shells
581
A. Data input The amount of data inputs required by the prop. largely depends on the mode of operation. Either a production mode or a game mode may be selected. Both modes require the user to specify values of the following control parameters: NO= 0, or 1 to select the mode as either game or production, respectively. Nl = 0, 1, or 2 to allow selection of the interpolation used to best fit the shell surface by either linear, parabolic, or qu~~tic function. Ft = -1, 0. 1, or 2 to allow selection of the interpolation used to best fit the deformations as either lineal, linear, parabolic, or quadratic. N2=The total number of ne~h~~ng nodes in the node set, excluding the current node of interest. A maximum of N2 = 20 is currently imposed. In the game mode, the location of node no, and the distances of the neighboring nodes from no are determined randomly by using a random num~r generator. N3 = 0, or positive integer to be used by the program B = BG’WF (12) as the maximum number of iterations to improve upon values of the matrix of direction cosines L defining the in which each of 111 and F contain five columns b_r*hz,... k, and trt b . . . f5 co~esponding to f = I?:, dY, local coordinate system at the current node. If N3 = 0, d:, B:, and e’;, respectively. independent of the fitting the user should provide L. N4 = 0, or 1 to select the type of weight assignment for exploiting eqn (41, the so called linear fitting [8,91 is also each node of the node set during interpolation. N4 = 0 incorporated in the program for comparison purposes. implies equal weights, while N4= 1 assumes weights Computed strains, curvatures, and stresses. Once the inversely prounion to the distance between the curfit parameter matrix H has been dete~ined from eqn (121, the middle surface strains and curvatures can be rent.node and each neighbor. NS = 0, or 1 indicating whether or not the membrane found[7], stress resultants are desired at the current node. -, N6 = 0, or 1 indicating whether or not the stress couru = d xs’ - hX1 -t ples are desired at the current node. 6.Y” dy,v = hzt N7 = 0, or 1 indicating whether or not the transverse shear stress resultants are desired at the current node. yxy= 6:.,. + &_r = hz, + h,z N8=0, or 1 depending upon whether or not the disyxz = d:.,. c 6:,~ = h,, - &O placements are available at members of the node set, for yyz = 6:,* 4 rTy,r= h23- ($0 (13) stress compu~tion. -, N9 = 0, or 1 depending upon whether or not the rotaXXX = BYS’- h,s tions are available at members of the node set, for stress -2 computation. Xyy= - 0 r.g’= - h24 MO= 1, or 0 indicating whether the computed local & = -h,, + h25. XV = -&+ coordinates at the current node are principal or nonprincipal. If the rotations Ji, iY are not available, the curvatures R = 0, or positive constant >O.Ol indicating whether in eqn (13) are replaced by the randomly generated distances between no and its neighbors I should be used as they are, or scaled by R. Note that R = 0 is treated as if F2 = 1. The value of R (141 simulates a mesh refinement. If the production mode is selected (NO= 11, the following data must also be supplied by the user. AU data inputs must refer to a rectangular overall coordinate The membrane stress resultants h$, NY,N,, the stress system: couples M,, MY, MrY and transverse shears QX, Q, are 1. Elements of a 3 x 3 matrix L, expressing the overall determined from the usual stress resultants-strain relax, y, z components of a vector parallel to each of the tionships~l0~. three o~hogona1 local axes x’, y’, z’; required only if N3=0. 2. Values of the elastic modulus, shear modulus, and The program is written in the basic languaget which is shell thickness in the nei~~rho~ of the current node. available with some variance on most computer in3. Description of the position vectors of all the nodes stallations, especially the so-called mini-computers, in the node set, in the overall xyr coordinate system. 4. Components of the displacement vector at each tThe basii laqpuage is a hiigh-level inUmu%ivecomputer node in the node set, irr the overall xyz coordinate language. It combines English and simple slgebra to $ve in- system. if the requested stress computation should utilize structio~s to the computer. The rcsu!t is shorter programs of such data. simple S~MWCwhich are easy SOwrite and anderstaad. 5. Components of the rotation vector at each node in
from which the principal curvatures and their directions can be corn~t~~?~. In thin shells, many loading conditions result in displacements in the normal directions, at least one order of magnitude greater than the tangential ones. Slight errors in the o~en~tion of n, i.e. hot and ltc~ in eqn (101, can produce signilkant tangential component of the normal displacement. When added to the usually small but true tangential displacements, such error-produced displacement can easily affect the computed membrane forces. Likewise, errors in the o~en~tion of the local normal will affect the accuracy of the computed stress couples and transverse shears. This difficulty is treated here by iteratively improving upon the definition of the local normal (Section 3). Shell deformation fitting. It is performed separately for each of the three displacement components k., 6:, 6:, and for each of the rotation components 8: and ii. Thus, in the general case, eqn (9) is rewritten as
582
M.!~.AMA and S. UTKU
the node set, in the overall xyz coordinate system, if the requested stress computation should utilize such data. In the game mode (NO= 01, the stress simulation is tested by a specific shell example whose explicit solution for displacements and rotations are used by the program to compute automatically the data described in items l-5 above. 3. Output of results Depending upon the values assigned by the user to the thirteen control parameters, some or all of the following quantities may be output by the program: I. An echo of the latest values of the problem parameters and data items described in l-5 above. whether these were program generated or user supplied. 2. Values of principal curvature directions and principal radii at the current node-as numerically computed by the program. 3. The desired components of the numerically computed stress resultants and/or couples at the current node, in the selected local coordinate system. 4. An estimate of the total number of multiplications used in performing the stress simulation, augmented by a numerical overhead simulating the cost of mesh refinement due to F2 (Section 4). 5. If the game mode is selected (NO= 01, the analytical and computed components of the stress resultants and couples are compared and the magnitude of an error function (Section 4) is computed and printed out. In addition to the above output, helpful diagnostic and error messages are printed out as necessary. C. Functions of elements of the program The program is divided to eight subprograms, the iirst of which is called PGI acts as the driving main program. The remaining ones, PG2, PG3, PG4, PCS, PG6, PG7, and PGS are called upon to perform special functions when needed. A brief description is given next. 1. Subprogram PGl. As the main program, it begins with a brief statement of the program objective, assumptions. and method of solutions. If desired, this introduction may be skipped. If the inputs have been prepared previously (i.e. not tutored), the program proceeds with reading appended values of the control parameters, and each of the five items listed in Section 3-A. The data is checked for completeness and consistency with the appropriate messages printed out. When in the production mode. if the user elects to be tutored in the data preparation, PGl will call upon PG2 and PG4 to perform the task. When in the game mode, if the user so desires, PG2 will be called upon to tutor him in inputing the parameters. The control will then be transferred to PG3 to compute internally items 1-5 of the input data for the game example. To print out the data provided so far, PGl will call upon PGS to do the job. Next, PGl will make the final preparation prior to computing the stresses, by computing values of all the position, displacement, and rotation vectors for all the neighboring nodes relative to the current node. If the weights are not uniform, PGI computes their values from eqn (7). Generating the least square equations is done by PC6 while the stress sedation, cost, and error evaluation is performed by PG7 or PGg. Further preparations for the next run are made by PGI upon returning from PG7 or PG8. 2. Subprogram PG2. is called by PGl to tutor the user in preparing inputs of the required control parameters,
one at a time. When N3 = 0, PG2 will additionally tutor the user, if desired, in inputing the required matrix, L. In this case, before returning control to PGl, the subprogram will check the inputs for consistency and validity so that ortbogonality of the local axes is assured. 3. Program PG3. When in the game mode (NO= 01, PGl will call upon this subpro~ to provide a short description of the example problem and its data. The data includes the material properties, shell thickness, and a set of randomly generated position vectors for the current node and its neighbors. For these nodes the program computes the true displacements and rotations according to the control parameters furnished earlier. It also computes and outputs for the current node, the matrix of direction cosines L of the local axes, the matrix of direction cosines Q of the axes in which explicit solutions are expressed according to Ref. [lo], and the appropriate matrices of stress resultants and couples. All are computed from explicit expressions for the randomly generated nodes. The deflections generated here are used in PG6, and the stresses are used in PG7 to compute the errors. 4. Subprogram PGA. This subprogram is called by PG 1 to tutor the user in preparing input data Items 2-5 of Section 3-A. When done, control is returned to PGl. 5. Sub~~g~ PC& The function of the subprogram is to print out, if desired, all the data prior to starting the stress computation. Such data may have been input by the user, or internally generated. 6. Subpnrgram PG6. Much of the computations needed to establish and improve upon the local coordinate definition i, and to evaluate B and G of eqns (9) and (8) are done in this subprogram according to the following iterative steps: (a) As a fust ~prox~ation, L is computed such that 5, x & gives the direction of 1,. (b) Perform the transformation i?= Lr? to obtain results of eqn (3). (c) Using the current local coordinates, compute B, G, 0 and 4, associated with the surface description for the type of fitting indicated by N 1. . (d) The local normal is computed from eqn (lo), and the new matrix L is establish: L = IIt 12 131 where Ii, I*, and I3 are unit vectors along the lirst, second, and third local axis, respectively. 1, = %/l&f, f3= ollnf, and I2= I3x I,. (e) The new and old I3 are compared. If the same, or if the specified iterations N3 have been cons~ed, go to Step f, otherwise go to Step b. (f) Compute the principal directions and their curvatures if MO= I. (g) Recompute B and G for the latest L, with the proper W taken into account. Steps a-g above are similar to ELAS programl91, except in ELAS, two restrictions are made: N3 = 3, and w= 1. 7. Subprogram FG7. This subprogram is called by PGI to compute the stresses (when Fi > -1) and-cost. First. the overall displacements dr, and rotations & are transformed to &, and & according to the latest L. Also, F is established and B is computed from eqn (12). With these, the appropriate strains of eqn (13) and (14) and the
583
Simulationof the stresscomputationin shells associated stress resultants and couples are obtained. Finally, the computational cost, and in the case of NO= 0, the errors in the stresses are also computed and printed out. 8. Subprogrum PG8. This program is caRed by PGl to compute the stresses when Fl = -1. The computation is identical to that of Ref. [9]. For the cost computation, PG7 is called upon to do the job. 4 ElmcATroNAL Ann PRotluclloN FUTuRlls A. Educationalfeatures As an educational tool, the program can be used in a game mode or pr~uctio~ mode. In the game mode, most of the required data such as data items 2-5 (!&&on 3-A) are generated in the program for the specific shell problem of a uniform thickness right circular cylindrical tank filled with liquid[lO]. The axis of the cylinder is coincident with the overall Z-axis, with the gravity vector in the negative Z-direction. The cylinder is clamped at the bottom Z-O, and is free at the top. The weight of the shell itself is ignored. Its thickness, radius, length, and fluid density are nomi~y taken as Loin., 4O.Oin., 100.0in. and 0.49 Iblin’, respectively. If the user desires, he may change these values. The elastic modulus and shear modulus are respectively, 1.0x lO*psi and 0.4x 106psi. The location of the current node no, and the distances of the neighboring nodes ni from no are determined randomly by using a random number generator. The parameter F2 enables the user to scale the distances between nt, and ni to simulate the mesh re~nement for the game. Analytical expressions for the membrane and bending stress resuhants; N,, N2, Nt2r MS, Mzr 62. C?,, and Q are included in the program so that they can be c_ompged with their n~meri+ly c_omputed_counterparts; Ni, N2, N12, h4,, &, h&2. Q1, and Q2. An error function, Q, of the following form is computed and printed for the node of interest;
Q=
CS = 100 * ( l/F2f3.
Here CS represents the overhead associated with the mesh refinement. Aside from the arbitrary multiplier (lOO), CS above is a plausible estimate of the overhead for each node. This is because the number of nodes is inversely proportional to the mesh size squared, and the cost of displacement analysis is proportional to the square of the number of the entire shell model11 11. Using the cost and/or error function above, the program can be used as a self-teaching tool by a p~~titioner having a first time exposure to the subject of numerical analysis of discntized shell structures. However, one need not limit himself to the built-in example above. The production mode can also be effectively used as a selfteaching tool. To begin with, a network of nodes (mesh) representing a current node and its neighbors, all lying on the midsurface of the shell structure should be established. Having dewing the necessary desc~ptions of position, displacements, and rotations at these nodes, one can create an educational game by allowing five of the control parameters to take on a number of possible values, each representing a degree of approximation. For example, one may assume two possible values for Nl, parabolic and quadratic; two possible values for Fl, linear and parabolic; three possible node sets for N2, perhaps differing in their total number or in their proximity to the current node of interest; three different values for the number of iterations N3 to improve upon the local coordinate definition; and two ways to assign the degree of importance to the neighboring nodes, N4 = 0, N4 = 1. Even with such a modest number of variations on the five parameters, as many as seventy two different combinations can be found, each giving rise to a unique cumpufational srrategy that simulates the state of stress
ANa + AN2 + ANtz + AQ, + AQZ+b(AM, + AM2+ A~*2)/r N,+N,+N,,+Q,+(h+6(Ml+Mz+M,t)/t
where AN, = N, - N,, AN2 = N2 - %, . . . , etc. A reasonable indication of the computational cost is the total number of multiplications (C) performed during; surface Wing and Iocal coordinate definition (Cl), principal axes de~ition (C2), ~fo~ation Wing and strain ~ompu~tion (C3), and stress resultant evaluation (C4). In the following estimate, each square rooting is assumed about 10 multiplications; Cl=(N3+2)*[60+9yN2+E2*(2+E2*E2/3) tE2*N2*(2+10*N4+E2/2)1 C2=30*MO*(MO+l) C3=N2*~~*(N8+Nq)~~*Zl] +E2*[2”21+3*N7 +3*(1-N7)*N6*(1-N9)+5*(NS+N6+N7)] C4= 10+3*NS+6*N6+2’N7 C=Cl+C2+C3+C4+CS. All terms above have been previously defined except 21 and C5;
+N7*N8*N9 (16)
(15)
at the current node in question. The objective of a game may then be stated as follows: given the above seventy two computational strategies, a participant would select one that results in the least cost (number of multipli~tions), and/or least error in the c~put~-s~esses. After having equal number of chances to improve upon their results, the participant with the least cost and error accumulated from all trials is declared a winner. Such a game may be played by one participant alone, or several participants against each other. Giving each participant more than one chance to modify or improve upon his initially selected strategy, has the advantage of allowing him to learn about sources of errors in the ~ornpu~ stresses, how much ~provement can be achieved by selecting more suitable set of parameters, and how much accuracy can be reafizcd for a reasonabie computational cost. The teacher may have as important a role in the game as he wishes, but the game can be conducted without his involvement. His participation may include any of the following: 1. Instead of the built-in shell example previously mentioned. he may ask the students to write a small
584
hi. SALAMA
program for another shell example to be used as a subject for the game. 2. As a coach to the participants, he may wish to modify or redefine the rules for the game and/or the scoring system. 3. Grganize his class lectures to revolve around and take advantage of the need to know situation created by questions that arise as the participants begin to interpret the results of various strategies. 4. Use class discussions to develop new ideas for the participants to follow-on, and perhaps program as alternatives to some of the options already available in the program. This would give the participants an opportunity to gain a modest experience in problems associated with the numerical and pro~~ming aspects of the stress simulation in shells. B. Production features In addition to its usage as an educational tool, the program can be used in the production mode for the purpose of evaluating the effect of various variables on the numerical accuracy of the computed stresses in a given shell structure. Even when the numerical accuracy has been established. the program can be used to compute the stresses at criticaI regions of a thin shell structure if the region of interest is not amenable to analytical methods. or if computing the stresses in another program at limited number of nodes means having to compute them everywhere on the structure. Because the stresses are computed here directly from the displacements and rotations at the nodes, with no requirement for element definition or stiiTness matrix input, the question of data transfer between other finite element or finite difference computer programs and this program presents no difBculty. 5. EYtw Three examples are given to demonstrate some of the options for input, output and program usage. The first two examples utilize the game mode in which the built-in shell problem described in Section 4 is used. In the third exampte the pr~uction mode is exploited. All inputs and outputs are listed in the Appendix. Whether using the game or production mode, one may make the required inputs by choosing to be tutored or by appending the data. To be tutored the user starts by typing: GET-PC?i RUN Each command is followed by a carriage return. From there on, one makes input only as requested-following the question mark. A carriage return should follow each completed input. This is shown in Examples 1 and 2. In the appended form, the user prepares the required data prior to starting the game or production run. This is shown by the data set in Tl at the end of the Appendix. The run is started by typing: GET-PG 1 APP-Tl RUN This form is illustrated by the third example.
and s. UTKU
Example 1. The inputs and outputs for this case are given in the Appendix. Several intended input mistakes, their corrections and the program comments are included in the definition of the thirteen control parameters. A printout of the Sinai values of these parameters, and a listing of the program generated random mesh, the associated displacements and rotations are given at the current node and its ten neighbors. The results begin with a rapidly converging definition of the orientation of the local normal, followed by a comparison between the actual (exact) and the numerically computed components of the membrane stress resultants. and couples. The error ratio and computational cost are also given, all for F2 = 1.0, which is labeled as subcase (1 - a). As mentions previously, smaller values of the scalar parameter F2 simulate tiner networks of nodes, and therefore higher compu~tion~ cost. As expected, this atso means better accuracy. This particular aspect is explored in the subsequent subcases (I - 6) and (I- c), for which F2=0.3 and F2 =O.l, respectively. By examining the remarkable accuracy that can be achieved by selecting a fine network of nodes and by observing the associated very high computational costs, one is reminded of the power and limitations of tinite element and finite difference methods. Example 2. Several options of the program, not included in the previous example, are illustrated here. These include extended printout of the program prologue. built-in shell problem description, quadratic fit for the shell surface, parabolic fit for the deformation, and weight assignment inversely proportional to the distance of generated nodes for F2 = I. The associated dispiacements and rotations, local normal orien~tion and the local axesdefini~on are listed. The resulting comparison between the actual and the numerically computed stress resultants are also given, first for F2 = 1.0 in subcase (2 - a), and subsequently by allowing smaller values of F2 = 0.3 and F2 = 0.1 in subcases (2- b) and (2 - cl, respectively. Unlike the previous example where a linear fit for deformations was used. the current parabolic fit for deformations allows the program to compute the transverse shearing forces. Notice that a comp~ison between the results of this example and the previous one is not valid because the location of the current node is not the same. However, a comparison between results for the different values of F2 in this example lead one to the same general conclusions regarding the computational accuracy and cost. A final subcase (2- d) is given in F2 = 0.3 for the same node set except the ELAS-type Lineal fit for deflection is used. With all other data identical to subcase (2 - 6). it is interesting to note that the lineal deflection tit produces comparable stress accuracy as the parabolic fit. Like the tinear fit, however, it is incapable of computing the transverse shears. This may not be regarded as a serious drawback in regions of thin shells where the membrane state is predominant. Example 3. The appended form of data input (previously prepared in Tl), is demonstrated here in conjunction with the production mode. Included is an echo of the selected values of the thirteen control parameters and a listing of the problem data. Notice that values of F2 other than zero or unity are not permissible in a production mode since the deformations cannot be correctly scaled by the same factor as the relative position of the node. The previously computed node locations and associated deformations must always be input. In the
Simulation of the stress computation in shells present example, data generated from the cylinder problem of Examples I and 2 are used to demonstrate the program input and output for the pr~u&tion mode. With such data one can exercise a variety of options for the remaining twelve parameters. The parameters selected for this example are listed in the Appendix along with the numerically computed stress resultants, couples, and shears. Except for the transverse shears which are usually more difhcult to match, the results shown in the Appendix are in good agreement with the exact solution listed in the Table below.
Table I. Exact solution for Example 3
0.0
1467.2
0.0
+I2057
+30.14
0.0
4.2
0.0
Ac~~ow~edge~~rs-me authors wish to thank the staff of the Civil Engineering Department at Duke IJniversity.and the Structures and Dynamics Section at the Jet Propulsion Laboratory for supporting this work. In part, the work presents the results of one phase of research carried out at the Jet Propulsion Laboratory, Caliiomia Institute of Technology under Contract NAS7. 100, sponsored by NASA. It is currently being pursued using funds from the Computer Structural Analysis Fund of Duke University.
585 REFERENCFS
1. A. J. Morris, A deficiency in current finite elements for thin shell application. fnt. J. Sol&& Brrccrures #II), 331-346 (Sept. I973). 2. R. A. Nav and S: Utku. An alternative for the finite element method. broc. Inr. Conf. VuriationalMeth. Engng Vol. I, pp. 3162-3174University of Southhampton, England (Sept. 1972). D. R. Navaratna, Computation of stress resultants in finite element analysis. A.I.A.A. J. 4(II), 2058-2060 (1966). A. K. Noor and V. K. Khandelwal, Improved finite difference variant for the bending analysis of arbitrary cylindrical shells. UNCIV Rep. R-58 (1%9). 1. A. Strickhn and D. R. Navaratna. Improvements on the analysis of shells of revolution by the matrix displacement method. A.I.A.A. 1. Otll). 2869-2072 (1966). S. Utku, Computation of stresses in triangular finite elements. TR-32-948, Jet Propulsion Laboratory, California Institute of Technology (June 1%6). S. Utku, Stiffness matrices for thin triangular elements of non-zero gaussian curvature. A.I.A.A. J. S(9), 1659-1667 (Sept. 1967). 8. M. Salama and S. Utku, Stress computation in displacement methods for two-material elastic media. Camp. Me&. Appf.
Me&. Engng l&325-338 (1977f. 9. S. Utku. ELAS-A general purpose computer program for
the equilibrium problems of linear structures. Vol. 11, Documentation of the Program-Addendum, TR-32-1240, Jet Propulsion Laboratory. California Institute of Technology (Oct. 1969). IO. S. Timoshenko and S. Woinowsky-Krieger. Theory of Notes nnd Shells. 2nd Edn. McGraw-Hill, New York (1959). Il. C. H. Norris, J. B. Wilbur and S. Utku, Ekmentaly Smcrural Anofysis. 3rd Edn. McGraw-Hill.New York (19761.
586
M. SNMA and S. UTKU
.
..~-....*‘.........*.....~..*...........-.........................*........
............................................................................ ............................................................................
Simulation of the stress computation in shells
587
..rS:eEf.C~ ..:OC:?L.P? ..:C.L:t*::
..-::c:c.:? ...POti‘CC.C. ..COL;oL.PC ..rrc,>‘O‘.“, . ..C.:Cf.,O ..D‘:2cc.‘: +.,o::Cc.l;
(Receioed28 Nouaber 1977;‘receicedfor publication 20 Aptif 1978) Abstract-A self-teaching computer program is described, whereby the stresses in thin shells can be computed with good accuracy using the best tit approach. The program is designed for USC in interactive game mode to aiiow the structural engineer to learn about (If the major sources of ~~c~~s and associated errors in the compumtion of stresses in thin she& (2) ~ssib~e ways to reduce the errors, and (3) trade-off between com~tationa~ cost and accuracy. Included are derivation of the compumtion~ approach, program description. and several examples illustrating the programusage.
1. INTRoDucrIO~ Accurate stress computation in the analysis of thin shells by displacement methods requires special attention, even when the displacements have been accurately computed. The stresses are more sensitive than the displacements to numerical i~cc~acies, modeling ~scon~uit~s and approximations of the shell behavior and geometry. Although several authors have dealt with ways to improve upon the computed stresses[ 1-71, it is unfortunate that this particular phase of the discrete shell analysis is least appreciated by practitioners who have limited acquaintance with the subject. One reason for this is that some of the sources of inaccuracies are subtle and can be easiiy overlooked. For example, it is customary to refer the stresses in shells to local o~hogon~ coordinate systems, oriented such that the local third axis coincides with the normal to the middle surface at the node of interest. Slight modeling inaccuracies in the orientation of these coordinates cause gross errors in the computed stresses. Another reason is that solutions found in the literature are often dedicated to dealing with a specific finite element or a specific finite diierence scheme. Although perfectly valid, such so~u~ons do not offer much help to someone seeking a basic unders~ding of the subject. Motivated by these observations, this paper presents a self-teaching computer program by which the stresses in thin shells can be simulated with good accuracy. The program is designed for use in an interactive game mode that allows the user to gain a first hand experience in; (1) the major sources of difficulties and associated errors in the stresses, (2) possible ways to reduce the errors, and 8) the trade off between the compu~tion~ accuracy vs computational cost. The program can also be used in a production mode to compute the stresses at the nodes of a given discrete shell model whose deformations are known. The deformations may have been previously obtained for a finite Member of the Technical Staff. SProfessor of Civil Engineering.
element program such as ELAS, NASTRAN, or SAP, or from a finite difference scheme, or even analytically. With the deformations computed in another program, one may have the stresses computed in this program for checking purposes or to obtain more accurate stresses at certain critical areas in the model of interest. Here, the user is given control over values of several parameters that determine the degree of accuracy of the computed stresses. Z MEl-HODOF %XUTlON A. Preliminaries
To minimize the program size and cost, the following ~s~ptions and i~i~t~ns are imposed; (I) only isotropic elastic thin shells with uniform thickness, at least in the neighborhood of interest, are assumed, (2) no initial strains (due to temperature or otherwise) are allowed, and (3) one shell example is provided for checking the accuracy of the computed quantities. Other assumptions and limitations’ will be discussed as they arise in subsequent sections. Because this paper is concerned only with the computational aspects of the stresses in shells, it is assumed that the shell geometry and deformations ~displa~ments and rotations) are already known with suflicient accuracy at a node set. A node set, a, (i = 0,1, . . . N2) is defined here as the collection of nodal points lying on the shell middle surface. These include the current node, no, at which the stresses are desired, and a set of neighboring nodes, nj, (i=1,2 ,... N2), located in the immediate vicinity of the current node. A right handed Cartesian overall coordinate system xyz (Fig. 1) is used to describe components of the position vector gi, displuc~~t vector &;,and rotation vector &, for any node ni in the node set; Cj+
=
[Xi
_Yi
Zi]
d:=[drx
di,
&‘=[&
Bjy 6:].
diz] (i=O,l,...N2)
(1)
With respect to node no, one may form the relative
580
M. SALAMAand S. UTKU
Both the linear and parabolic fit are special cases of eqn (4). In the linear fit, h,, and hz are the only unknown parameter (E2 = 2). and hS = . . . hs = 0. In the parabolic fit, h,-h5 are unknown (E2 = 5) and h6 = h, = h,, = 0.
Assuming that a representation of the type in eqn (4) is valid, at least in the neighborhood of the node set. eqn
(4) may be evaluated at each neighboring node to obtain N?-condition equations : V’h=f
(5)
in which f’=[f,fi...fN*l v = (VI v2
. .
. VNZ]
vi’ = [a: yf (f:)’ . . . (i’,)‘], - i = 1,. . . N?. Fig. I.
vector of position. c, (Fig. 1). displacements &. and rotations 8:
I%~=[&
4,
(1,J (i=l,2,...N2)
&‘=[&
e;,
e;,]
(2)
where .&=
xi - ~0, . . , &
=
dir - do,, . . . , 8, = &x - @ix, . . ..
The vectors in eqn (2) may also be expressed in a right handed Cartesian local coordinate system, x’y’:‘; whose origin is located at node no, its z’-axis is coincident with the local normal n to the shell’at node no; and its x’-axis is (1) either coincident with the principal direction having the algebraically smaller curvature, or (2) in the plane containing n and the vector joining the current node no to its first neighbor n ,. Thus, eqn (2) may be alternatively written in terms of the local coordinates x’v’z’; PI= [f: &=[d:,
$
Z’] 6:,
e’: = [e’;, $.
In evaluating the state of stress at node no from known information at members of the node set, one is justified in placing different weight (emphasis) on information associated with different nodes. Other than the usual equal weight distribution, a weight assignment scheme is provided in which the weight at node ni is inversely proportional to its distance from’ node no. This scheme has the advantage of preserving the local essential characteristics of the computed stresses and their true variation. Thus, the weighted condition equations are modified by the diagonal weight matrix W so that.
dj:]
(i = I, 2, . . . N2)
(3)
W-VTh= Wf
(6)
W = diag ( Usi), and si = (ZiTZi)“*.
(7)
in which
When N2 = E2. unique solution of h can be found from h = G-‘Wf
(8)
where G = WV? However, when N2> E2, the fit’parameter vector h is found by least squares from
&,I.
h=BGTWf
According to the usual thin shell theory ,& = 0. B. Best fit approach Instead of the discrete description of eqn (3). smooth functional representation may be achieved by fitting either a linear, parabolic, or quadratic continuous function, f, to the Bosition El, each of the displacement components &, d&, a,, and each of the rotation components A, t%,.The number of unknown parameters, E2, in each function determines the minimum number of neighboring nodes N2 that must be used in the computation. For example, in a quadratic fit (E2 = 8);
(9)
where B = (GTG)-‘. Equation (9) is valid whether f represents the shell surface geometry or the shell deformations. In any case. the solution accuracy for the vector of fit parameters h has special significance because it directly affects the accuracy of the computed stresses. Shell surface fitting. The.quadratic form of ur in eqn (4) is appropriately interpreted here by taking f, = I: in eqn (6). For later reference, the vector of fit parameters for this case will be termed ho, and the associated function fo. Notice that elements of ho are related to the direction-cosines of the local normal n at node no by: II= = [Y**Yi*- l]/[(r:,*)* + (fly.)2+ 11”’
vTh=f where VT = [a’, j’, (n’)Z.,‘,‘, (y’)2.f’i’, j’f’. (#I and the vector of fit parameters hr = (h, hz . h,).
(4)
(10)
where & = hoI and 2;. = ho2. Also, the normal and cross curvatures at node no, and f’ and .v’directions are: r;,,, = hoa;
f;l.C. = h,;
?,_:= ho5
(11)
Simulationof the stress computationin shells
581
A. Data input The amount of data inputs required by the prop. largely depends on the mode of operation. Either a production mode or a game mode may be selected. Both modes require the user to specify values of the following control parameters: NO= 0, or 1 to select the mode as either game or production, respectively. Nl = 0, 1, or 2 to allow selection of the interpolation used to best fit the shell surface by either linear, parabolic, or qu~~tic function. Ft = -1, 0. 1, or 2 to allow selection of the interpolation used to best fit the deformations as either lineal, linear, parabolic, or quadratic. N2=The total number of ne~h~~ng nodes in the node set, excluding the current node of interest. A maximum of N2 = 20 is currently imposed. In the game mode, the location of node no, and the distances of the neighboring nodes from no are determined randomly by using a random num~r generator. N3 = 0, or positive integer to be used by the program B = BG’WF (12) as the maximum number of iterations to improve upon values of the matrix of direction cosines L defining the in which each of 111 and F contain five columns b_r*hz,... k, and trt b . . . f5 co~esponding to f = I?:, dY, local coordinate system at the current node. If N3 = 0, d:, B:, and e’;, respectively. independent of the fitting the user should provide L. N4 = 0, or 1 to select the type of weight assignment for exploiting eqn (41, the so called linear fitting [8,91 is also each node of the node set during interpolation. N4 = 0 incorporated in the program for comparison purposes. implies equal weights, while N4= 1 assumes weights Computed strains, curvatures, and stresses. Once the inversely prounion to the distance between the curfit parameter matrix H has been dete~ined from eqn (121, the middle surface strains and curvatures can be rent.node and each neighbor. NS = 0, or 1 indicating whether or not the membrane found[7], stress resultants are desired at the current node. -, N6 = 0, or 1 indicating whether or not the stress couru = d xs’ - hX1 -t ples are desired at the current node. 6.Y” dy,v = hzt N7 = 0, or 1 indicating whether or not the transverse shear stress resultants are desired at the current node. yxy= 6:.,. + &_r = hz, + h,z N8=0, or 1 depending upon whether or not the disyxz = d:.,. c 6:,~ = h,, - &O placements are available at members of the node set, for yyz = 6:,* 4 rTy,r= h23- ($0 (13) stress compu~tion. -, N9 = 0, or 1 depending upon whether or not the rotaXXX = BYS’- h,s tions are available at members of the node set, for stress -2 computation. Xyy= - 0 r.g’= - h24 MO= 1, or 0 indicating whether the computed local & = -h,, + h25. XV = -&+ coordinates at the current node are principal or nonprincipal. If the rotations Ji, iY are not available, the curvatures R = 0, or positive constant >O.Ol indicating whether in eqn (13) are replaced by the randomly generated distances between no and its neighbors I should be used as they are, or scaled by R. Note that R = 0 is treated as if F2 = 1. The value of R (141 simulates a mesh refinement. If the production mode is selected (NO= 11, the following data must also be supplied by the user. AU data inputs must refer to a rectangular overall coordinate The membrane stress resultants h$, NY,N,, the stress system: couples M,, MY, MrY and transverse shears QX, Q, are 1. Elements of a 3 x 3 matrix L, expressing the overall determined from the usual stress resultants-strain relax, y, z components of a vector parallel to each of the tionships~l0~. three o~hogona1 local axes x’, y’, z’; required only if N3=0. 2. Values of the elastic modulus, shear modulus, and The program is written in the basic languaget which is shell thickness in the nei~~rho~ of the current node. available with some variance on most computer in3. Description of the position vectors of all the nodes stallations, especially the so-called mini-computers, in the node set, in the overall xyr coordinate system. 4. Components of the displacement vector at each tThe basii laqpuage is a hiigh-level inUmu%ivecomputer node in the node set, irr the overall xyz coordinate language. It combines English and simple slgebra to $ve in- system. if the requested stress computation should utilize structio~s to the computer. The rcsu!t is shorter programs of such data. simple S~MWCwhich are easy SOwrite and anderstaad. 5. Components of the rotation vector at each node in
from which the principal curvatures and their directions can be corn~t~~?~. In thin shells, many loading conditions result in displacements in the normal directions, at least one order of magnitude greater than the tangential ones. Slight errors in the o~en~tion of n, i.e. hot and ltc~ in eqn (101, can produce signilkant tangential component of the normal displacement. When added to the usually small but true tangential displacements, such error-produced displacement can easily affect the computed membrane forces. Likewise, errors in the o~en~tion of the local normal will affect the accuracy of the computed stress couples and transverse shears. This difficulty is treated here by iteratively improving upon the definition of the local normal (Section 3). Shell deformation fitting. It is performed separately for each of the three displacement components k., 6:, 6:, and for each of the rotation components 8: and ii. Thus, in the general case, eqn (9) is rewritten as
582
M.!~.AMA and S. UTKU
the node set, in the overall xyz coordinate system, if the requested stress computation should utilize such data. In the game mode (NO= 01, the stress simulation is tested by a specific shell example whose explicit solution for displacements and rotations are used by the program to compute automatically the data described in items l-5 above. 3. Output of results Depending upon the values assigned by the user to the thirteen control parameters, some or all of the following quantities may be output by the program: I. An echo of the latest values of the problem parameters and data items described in l-5 above. whether these were program generated or user supplied. 2. Values of principal curvature directions and principal radii at the current node-as numerically computed by the program. 3. The desired components of the numerically computed stress resultants and/or couples at the current node, in the selected local coordinate system. 4. An estimate of the total number of multiplications used in performing the stress simulation, augmented by a numerical overhead simulating the cost of mesh refinement due to F2 (Section 4). 5. If the game mode is selected (NO= 01, the analytical and computed components of the stress resultants and couples are compared and the magnitude of an error function (Section 4) is computed and printed out. In addition to the above output, helpful diagnostic and error messages are printed out as necessary. C. Functions of elements of the program The program is divided to eight subprograms, the iirst of which is called PGI acts as the driving main program. The remaining ones, PG2, PG3, PG4, PCS, PG6, PG7, and PGS are called upon to perform special functions when needed. A brief description is given next. 1. Subprogram PGl. As the main program, it begins with a brief statement of the program objective, assumptions. and method of solutions. If desired, this introduction may be skipped. If the inputs have been prepared previously (i.e. not tutored), the program proceeds with reading appended values of the control parameters, and each of the five items listed in Section 3-A. The data is checked for completeness and consistency with the appropriate messages printed out. When in the production mode. if the user elects to be tutored in the data preparation, PGl will call upon PG2 and PG4 to perform the task. When in the game mode, if the user so desires, PG2 will be called upon to tutor him in inputing the parameters. The control will then be transferred to PG3 to compute internally items 1-5 of the input data for the game example. To print out the data provided so far, PGl will call upon PGS to do the job. Next, PGl will make the final preparation prior to computing the stresses, by computing values of all the position, displacement, and rotation vectors for all the neighboring nodes relative to the current node. If the weights are not uniform, PGI computes their values from eqn (7). Generating the least square equations is done by PC6 while the stress sedation, cost, and error evaluation is performed by PG7 or PGg. Further preparations for the next run are made by PGI upon returning from PG7 or PG8. 2. Subprogram PG2. is called by PGl to tutor the user in preparing inputs of the required control parameters,
one at a time. When N3 = 0, PG2 will additionally tutor the user, if desired, in inputing the required matrix, L. In this case, before returning control to PGl, the subprogram will check the inputs for consistency and validity so that ortbogonality of the local axes is assured. 3. Program PG3. When in the game mode (NO= 01, PGl will call upon this subpro~ to provide a short description of the example problem and its data. The data includes the material properties, shell thickness, and a set of randomly generated position vectors for the current node and its neighbors. For these nodes the program computes the true displacements and rotations according to the control parameters furnished earlier. It also computes and outputs for the current node, the matrix of direction cosines L of the local axes, the matrix of direction cosines Q of the axes in which explicit solutions are expressed according to Ref. [lo], and the appropriate matrices of stress resultants and couples. All are computed from explicit expressions for the randomly generated nodes. The deflections generated here are used in PG6, and the stresses are used in PG7 to compute the errors. 4. Subprogram PGA. This subprogram is called by PG 1 to tutor the user in preparing input data Items 2-5 of Section 3-A. When done, control is returned to PGl. 5. Sub~~g~ PC& The function of the subprogram is to print out, if desired, all the data prior to starting the stress computation. Such data may have been input by the user, or internally generated. 6. Subpnrgram PG6. Much of the computations needed to establish and improve upon the local coordinate definition i, and to evaluate B and G of eqns (9) and (8) are done in this subprogram according to the following iterative steps: (a) As a fust ~prox~ation, L is computed such that 5, x & gives the direction of 1,. (b) Perform the transformation i?= Lr? to obtain results of eqn (3). (c) Using the current local coordinates, compute B, G, 0 and 4, associated with the surface description for the type of fitting indicated by N 1. . (d) The local normal is computed from eqn (lo), and the new matrix L is establish: L = IIt 12 131 where Ii, I*, and I3 are unit vectors along the lirst, second, and third local axis, respectively. 1, = %/l&f, f3= ollnf, and I2= I3x I,. (e) The new and old I3 are compared. If the same, or if the specified iterations N3 have been cons~ed, go to Step f, otherwise go to Step b. (f) Compute the principal directions and their curvatures if MO= I. (g) Recompute B and G for the latest L, with the proper W taken into account. Steps a-g above are similar to ELAS programl91, except in ELAS, two restrictions are made: N3 = 3, and w= 1. 7. Subprogram FG7. This subprogram is called by PGI to compute the stresses (when Fi > -1) and-cost. First. the overall displacements dr, and rotations & are transformed to &, and & according to the latest L. Also, F is established and B is computed from eqn (12). With these, the appropriate strains of eqn (13) and (14) and the
583
Simulationof the stresscomputationin shells associated stress resultants and couples are obtained. Finally, the computational cost, and in the case of NO= 0, the errors in the stresses are also computed and printed out. 8. Subprogrum PG8. This program is caRed by PGl to compute the stresses when Fl = -1. The computation is identical to that of Ref. [9]. For the cost computation, PG7 is called upon to do the job. 4 ElmcATroNAL Ann PRotluclloN FUTuRlls A. Educationalfeatures As an educational tool, the program can be used in a game mode or pr~uctio~ mode. In the game mode, most of the required data such as data items 2-5 (!&&on 3-A) are generated in the program for the specific shell problem of a uniform thickness right circular cylindrical tank filled with liquid[lO]. The axis of the cylinder is coincident with the overall Z-axis, with the gravity vector in the negative Z-direction. The cylinder is clamped at the bottom Z-O, and is free at the top. The weight of the shell itself is ignored. Its thickness, radius, length, and fluid density are nomi~y taken as Loin., 4O.Oin., 100.0in. and 0.49 Iblin’, respectively. If the user desires, he may change these values. The elastic modulus and shear modulus are respectively, 1.0x lO*psi and 0.4x 106psi. The location of the current node no, and the distances of the neighboring nodes ni from no are determined randomly by using a random number generator. The parameter F2 enables the user to scale the distances between nt, and ni to simulate the mesh re~nement for the game. Analytical expressions for the membrane and bending stress resuhants; N,, N2, Nt2r MS, Mzr 62. C?,, and Q are included in the program so that they can be c_ompged with their n~meri+ly c_omputed_counterparts; Ni, N2, N12, h4,, &, h&2. Q1, and Q2. An error function, Q, of the following form is computed and printed for the node of interest;
Q=
CS = 100 * ( l/F2f3.
Here CS represents the overhead associated with the mesh refinement. Aside from the arbitrary multiplier (lOO), CS above is a plausible estimate of the overhead for each node. This is because the number of nodes is inversely proportional to the mesh size squared, and the cost of displacement analysis is proportional to the square of the number of the entire shell model11 11. Using the cost and/or error function above, the program can be used as a self-teaching tool by a p~~titioner having a first time exposure to the subject of numerical analysis of discntized shell structures. However, one need not limit himself to the built-in example above. The production mode can also be effectively used as a selfteaching tool. To begin with, a network of nodes (mesh) representing a current node and its neighbors, all lying on the midsurface of the shell structure should be established. Having dewing the necessary desc~ptions of position, displacements, and rotations at these nodes, one can create an educational game by allowing five of the control parameters to take on a number of possible values, each representing a degree of approximation. For example, one may assume two possible values for Nl, parabolic and quadratic; two possible values for Fl, linear and parabolic; three possible node sets for N2, perhaps differing in their total number or in their proximity to the current node of interest; three different values for the number of iterations N3 to improve upon the local coordinate definition; and two ways to assign the degree of importance to the neighboring nodes, N4 = 0, N4 = 1. Even with such a modest number of variations on the five parameters, as many as seventy two different combinations can be found, each giving rise to a unique cumpufational srrategy that simulates the state of stress
ANa + AN2 + ANtz + AQ, + AQZ+b(AM, + AM2+ A~*2)/r N,+N,+N,,+Q,+(h+6(Ml+Mz+M,t)/t
where AN, = N, - N,, AN2 = N2 - %, . . . , etc. A reasonable indication of the computational cost is the total number of multiplications (C) performed during; surface Wing and Iocal coordinate definition (Cl), principal axes de~ition (C2), ~fo~ation Wing and strain ~ompu~tion (C3), and stress resultant evaluation (C4). In the following estimate, each square rooting is assumed about 10 multiplications; Cl=(N3+2)*[60+9yN2+E2*(2+E2*E2/3) tE2*N2*(2+10*N4+E2/2)1 C2=30*MO*(MO+l) C3=N2*~~*(N8+Nq)~~*Zl] +E2*[2”21+3*N7 +3*(1-N7)*N6*(1-N9)+5*(NS+N6+N7)] C4= 10+3*NS+6*N6+2’N7 C=Cl+C2+C3+C4+CS. All terms above have been previously defined except 21 and C5;
+N7*N8*N9 (16)
(15)
at the current node in question. The objective of a game may then be stated as follows: given the above seventy two computational strategies, a participant would select one that results in the least cost (number of multipli~tions), and/or least error in the c~put~-s~esses. After having equal number of chances to improve upon their results, the participant with the least cost and error accumulated from all trials is declared a winner. Such a game may be played by one participant alone, or several participants against each other. Giving each participant more than one chance to modify or improve upon his initially selected strategy, has the advantage of allowing him to learn about sources of errors in the ~ornpu~ stresses, how much ~provement can be achieved by selecting more suitable set of parameters, and how much accuracy can be reafizcd for a reasonabie computational cost. The teacher may have as important a role in the game as he wishes, but the game can be conducted without his involvement. His participation may include any of the following: 1. Instead of the built-in shell example previously mentioned. he may ask the students to write a small
584
hi. SALAMA
program for another shell example to be used as a subject for the game. 2. As a coach to the participants, he may wish to modify or redefine the rules for the game and/or the scoring system. 3. Grganize his class lectures to revolve around and take advantage of the need to know situation created by questions that arise as the participants begin to interpret the results of various strategies. 4. Use class discussions to develop new ideas for the participants to follow-on, and perhaps program as alternatives to some of the options already available in the program. This would give the participants an opportunity to gain a modest experience in problems associated with the numerical and pro~~ming aspects of the stress simulation in shells. B. Production features In addition to its usage as an educational tool, the program can be used in the production mode for the purpose of evaluating the effect of various variables on the numerical accuracy of the computed stresses in a given shell structure. Even when the numerical accuracy has been established. the program can be used to compute the stresses at criticaI regions of a thin shell structure if the region of interest is not amenable to analytical methods. or if computing the stresses in another program at limited number of nodes means having to compute them everywhere on the structure. Because the stresses are computed here directly from the displacements and rotations at the nodes, with no requirement for element definition or stiiTness matrix input, the question of data transfer between other finite element or finite difference computer programs and this program presents no difBculty. 5. EYtw Three examples are given to demonstrate some of the options for input, output and program usage. The first two examples utilize the game mode in which the built-in shell problem described in Section 4 is used. In the third exampte the pr~uction mode is exploited. All inputs and outputs are listed in the Appendix. Whether using the game or production mode, one may make the required inputs by choosing to be tutored or by appending the data. To be tutored the user starts by typing: GET-PC?i RUN Each command is followed by a carriage return. From there on, one makes input only as requested-following the question mark. A carriage return should follow each completed input. This is shown in Examples 1 and 2. In the appended form, the user prepares the required data prior to starting the game or production run. This is shown by the data set in Tl at the end of the Appendix. The run is started by typing: GET-PG 1 APP-Tl RUN This form is illustrated by the third example.
and s. UTKU
Example 1. The inputs and outputs for this case are given in the Appendix. Several intended input mistakes, their corrections and the program comments are included in the definition of the thirteen control parameters. A printout of the Sinai values of these parameters, and a listing of the program generated random mesh, the associated displacements and rotations are given at the current node and its ten neighbors. The results begin with a rapidly converging definition of the orientation of the local normal, followed by a comparison between the actual (exact) and the numerically computed components of the membrane stress resultants. and couples. The error ratio and computational cost are also given, all for F2 = 1.0, which is labeled as subcase (1 - a). As mentions previously, smaller values of the scalar parameter F2 simulate tiner networks of nodes, and therefore higher compu~tion~ cost. As expected, this atso means better accuracy. This particular aspect is explored in the subsequent subcases (I - 6) and (I- c), for which F2=0.3 and F2 =O.l, respectively. By examining the remarkable accuracy that can be achieved by selecting a fine network of nodes and by observing the associated very high computational costs, one is reminded of the power and limitations of tinite element and finite difference methods. Example 2. Several options of the program, not included in the previous example, are illustrated here. These include extended printout of the program prologue. built-in shell problem description, quadratic fit for the shell surface, parabolic fit for the deformation, and weight assignment inversely proportional to the distance of generated nodes for F2 = I. The associated dispiacements and rotations, local normal orien~tion and the local axesdefini~on are listed. The resulting comparison between the actual and the numerically computed stress resultants are also given, first for F2 = 1.0 in subcase (2 - a), and subsequently by allowing smaller values of F2 = 0.3 and F2 = 0.1 in subcases (2- b) and (2 - cl, respectively. Unlike the previous example where a linear fit for deformations was used. the current parabolic fit for deformations allows the program to compute the transverse shearing forces. Notice that a comp~ison between the results of this example and the previous one is not valid because the location of the current node is not the same. However, a comparison between results for the different values of F2 in this example lead one to the same general conclusions regarding the computational accuracy and cost. A final subcase (2- d) is given in F2 = 0.3 for the same node set except the ELAS-type Lineal fit for deflection is used. With all other data identical to subcase (2 - 6). it is interesting to note that the lineal deflection tit produces comparable stress accuracy as the parabolic fit. Like the tinear fit, however, it is incapable of computing the transverse shears. This may not be regarded as a serious drawback in regions of thin shells where the membrane state is predominant. Example 3. The appended form of data input (previously prepared in Tl), is demonstrated here in conjunction with the production mode. Included is an echo of the selected values of the thirteen control parameters and a listing of the problem data. Notice that values of F2 other than zero or unity are not permissible in a production mode since the deformations cannot be correctly scaled by the same factor as the relative position of the node. The previously computed node locations and associated deformations must always be input. In the
Simulation of the stress computation in shells present example, data generated from the cylinder problem of Examples I and 2 are used to demonstrate the program input and output for the pr~u&tion mode. With such data one can exercise a variety of options for the remaining twelve parameters. The parameters selected for this example are listed in the Appendix along with the numerically computed stress resultants, couples, and shears. Except for the transverse shears which are usually more difhcult to match, the results shown in the Appendix are in good agreement with the exact solution listed in the Table below.
Table I. Exact solution for Example 3
0.0
1467.2
0.0
+I2057
+30.14
0.0
4.2
0.0
Ac~~ow~edge~~rs-me authors wish to thank the staff of the Civil Engineering Department at Duke IJniversity.and the Structures and Dynamics Section at the Jet Propulsion Laboratory for supporting this work. In part, the work presents the results of one phase of research carried out at the Jet Propulsion Laboratory, Caliiomia Institute of Technology under Contract NAS7. 100, sponsored by NASA. It is currently being pursued using funds from the Computer Structural Analysis Fund of Duke University.
585 REFERENCFS
1. A. J. Morris, A deficiency in current finite elements for thin shell application. fnt. J. Sol&& Brrccrures #II), 331-346 (Sept. I973). 2. R. A. Nav and S: Utku. An alternative for the finite element method. broc. Inr. Conf. VuriationalMeth. Engng Vol. I, pp. 3162-3174University of Southhampton, England (Sept. 1972). D. R. Navaratna, Computation of stress resultants in finite element analysis. A.I.A.A. J. 4(II), 2058-2060 (1966). A. K. Noor and V. K. Khandelwal, Improved finite difference variant for the bending analysis of arbitrary cylindrical shells. UNCIV Rep. R-58 (1%9). 1. A. Strickhn and D. R. Navaratna. Improvements on the analysis of shells of revolution by the matrix displacement method. A.I.A.A. 1. Otll). 2869-2072 (1966). S. Utku, Computation of stresses in triangular finite elements. TR-32-948, Jet Propulsion Laboratory, California Institute of Technology (June 1%6). S. Utku, Stiffness matrices for thin triangular elements of non-zero gaussian curvature. A.I.A.A. J. S(9), 1659-1667 (Sept. 1967). 8. M. Salama and S. Utku, Stress computation in displacement methods for two-material elastic media. Camp. Me&. Appf.
Me&. Engng l&325-338 (1977f. 9. S. Utku. ELAS-A general purpose computer program for
the equilibrium problems of linear structures. Vol. 11, Documentation of the Program-Addendum, TR-32-1240, Jet Propulsion Laboratory. California Institute of Technology (Oct. 1969). IO. S. Timoshenko and S. Woinowsky-Krieger. Theory of Notes nnd Shells. 2nd Edn. McGraw-Hill, New York (1959). Il. C. H. Norris, J. B. Wilbur and S. Utku, Ekmentaly Smcrural Anofysis. 3rd Edn. McGraw-Hill.New York (19761.
586
M. SNMA and S. UTKU
.
..~-....*‘.........*.....~..*...........-.........................*........
............................................................................ ............................................................................
Simulation of the stress computation in shells
587
..rS:eEf.C~ ..:OC:?L.P? ..:C.L:t*::
..-::c:c.:? ...POti‘CC.C. ..COL;oL.PC ..rrc,>‘O‘.“, . ..C.:Cf.,O ..D‘:2cc.‘: +.,o::Cc.l;