Stiffness version of finite element method used for numerical analysis of small elementary regions of solid propellant grain

Key, S. W.

English

NASA Technical Reports Server

MSC R e p o r t  N o .  66-21-1 
FLEXIBLE CASE ANALYSIS 
FOR COMPRESSIBLE 
SOLID PROPEUANT GRAIN MOTORS 
GPO PRICE s 
CFSTI PRICE(S) $ 
S t a t u s  R e p o r t  
C o n t r a c t  N o .  NAS 7-242 
A p r i l ,  1966 
Hard czpy (HC) /e / 
Mfcrofiche 1MF) ,m 
ff 653 July 65 
- 
to  
NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 
P r e p a r e d  by 
S a m u e l  W .  Key 
MATHEMATICAL SCIENCES CORPORATION 
S e a t t l e  Office 
1 1 0 7  N o r t h e a s t  45th S t r e e t  
S e a t t l e ,  Washington 
https://ntrs.nasa.gov/search.jsp?R=19660021325 2020-03-16T20:47:52+00:00Z
TABLE OF CONTENTS 
1 . INTRODUCTION . . . . . . . . . . . . . . . . . . . . .  1 
2 . STIFFNESS METHOD OUTLINE . . . . . . . . . . . . . . .  2 
3 . PLANE STRAIN PROBLEM B O W E D  BY A FLEXIBLE CASE . . .  3 
4 . AXISYMMETRIC PROBLEM BOUNDED BY A FLEXIBLE CASE . . .  9 
5 . REVIEW OF CURRENT STATUS . . . . . . . . . . . . . . .  14 
6 . RECOMMENDATIONS FOR FUTURE WORK . . . . . . . . . . .  19 
1. INTRODUCTION 
Since the propellant grain must support itself, it is a structural 
member. Moreover, if the grain cracks before or during firing, then there 
results an unknown thrust-time profile and perhaps a catastropic failure. 
Part of the design process must therefore be an investigation of the struc- 
tural integrity of the grain. A fundamental part of this investigation 
is a stress-strain analysis of the propellant grain. On the basis of the 
thermal, pressure, and acceleration loads and the support conditions, the 
stress and strain distribution in the grain is to be determined. 
Due to the complex grain geometries and complex loads and supports, 
a rwmerical solution must be sought. A numerical method of sufficient 
versatility to handle these complexities is the finite element method. 
In particular for continuum problems, the stiffness version of the finite 
element method has proven to be the most successful. 
-2- 
2. STIFFNESS METHOD OUTLINE 
The essential feature of the stiffness method is the ability to 
analyze small elementary regions of the propellant grain independent of 
the total configuration and then combine these regions together to find 
the response of the total configuration. 
The stiffness method has similarities to the Ritz method. The 
equilibrium equations and stress boundary conditions are replaced by an 
equivalent variational principle, 87T(uJ = 0. The grain is divided into 
elements within which the displacement has approximately a specified form. 
For a plane strain problem a n d  an axisymmetric elasticity problem, tri- 
angular rings are used; in the axially symmetric shell problems, line 
elements are used. The unknowns are the displacements in the elements. 
They are expressed in terms of the node displacements which are then 
determined so that the functional T(u) is an extremum. 
While the theory underlying the stiffness method is sufficiently 
general to include three-dimensional problems, the programing and computer 
capacity needed to carry out practical calculations is not generally 
available. The current practice is to utilize two-dimensional programs. 
The plane strain program and the axisymmetric program together provide 
information about the state of stress which may be used for design. The 
programs of the Mathematical Sciences Corporation have the ability to 
include the effect of a flexible case bounding the grain. 
-3 - 
3. PLANE STRAIN PROBLEM BOUNDED BY A FLEXIBLE CASE 
The plane strain program provides a stiffness method solution to the 
following problem. 
by a circular shell. 
In Figure 1 is a finite two-dimensional region bounded 
Y 
Figure 1 
This is the typical situation encountered in a solid propellant grain 
analysis. 
The solid region is governed by the following system of equations, 
referred to rectangular Cartesian coordinates. 
Equilibrium: c T d P , p  + f C X = O  , (1) 
+ u  1 s (2) E,@ = %I, P, a 
Strain-Displacement: 
where 
Cap i s  the  s t ress  tensor,  
i s  the s t r a in  tensor,  
U Q  i s  the displacement vector,  
Cap'' 
a Q P  
T i s  the temperature change, 
f a  
are  the coeff ic ients  of e l a s t i c i t y ,  
a re  the coefficients of thermal s t r e s s ,  
is  the body force vector.  
The comma denotes p a r t i a l  d i f fe ren t ia t ion ,  and the  repeated indices are  
summed over the range 1, 2. The index 1 i s  associated with the x 
direct ion;  the index 2 with the y direct ion.  
A t  each point of the boundary of the sol id  with unit  normal n m ,  
e i the r  the displacement vector ug( or  the surface force 
must be given. 
the  body. The coeff ic ients  of e l a s t i c i t y  have the following symmetry: 
The temperature and body forces must be given throughout 
The coeff ic ients  aaP are symmetric. 
The she l l  i s  governed by the following equations. 
-5- 
Equilibrium : 
[ dx- N - O  , 
\ g + a + p = o  a dx , 
dv w E = - + -  
dx a ’ Strain -Displacement : 
N = D E  + N o  , r 
where 
N is the in-plane stress resultant, 
Q is the transverse shear, 
M is the moment resultant, 
p 
D is the shell membrane stiffness, 
K is the shell bending stiffness, 
No 
w is the normal displacement, 
v is the tangential displacement, 
a is the radius of the shell, 
x is a coordinate tangent to the shell. 
is the shell pressure distribution, 
is the thermal stress resultant, 
(5) 
-6- 
The stiffness method is applied to the solid using triangular elements 
and linear variation of displacement in each element. In Figure 2 is a 
typical element. 
Y 
Figure 2 
Within the element the displacements u and v are assumed to be of 
the form 
u = a x + b y + c  , 
v = d x + e y + f  . 
The six constants may be expressed in terms of the displacement at the 
three nodes i, j, k ; 
Note that the indices i, j, k are not free indices here. They denote 
definite node numbers. This process is repeated for all the elements. 
-7 - 
The result is a polyhedral displacement field defined by the values 
u v  at nodes. Putting this expression in the potential energy 
expression and minimizing results in 2n linear equations in the 2n 
unknown displacements at the n nodes; the coefficients are the stiffness 
matrix for the solid. 
i’ i’ ‘ 
The stiffness matrix is applied to the shell using line elements 
and the exact solution to the equations (5) through (7) for concentrated 
loads at the ends. In Figure 3 is a typical shell element. 
Figure 3 
The result is an element stiffness expression, 
-8 - 
where the loads F to F are associated with the six displacements; 
is a 6 x 6 element stiffness matrix. The process is repeated for the 
1 6 
whole shell. The element stiffnesses are summed together to obtain a 
stiffness matrix [k3 for the complete shell. 
Theoretically the stiffness matrix for the shell may be combined 
with that of the solid and inverted to obtain the displacement solution. 
However, in practice the result is a very poorly conditioned matrix that 
is difficult to invert. Instead the solution is obtained by an iterative 
scheme. 
The solution is obtained for the solid using zero values for shell 
displacements, and the reactions of the solid on the shell are calculated. 
The deflections of the shell are then determined using these reactions 
as loads. The solution is next obtained for the solid using the corrected 
values of the shell displacement to obtain new reactions. The procedure 
is then repeated until the error is as small as desired. 
The computer program used to carry out the numerical computations 
can handle up to 280 nodal points in a triangular breakdown of the solid 
and 21 nodal points for the shell. The program accounts for body force, 
temperature, and surface traction loads for the solid and temperature and 
surface traction loads for the shell. The material properties may be 
orthotropic and nonhomogeneous for both the shell and the solid. 
Using the given loads, nodal point coordinates, and the node groups 
bounding the elements, the program generates the stiffness matrix for the 
solid and the shell, solves the system of linear equations to find the 
displacement solution, and calculates the stresses. 
-9 - 
4 .  AXISYMMETRIC PROBLEM BOUNDED BY A nEXIBLE CASE 
The axisyrmnetric program provides a stiffness method solution to the 
following problem. In Figure 4 is a cross section of an axisymmetric body 
bounded by a flexible case. 
4 
I = r  
Figure 4 
This represents a problem often encountered in solid propellant grain 
analysis. 
The solid region is governed by the following system of equations. 
Equilibrium: 
+ u  1 Y (12) jyi Strain-Displacement : E = %(ui 
a 
where the def in i t ion  of the terms i s  the  same as before. The coma denotes 
covariant d i f fe ren t ia t ion  i n  the cyl indrical  coordinate system r ,  $, z ; 
the iQdices range from 1 t o  3.  I n  t h i s  ax isymetr ic  problem, there i s  no 
8 dependence i n  any of the variables,  and only the displacemects LI aid , 
w i n  the r and z directions respectively are nonzero. Circumferen- 
t i a l  or  hoop s t r a i n s  may ex i s t :  
is in general nonzero, 
A t  each point of the boundary (of the sol id  with uni t  normal n ) i 
e i the r  the displacement vector u o r  the surface force i 
must be specified.  Both u2 = v and F2 = F’ must be zero for  t h i s  I 
j 
axisymetr ic  problem. The temperature and body forces must be given I 
throughout the body. The coeff ic ients  of e l a s t i c i t y  have the same sym- 
metry as before. The axisymmetric s h e l l  i s  f i r s t  approximated by conical 
elements. In Figure 5 i s  a two-element approximation. 
-11- 
t 
I r  1 r  
Figure 5 
The conical fruatrum are g~verned by the equations obtained fron! the 
variational principle 
S T = O ,  
where 
JO 
a and p range over 1 and 2 ,  and repeated indices are summed. 
Z U n  
I + r  
Figure 6 
-12- 
a 
Using the notation of Figure 6 ,  the strains E, and X a  are given by 
du 
E 1 = x  Y 
2 d w  
r r Y (172 
u sina + w cosa E2 = 
sina dw 
x 2  = r ds . 
The other quantities are defined as follows: 
Cap = the membrane coefficients of elasticity, 
Da@ 
T P Q  = the thermal stresses, and 
pa = the surface tractions. 
= the bending coefficients of elasticity, 
For the stiffness method, displacements of the form 
S S u = (1 - -)u + - II 1 LU2 
and 
w = [ 1 - 3($ + $73 w1 + [ 3 ( i ) 2  - 2(;,?] w2 
are assumed for each element. 
stiffness matrices may be obtained. 
From this and equation (16) the element 
A s  in the plane strain program the shell and solid stiffness matrices 
are not combined. 
previously described. 
The solution is obtained by the iteration technique 
-13- 
e 
The computer program used to carry out the numerical computations 
ca? handle up to 280 nodal points in a triangular breakdown of the s o l i d  
and 100 nodal points for the shell. The program accounts for body forces, 
temperature, and surface traction loads for the solid and temperature and 
surface traction loads for the shell. The material properties may be 
orthotropic and nonhomogeneous for both the shell and the solid. 
Using the given loads, nodal point coordinates, atld the node groups 
bounding the elements, the progran generates the stiffness matrix for the 
solid and the shell, solves the system of linear equations to find the 
displacement solution, and calculates the stresses. 
-14- 
5 .  REVIEW OF CURRENT STATUS 
The programs as they now stand are sa t i s fac tory  for  the majority of 
problems encountered. However, they cannot handle materials tha t  a re  
incompressible o r  nearly incompressible nor motor cases tha t  are r e l a -  
t ive ly  f lex ib le  because the i t e r a t ion  solution technique i s  poorly 
behaved. 
The f i r s t  d i f f i c u l t y  l i e s  i n  the equations of e l a s t i c i t y  ra ther  than 
the s t i f f n e s s  method. 
Got va l id  for  incompressible materials ( V =  0.5). With displacements 
as unknowns, the  equations contain coef f ic ien ts  tha t  approach h f i n i t y  
as Poisson's r a t i o  V approaches 0.5. This behavior r e s u l t s  i n  a 
de te r iora t ion  of the numerical method as Poisson's r a t i o  approaches 0.5. 
The equations of e l a s t i c i t y  as  usually s ta ted  are  
A formulation su i tab le  for incompressible i so t ropic  nonhomogeneous 
materials i s  the following. 
(21.) 
i 
Equilibrium; Gij + f  = o  , 
' j  
Strain,-Displacement: E = %(ui, + u ) , ( 2 2 )  ij j , i  
P 
Mean Strain-Mean Stress: E - 3aT = (24) 
-15- 
a 
where 
a 
E i j  
i 
U 
a 
T 
f a  
P 
i s  the s t r e s s  tensor , 
i s  the s t r a i n  tensor,  
i s  the displaceaent vector ,  
a re  the Lame constants,  
i s  t he  coeff ic ient  of thermal expansion 
i s  the temperature change, 
i s  the body force vector ,  
i s  the  mean stress o r  hydrostat ic  pressure. 
/ 
The comma denotes covariant d i f fe ren t ia t ion ,  and t h e  repeated indices a re  
summed. The indices range over the  values 1, 2 ,  and 3. 
A t  each point of the boundary of the  sol id  with uni t  normal ni 
e i t h e r  the displacement vector ui or  the surface force 
(25) 
i s  specif ied.  The temperature and body forces a re  given throughour. t5e 
body. 
This f o m u l a t i o i  i s  sui table  for  a l l  values of Poisson's r a t i o  i2  
the in te rva l  0 < v < 0.5. It does not exhibi t  the i n f i n i t e  coef f i -  
c i en t s  a t  V = 0.5 tha t  the  usual formulation does. 
The s t i f fnes s  method has been applied t o  t h i s  set of equations by 
Mathematical Sciences Corporation. 
i n  MSC Report N o .  65-21-3. 
The resu l t ing  program i s  documented 
-16- 
a 
The second difficulty occurring with the relatively flexible cases 
lies in the implementation of the iterative solution technique. In 
Figure 7 the solid and the shell are shown as two separate 
the way they are considered by the program. Also shown is 
pressure load p . 
problems, much 
an internal 
-r 
Figure 7 
I =rr 
Consider the point P in Figure 7 and its radial displacemect u . If 
the boundary of the solid is held rigidly, then there is no displaceaent 
but a reaction stress CT . If the boundary of the solid is free, ther! 
there is a displacement u but no stress 0- . If the case experiences 
a stress 0- , then there is a displacement u . If there is no stress 
on the case, then there is no displacement. 
graphically in Figure 8.  
These statements are repeated 
-17 - 
a- 
A 
Solid Yse 
\ 
\ 
* U  
Figure 8 
The solution t o  the combined problem i s  the intersect ion of case and sol id  
curves. The i t e r a t i v e  solution proceeds as  follows. A solution for the 
sol id  is obtained using zero values for the she l l  displacements, point 1. 
This r e su l t s  i n  reaction loads Cl on the  she l l .  The deflectiozs of 
the she l l  a r e  then obtained due t o  these reactions,  point 2. The so l id  
i s  held a t  the corrected displacements 
point 3 .  The procedure i s  then repeated u n t i l  the error  i s  as small as 
desired. 
u2 and a solution obtained, 
This i t e r a t i v e  technique works very well for re la t ive ly  s t i f f  cases. 
A very high slope of the case curve i n  Figure 8 i s  a re la t ive ly  s t i f f  
case;  a f l a t t e r  slope i s  a re la t ively f lex ib le  case. I f  t h i s  process i s  
-18- 
repeated for f l a t t e r  and f l a t t e r  slopes of the case curve, a point w i l l  
be reached where the process no longer sp i r a l s  inward but diverges. 
By the t i m e  the solution a t  point 4 i s  established, there  i s  su f f i -  
c iegt  information available to  estimate the intersect ion of the sol id  and 
case curves. The program does t h i s  and continues. The estimation of the 
intersect ion is  repeated every t i m e  suf f ic ien t  information i s  available.  
For re la t ive ly  s t i f f  cases t h i s  process works very w e l l .  Below a 
cer ta in  case f l e x i b i l i t y  the i t e r a t ive  process begins t o  de te r iora te ,  
and the e r ror  cannot be reduced t o  a sat isfactory level .  
Wben the sol id  and the case s t i f fnes s  matrices a re  combined and a 
solution sought d i r ec t ly ,  j u s t  the opposite occurs. The re la t ive ly  s t i f f  
cases lead t o  numerical d i f f i c u l t i e s ,  while the re la t ive ly  f lex ib le  cases 
a r e  the eas ies t  t o  handle. 
-19- 
6. RECOMMENDATIONS FOR FUTURE WORK 
The application of the s t i f fnes s  method t o  t h e  incompressible problem 
i s  w e l l  on i t s  way t o  being an accomplished fac t .  
The i t e r a t i v e  technique has proven, i n  many problems, t o  be a very 
e f fec t ive  means t o  incorporate the f lex ib le  case i n t o  the analysis.  For 
the r e l a t ive ly  f l ex ib l e  case,  more work needs t o  be done with the i t e r a t i v e  
technique i n  order t o  expand the range of case f l e x i b i l i t i e s  t ha t  may be 
handled. It i s  recommended that  two avenues be followed. Only elementary 
numerical means have been used t o  find the  in te rsec t ion  of the case and 
so l id  curves i n  Figure 8 .  Continued e f f o r t  i n  t h i s  area should lead t o  
useful improvements. Second, a t  some point it w i l l  become more prof i tab le  
i n  terms of machine t i m e  and accuracy t o  combine the she l l  and case s t i f f -  
ness matrices and seek the solution d i rec t ly .  This should be included a s  
a program option. 
Although the application of the s t i f fnes s  method t o  the incompressible 
equations and t h e  considerations of a f lex ib le  case are  separate ideas, 
t h e i r  incorporation in to  a computer program cannot proceed independently. 
As  a r e s u l t ,  fur ther  study i s  required i n  order t o  incorporate the f l ex ib l e  
cases in to  the incompressible programs. 


Flexible case analysis for compressible solid propellant grain motors Status report

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server