Development of a three-dimensional nonlinear viscoelastic constitutive model of solid propellant

Development of a Three-Dimensional Nonlinear Viscoelastic Constitutive Model of Solid Propellant

Gyoo-Dong Jung

T&R Center, Agency for Defense Development Yusung. P. O. Box 35-5. Taejon. 305-600. Korea

Sung-Kie Youn

Bong-Kyu Kim

Department of Mechanical Engineering, Korea Advanced Institute of Science and Technology. 373-1, Kusung. Yusung. Taejon, 305-701. Korea

A three dimensional nonlinear viscoelastic constitutive model for the solid propellant is developed. In their earlier work, the authors have developed an isotropic constitutive model and verified it for one dimensional case. In the present work, the validity of the model is extended to three-dimensional cases. Large deformation, dewetting and cyclic loading effects are treated as the main sources of nonlinear behavior of the solid propellant. Viscoelastic dewetting criteria is used and the softening of the solid propellant due to dewetting is treated by the modulus decrease. The nonlinearities during cyclic loading are accounted for by the functions of the octahedral shear strain measure. The constitutive equation is implemented into a finite element code for the analysis of propellant grains. A commercial finite element package ‘ABAQUS’ is used for the analysis and the model is introduced into the code through a user subroutine. The model is evaluated with different loading conditions and the predicted values are in good agreement with the measured ones. The resulting model applied to analyze a solid propellant grain for the thermal cycling load.

Keywords: Solid propellant, viscoelastic constitutive model, dewetting, damage, cyclic load effects

Introduction

Composite solid propellants are considered as lightly cross-linked, long chain polymers filled with solid particles. These materials exhibit highly nonlinear viscoelastic response due to the damage process such as the Mullin’s effect, dewetting, nonlinear time-temperature effect and large deformation (Farris and Schapery, 1973). These microscopic changes cause the solid propellant to lose its initial stiffness and change its bulk behavior from incompressible to compressible behavior. The very complicated nonlinear behavior of solid propellants also comes from the features associated with time-temperature effects (Ravichandran and Liu, 1995). Recently, Jung and Youn (1999) proposed a constitutive model, which includes the effect of the large deformation, stress softening by dewetting in micro structure and nonlinearities in the cyclic loading on the overall stress-strain behavior. This model was experimentally verified on the various uniaxial loading conditions. In the present study, the model is extended to a general three dimensional state for the stress analysis of propellant grains. In extending the existing model, the cyclic loading damage functions are modified to improve stress prediction accuracy. A commercial finite element package ‘ABAQUS’ is adopted for analysis and the developed constitutive equation is implemented through a user subroutine, UMAT. The validity is tested using the strip biaxial specimens with circular holes subjected to different loading conditions.

Also, the model is used to analyze a solid rocket motor grain during thermal cycling after curing. The results are compared to those obtained from other geometrically nonlinear constitutive model.

Nomenclature

a_T= Time-temperature shift factor

c = Current filler volume fraction

C = Right Cauchy-Green deformation tensor

D_d = Deviatoric part of damage

D_v = Volumetric part of damage

E_ij = Green strain

E_rel = Stress relaxation modulus

F_ij =Deformationgradient

= Volume preserving deformation gradient

G = Shear modulus of composite or propellant

G₀ = Glassy shear modulus

G_rel = Shear stress relaxation modulus

G_c= Proportional constant, adhesion energy

= Octahedral shear strain

= First invariant of volume preserving right Cauchy-Green deformation tensor

= Second invariant of volume preserving right Cauchy-Green deformation tensor

J = Determinant of ^F

K = Bulk modulus of composite propellant

K₀ = Glassy bulk modulus

K_rel = Bulk stress relaxation modulus

R = Particle radius

S_ij = Second Piola-Kirchhoff stress

= Deviatoric stress

t = Time

T = Temperature

=Volumetric part of the initial elastic energy

V_o = Volume

e_ij = Infinitesimal strain

= Poisson’s ratio of composite

x = Time

s_ij = Cauchy stress

= Deviatoric part of the initial elastic energy

= Tangent modulus

The Constitutive Model Formulation

The Developed Constitutive Model

In the constitutive model developed by Jung and Youn (1999), it is assumed that the time characteristics of the free energy of solid propellant do not change by damage and the propellant is homogeneous and isotropic. In this case, the free energy function of the propellant can be regarded as the free energy of the undamaged propellant multiplied by a damage function. Damages considered are the dewetting phenomena and the nonlinearities inherent in the load-unload cycles. A viscoelastic dewetting criterion is derived from the law of energy conservation. The damage due to filler dewetting causes decrease in moduli of the propellant. The reduction in moduli is regarded as the results of the decrease in the effective filler volume fraction. The function of generalized octahedral shear strain is introduced as the medium of cyclic load effect. The large strain capability of solid propellant requires that the stress-strain relation be formulated in terms of the appropriate stress and strain tensors. For this, Simo’s model (1987) for stress-strain behavior is adopted and modified. Simo’s constitutive relation is modified by generalizing the elastic bulk response to viscoelastic one and introducing the volumetric and deviatoric damage functions, D_v D_d. The model was verified in the uniaxial stress field. While some models require many phenomenological nonlinear functions (Özüpek and Becker, 1996, Park and Schapery, 1997), this model has the advantages in that it reflects the realistic viscoelastic dewetting behavior as well as the effects of the properties of the individual constituent in the solid propellant. The resulting constitutive equation is of the following form.

, (2.1)

where

. (2.2)

, (2.3)

S(t) is the second Piola-Kirchhoff stress, J the determinant of the deformation gradient, the volume preserving part of Green strain E, the bulk relaxation modulus, the shear relaxation modulus, ,, and C the Cauchy-Green tensor. Also, and are the uncoupled volumetric and deviatoric parts of the initial elastic stored energy, D_vand D_d the volumetric and deviatoric damage functions respectively. Also

, (2.4)

are the reduced times and a_T the time-temperature shift factor. The nonlinearities during cyclic loading are accounted for by the functions of the octahedral shear strain measure(Özüpek and Becker, 1996). The damage functions are defined as follows:

. (2.5)

, (2.6)

where K_D and G_D are the decreased propellant moduli due to dewetting and f the function representing the nonlinearities during cyclic loading. So, K_D/K and G_D/G are the damages representing respective modulus decrease by dewetting. To calculate K_D and G_D, we need to know when debonding will occur and how much reduction of the moduli caused by dewetting. Viscoelastic dewetting criteria is developed to calculate the critical stress at which the particle will debond (Jung and Youn, 1999).

. (2.7)

Where c is the current filler volume fraction, S’_ij and s_ii the deviatoric and volumetric stresses, R the filler radius, G_c the adhesion energy between the filler and binder of solid propellant. From Eq. (2.7), we could say that large particle debond first under the same stress level. In the implementation, the filler particles of solid propellent are divided into several groups with different radii and statistical distributions. The debonded particles are eliminated from the propellant and replaced with voids of equal sizes. The dewetting causes loss of reinforcement in the propellant. Modulus changes based on filler volume fraction are calculated using the Farber-Farris equations(Farber-Farris 1987):

(2.8)

where G is the shear modulus of propellant, Poisson’s ratio of propellant, K the bulk modulus of propellant, and K_i, G_i the moduli of filler. The damage caused by dewetting can be evaluated by the following procedure. At first, stress is calculated for given deformation history. And then dewetting occurrence is determined by Eq. (2.7) for current group of undewetted filler particles with largest radius. And then if dewetting occurs, the dewetted particles in the current group are removed to yield new filler volume fraction. Then K_D and G_D are calculated using Eq. (2.8). Finally, damages due to dewetting are calculated by Eq. (2.5) and Eq. (2.6). Above procedure is repeated until no dewetting occurs.

The cyclic load effects are characterized as the rapid decrease of stress during unloading and the large amount of hysteresis in load-unload cycles. The cyclic load effects depend on current shear strain and maximum shear strain during the loading history, which is known through experiment (Swanson and Christensen, 1983). Thus the nonlinearities during cyclic loading are accounted by the following function ¦ which depend on the generalized octahedral shear strain measure (Özüpek and Becker, 1996).

(2.9)

where , and are the first and second invariant of . And f_u and f_r are the functions obtained from unloading and reloading test. Also, f_r/u = f_r / f_u, and is the value at the end of unloading.

Improved Cyclic Loading Damage Function

The unloading and reloading functions, f_u and f_r are shown in Fig. 1. These functions are obtained from the uniaxial cycling test with constant strain amplitude (0~20% strain) at 50%/min, 20⁰C by comparing the measured and predicted stresses(Jung and Youn, 1999). The predicted stresses were well matched with the measured ones for the uniaxial loading-unloading test. For the reloading after partial unloading condition as in Fig. 2, the function f evaluated by Eq (2.9) is shown as the triangular symbols in Fig.3. In Fig. 3, the values of f turn out to be impractically small or large when the load direction is reversed. This does not reflect the real behavior of propellant and make the convergence rate slow when used in the finite element analysis. Therefore, the form of f in Eq. (2.9) is modified for the reloading or unloading from reloading conditions as follows:

(2.10)

In Eq. (2.10), the values of f is considered to vary lineary from f _u to f_g during . The final form of f is obtained as in Eq. (2.11).

(2.11)

where

(2.12)

This form of f results in smooth transition of f between f_u and f_r when the load direction is reversed. This is shown as the circular symbols in Fig. 3. When this improved form of f is applied for the uniaxial complex multiple load test conducted by Jung and Youn (1999), the comparisons of the predicted and observed stress-strain behaviors are shown in Fig. 4. The modified form of f, Eq. (2.11) shows more accurate stress predictions than the previous one, Eq. (2.9), at cyclic loading condition.

Computation of the Constitutive Equation

For the analysis of three-dimensional finite element analysis, the constitutive equation is implemented in the user subroutine UMAT of ABAQUS code. This can be done by updating the stresses at the end of the increment and providing the tangent stiffness using the presented constitutive model (ABAQUS user’s manual, 1996).

Stress Formulation

For finite strain applications, the interface for subroutine UMAT requires Cauchy stress components as stress measures. The second Piola-Kirchhoff stress in the constitutive equation Eq. (2.1) is converted to Cauchy stress by the relation

. The constitutive equation in terms of Cauchy stress becomes as follows:

(3.1)

where and is volume preserving deformation gradient. The calculation of the convolution integrals, H and P, is accomplished by using a numerical algorithm developed by Taylor et al.(1970). By using the Prony series representations of the relaxation moduli and approximating the integrals of each series terms, H and P are arranged in the following forms.

(3.2)

, (3.3)

where are modulus ratios, , and the right subscripts mean time steps, that is, means the value of at time t_n+1.

Tangent Stiffness Formulation

In ABAQUS(ABAQUS theory manual, 1996), Kirchhoff stress increment d_tt is represented by spin tensor W and tangent stiffness as

(3.4)

By letting as ,using the relation between 2nd Piola Kirchhoff stress S and Cauchy stress s and introducing time derivative of Kirchhoff stress, t, is represented by the following relation

, (3.5)

where

, (3.6)

is obtained by differentiating Eq. (2.1) by Green strain E as

, (3.7)

where L_dev is

(3.8)

and

(3.9)

, (3.10)

, (3.11)

. (3.12)

Therefore, substituting Eq. (3.7) for Eq. (3.6) and simplifying gives

, (3.13)

where is the identity tensor of order 4. _v and are represented as:

, (3.14)

(3.15)

where , and.

The unsymmetric feature of in Eq. (3.14) produces the unsymmetric stiffness matrix in the finite element analysis. For the sake of computational efficiency, symmetric approximation of , , is used. Substituting Eq. (3.13) for Eq. (3.5), the final form of tangent stiffness can be written as follows:

.(3.16)

In this study, Mooney-Rivlin hyper elastic model is used for the deviatoric free energy function and second degree polynomials of J for the volumetric free energy function.

. (3.17)

. (3.18)

Therefore and yield followings.

, (3.19)

. (3.20)

As stated above, the stress and tangent stiffness components are calculated in UMAT from Eq. (3.1), (3.16). The algorithm for the implementation of the model in a computer program to predict stress-strain behavior can be summarized as follows :

Initial propellant properties and configuration:

1. Specify mechanical properties of the solid propellant, matrix, filler, void and adhesion energy

2. Specify initial statistical distribution of filler particle size.

Deformation of solid propellant:

3. Calculate the cyclic loading damage function by Eq.(2.11)

4. Calculate the viscoelastic stress by Eq. (3.1)

5. Check, with Eq. (2.7), whether dewetting occurs at this stress for the largest particles

6. If dewetting do not occurs, calculate tangent modulus by Eq.(3.16) and go to ABAQUS main routine.

7. If dewetting occurs, reduce the dewetted filler volume fraction and increase the void volume fraction by the same amount.

8. Calculate the reduced relaxation modulus and damage functions by Eq.(2.8, 2.5, 2.6, 3.16).

9. Reduce the current filler volume fraction by the dewetted filler amount.

10. Go to STEP 4 and calculate the reduced nonlinear viscoelastic stress by Eq.(3.1).

Laboratory Experiment

The material used in this study is an HTPB solid propellant with 76% particle volume fraction of AP(ammonium perchlorate) and Al(aluminum) powder. The initial propellant properties and configuration required in the simulations were measured (Jung and Youn,1999). Moduli were obtained from uniaxial, shear relaxation tests and Pockerchip constant rate tests at –90~60^oC as shown in Fig. 5. Adhesion energy between binder and AP was obtained from 180^o peel test. Moduli of the void were determined by matching the uniaxial constant rate stress-strain curve at 20^oC, 50%/min. Statistical distribution of filler particle size was measured by a Malvern Series Particle sizer as shown in Fig. 6. The methods to obtain these parameters and the results are explained in detail at the previously paper(Jung and Youn,1999).

The biaxial strip tests are conducted to verify the validity of the three-dimensional constitutive model. The strip biaxial specimen is rectangular in shape and bonded to wooden-tab as illustrated in Fig. 7 (Schapery, et. al, 1975). The specimen contains a circular hole to produce severe stress gradients. A large block of the material is bonded to wooden tab with the polyurethane adhesive and specimens are mill machined out of the block. The specimens are stored in a desicator (RH<10%) to minimize the effects of humidity. Constant rate tests, complex load tests are conducted. Four specimens are tested at one test condition and stress responses are averaged. All tests are performed in an Instron 1122 and gas dilatometer in the humidity controlled room (RH<30%). Prior to testing, the specimens are conditioned at each test temperature for more than 1 hour .

Finite Element Analysis

Finite element mesh model of the specimen is shown in Fig. 8. One-eighth specimen is modelled using the symmetry conditions. Twenty noded, three-dimensional isoparametric elements are used to model the specimens. Finite element analyses are conducted for the loading conditions tested. The initial time increment equivalent to the 0.5% global strain is chosen for the analysis. Auto time increment and the convergence criteria in ABAQUS are applied.

Application Results and Discussion

The Constant Rate Test

The biaxial specimen is loaded at constant strain rates, 66.7%/min and 0.67%/min and temperatures, 60° C, 20° C, and -40° C. The comparisons between the predicted and measured responses are shown in Fig. 9, Fig. 10 and Fig. 11. Global stress and volume dilatation against global strain are plotted. Global strain means the displacement divided by the specimen height and global stress the load divided by the minimum cross section area of the specimen. A specimen with a hole creates a strain concentration around the perimeter of the hole, the crack initiate at either side of the hole at the global strain lower than the uniaxial test. The maximum values of each measured strain represent the crack initiation points. The finite element analysis is conducted to the global strain of 30%. For some loading conditions, the solutions do not converge up to the 30% global strain. These nonconvergent point corresponds to the one that the maximum stress in the specimen reaches the tensile strength of the propellant. In this case, the norm of tangent stiffness approaches zero and the material instability occurs. In reality, the crack initiates at this nonconvergent stress values. In general, the predicted stress responses are in good agreement with the measured ones for the various strain rates and temperatures. So, the current model is verified for the biaxial geometry in addition to the uniaxial one(Jung and Youn, 1999). This demonstrate the validity of the model which treat the softening of the propellant by the modulus decrease due to viscoelastic dewetting. The calculated dilatations are distributed somewhat lower than the measured values. However their trends are shown to be consistent.

Constant Cycling Test

Cyclic load is applied to the biaxial specimen at 20° C. The test is done with 5 cycles at global strain levels of 0~ 17.6% and strain rate 66.7%/min. There is a good agreement between the predicted and measured stress as shown in Fig. 12. Thus the unloading and reloading functions obtained from one uniaxial cycling test condition can be used for the biaxial geometry. The applicability of cyclic load damage function is verified.

Similitude Test

The biaxial specimen is loaded at strain rate 1.67 %/min to the global strain level 12.3 % and then allowed to relax for 3 hrs. The relaxation is repeated at the global strain of 3.3 % for 2 hrs and then the specimen is loaded to failure at 6.7 %/min. While the relaxation response is somewhat overpredicted, the loading portions are well predicted as shown in Fig. 13. This results show the applicability of the model for the relaxation and loading conditions.

Straining and Cooling Test

These tests are intended to verify the model for the loads that the solid propellant grains usually experience. The biaxial specimen is cyclic loaded with the temperature lowering at a constant rate. The first test is done with the increasing strain amplitude 0~ 22.3% at 1.67 %/min, while the temperature is changed at –1^oC/min from 38 ^oC as shown in Fig. 14. Another cyclic test with subsequent loading and unloading at different strain levels is also conducted with the same strain rate of 1.67 %/min while the temperature is changed at -1 ^oC /min from 27 ^oC as shown in Fig. 15. As shown in Fig. 14, 15, the model somewhat underpredicts the magnitude of the stress. But the predicted stresses are close to the measured ones. Also, as shown in Fig. 15, the stresses are predicted well for the reloading or unloading from reloading conditions. So, the proposed cyclic load damage functions are verified for the complex loading conditions and the model seems to be reasonable for the complex cyclic load and cooling test.

As mentioned above, the constitutive model predicts reasonably well the solid propellant behaviors under multi-axial and complex loading conditions which consist of constant loading, cycling loading, relaxation and straining and cooling. The model uses the damage functions developed from a few uniaxial tests(Jung and Youn, 1999). The accuracies of the biaxial stresses and dilatations are somewhat lower than the uniaxial cases(Jung and Youn, 1999). However the predicted values are very close to the measured ones. So the current model can be used for the stress analysis of the solid propellant grain.

Comparison with other Models

The comparisons with other models are impossible because the filler particle distributions and adhesion data needed in current model are not available in the references of other models. Also, the implementations of other models require some large amount of experiments and programming to characterize the damage functions. So, the comparisons are made for the results predicted by the constitutive model used in ABAQUS, Simo’s model used in TEXPAC code(Becker and Miller, 1989) and current model. The constitutive model in ABAQUS and Simo’s model suitably generalized the linear viscoelastic hereditary integral to finite strain formulation. These two models can be used for viscoelastic materials, which is homogeneous and capable of large deformation without damage. But it is inappropriate to apply these models to composite viscoelastic materials like solid propellant where damages are unavoidable when in use. For the loading conditions of constant cycling test(for the case of Fig. 12), the predicted and measured responses are shown in Fig. 16. The predicted results by ABAQUS and Simo’s model are almost the same. As the strain increases, the predicted stress response by ABAQUS and Simo’s model overestimates the measured values and the overestimation become more conspicuous for the unloading. Th reason is that material nonlinearities due to dewetting and cycling loading are not considered in the constitutive model used in ABAQUS and Simo’s model in TEXPAC. The original Simo’s model(1987) considers the material nonlinearity due to only strain softening. However, the model does not reflect the complicated nonlinear behavior of solid propellants as particulate composites and the one in TEXPAC considers only large deformation.

The Simulation of Solid Rocket Motor Grain

We shall consider the thermal cycling period after propellant curing and examine the response of the motor grain. A simplified 5-mammel shaped grain cross section is considered. The grain is assumed to be in plane strain condition and the web fraction is 0.57. The motor is thermal cycled from 60^oC, the curing temperature, to –40^oC, the minimum temperature the propellant is expected to experience. The cooling and heating occurs individually over a period of 24 hrs as shown in Fig. 17. The relaxation moduli of Fig. 5 are used. The liner is a nearly incompressible viscoelastic material and the casing is made of steel. The thermal expansion coefficient of the propellant is 1.5´10^-4oC^-1 and that of the steel case is 1.0´10^-5oC^-1. Mesh model and the maximum principal stress distribution at 24 hrs are shown in Fig. 18. One-tenth grain cross section is modelled using the symmetry conditions. Eight noded, plane-strain isoparametric elements are used to model the grain. The largest maximum principal strain and stress are obtained at point A. The comparisons are made for the maximum principal strain and stress at point A that are predicted by the Simo’s models used in TEXPAC code and current model as shown in Fig. 19 and 20. As shown in Fig. 19 and 20, dewetting occurs at the lower temperatures. At 24 hrs, the maximum strain predicted by current model is 15.6% and the one obtained from Simo’s model is 15.2. But, the maximum stress is reduced from 14.4 kg/cm2 to 13.6 kg/cm2 at the same time. These is due to the fact that the modulus at point A is decreased by dewetting. But the strain is caused mainly by the difference of the thermal expansion coefficients between the propellant and steel case. So, the stress decrease is much higher than the strain increase. After 24 hrs, cyclic loading effects begin to take effect. The stresses predicted by current model become lower than those predicted by Simo’s model at unloading and reloading conditions. So, the current model produces lower stresses and slightly higher strains than Simo’s model. These results show the importance of taking such material nonlinearities.

Conclusions

A three-dimensional nonlinear viscoelastic constitutive model of composite solid propellant is developed and tested for biaxial loading conditions. In the model, the cyclic loading damage function of Özüpek and Becker that is based on the octahedral shear strain measure is modified to improve the stress prediction accuracy. The constitutive equation is implemented into a finite element analysis code for the analysis of biaxial specimens. A commercial finite element package ‘ABAQUS’ is adopted for analysis and the developed constitutive equation is implemented through a user subroutine, UMAT. The constitutive model has been tested for biaxial complex loading conditions and the predictions have been compared with the experiments. The results show that the model predicts the broad range of the propellant behaviors with reasonable accuracies. Also, the importance of the material nonlinearities is demonstrated through the comparisons with the constitutive models in ABAQUS and Simo’s model for certain biaxial loading conditions. Also, the model is used to analyze a solid rocket motor grain during thermal cycling after curing. Therefore the developed constitutive model can be effectively used in the three dimensional finite element analysis for solid propellant grains.

Presented at DINAME 99 – 8th International Conference on Dynamics Problems in Mechanics, 4-8 January 1999, Rio de Janeiro. RJ. Brazil. Technical Editor: Hans Ingo Weber.

• Becker, E. B. and Miller, T., 1989, "User's manual for the TEXPAC computer code", MSI, Austin, Texas.
• Farber, J. N. and Farris, R. J., 1987, "Model for prediction of the elastic response of reinforced materials over wide ranges of concentration" J. Appl. Polym. Sci. Vol. 34, pp. 2093-2104.
• Farris, R. J. and Schapery, R. A., 1973, "Development of a solid rocket propellant nonlinear viscoelastic constitutive theory" AFRPL-TR-73-50.
• Hibbit, Karlsson & Sorensen, Inc., 1996, "ABAQUS user’s manual", Ver. 5.6.
• Hibbit, Karlsson and Sorenson Inc., 1996, "ABAQUS theory manual", Version 5.6.
• Jung, G. D. and Youn, S. K., 1998, "A nonlinear viscoelastic constitutive model of solid propellant", Int. J. Solids Structures, Vol. 36, pp. 3755-3777.
• Özüpek , S. and Becker, E. B., 1996, "Constitutive equations for solid propellants", Int. J. Engng Mat. Tech. Vol. 119, pp. 125-132.
• Park, S. W. and Schapery, R. A., 1997, "A viscoelastic constitutive model for particulate composites with growing damage", Int. J. Solids Structures, Vol. 34, pp. 931-947.
• Ravichandran, G. and Liu, C. T., 1995, "Modelling constitutive behavior of particulate composites undergoing Damage", Int. J. Solids Structures, Vol. 32, pp. 979-990.
• Taylor, R. L., Pister, K. S. and Goudreau G. L., 1970, "Thermomechanical analysis of viscoelastic solids", Int. J. Numer. Methods Eng., Vol. 2, pp. 45-59.
• Schapery, R. A. et. al., 1975, "Solid propellant mechanical behavior manual", CPIA Pub. 21.
• Simo, J. C., 1987, "On a fully three-dimensional finite strain viscoelastic damage model: formulation and computational aspects", Computer Methods in Applied Mechanics and Engineering, Vol., 60, pp. 153-173.
• Swanson, S. R., Christensen, L. W., 1983, "A constitutive formulation for high elongation propellant", J. Spacecraft, Vol. 20, pp. 559-566.

Publication Dates

• Publication in this collection
  18 Dec 2000
• Date of issue
  2000