Abaqus-umat Program Example

ABAQUS IMPLEMENTATION OF CREEP FAILURE IN POLYMER MATRIX COMPOSITES WITH TRANSVERSE ISOTROPY

A Thesis Presented to The Graduate Faculty of The University of Akron

In Partial Fulfillment of the Requirements for the Degree Master of Science

Fengxia Ouyang December, 2005

ABAQUS IMPLEMENTATION OF CREEP FAILURE IN POLYMER MATRIX COMPOSITES WITH TRANSVERSE ISOTROPY

Fengxia Ouyang

Thesis

Approved:

Accepted

________________________________ Advisor Wieslaw Binienda

_______________________________ Dean of the College George K. Haritos

________________________________ Faculty Reader Pizhong Qiao

______________________________ Dean of the Graduate School George R. Newkome

________________________________ Department of Chair Wieslaw Binienda

______________________________ Date

ii

ABSTRACT

Polymer Matrix Composites (PMC) are increasingly favored in structural applications for their light weight and durability. However, the numerical modeling of these materials poses several challenges. This is primarily due to the highly anisotropic nature of the creep exhibited by these materials above the glass transition temperature. Also, the damage and failure of the material is of particular interest to designers using the PMCs. Recently, a sustained effort has been to provide the designer with large computer codes containing comprehensive constitutive equations, often equipped with large number of internal variables and with the most general mathematical forms, for use in structural design analysis. Although elegant and useful, such constitutive laws are often expensive in implementation. Specially for early stages of the design, a quicker way of estimating complicated PMC behavior is needed. In this work, the constitutive material law by Robinson and Binienda (2001) [1,2] is utilized for such an approach. The model is successful in describing polymer matrix composite (PMC) materials having long or continuous reinforcement fibers embedded in a polymer matrix. Although the material law includes a single scalar parameter to describe the damage, it retains the essential material behavior. The material law is implemented computationally as a user defined subroutine (UMAT) in a commercially available FEA code (ABAQUS). The material parameters iii

are obtained from experiments of thin-walled tubular specimens reinforced with unidirectional, helical fibers at an angle θ = 90 o , 60 o and 45 o under tensile and shear loading . The model correctly predicts the relation between logarithmic creep rate and logarithmic stress. The user subroutine has robust convergence properties. The creep strain rate and the effect of damage on the creep strain rate are presented for the benchmark problem of a square plate with a circular hole at the center and pressure vessel. The effect of fiber orientation on the durability of the square plate and pressure vessel under damaging loads, is studied.

.

iv

TABLE OF CONTENTS

Page LIST OF TABLES…………………………………………………………….…...…..vii LIST OF FIGURES........................................................................................................viii CHAPTER I. INTRODUCTION…………………………………………………….……...…1 II. THE ANISOTROPIC VISCOELASTIC MODEL………………………..…....3 III. THE ISOTROPIC POWER LAW MODEL…………………………….……....9 IV. USER SUBROUTIN UMAT………………………………..……………...….11 4.1

Introduction of UMAT …………………………………………..…11

4.2

Implementation of UMAT………………………………….……....13 4.2.1

Jacobian matrix for plane stress element and shell element……………………….……..17

4.2.2

Newton Ralphson Method………………………...…...20

V. PLANE STRESS ELEMENT IMPLEMENTATION ……………………..…23 5.1

Single element test……………………………….…..………….…23

5.2

Square plate ……………………………………..…………….......34

5.3

Plate with a hole……………………………………………………37

v

VI. SHELL ELEMENT IMPLEMENTATION……………………………….…..49 6.1

Theory of shell element……..……………………………………...49

6.2

Shell element in ABAQUS……………………..………………….50

6.3

Implementation of shell element………..……………..…………..51

VII. CONCLUSION……………………………………………………...….……..62 REFERENCES………………………………………………………….………….….64 APPENDICES…………………………………………………………………….…...66 APPENDIX A UMAT………………………………………….………......67 APPENDIX B INPUT FILE FOR PRESSURE VESSEL………….…..….77

vi

LIST OF TABLES

Table

Page

1 Fiber orientation 90 deg under tension ………………………………….………..28 2 Fiber orientation 45 deg under tension.……………….…...…………….………..28

vii

LIST OF FIGURES

Figure

Page

2.1 Illustration of the isochronous damage function in the normal stress (y axis). P=1 is the linear form……………………………………….……5 2.2 Thin wall tube under tension and torsion loading. (a) thin wall tube with fiber orientation. (b) typical experimental creep data under axial and shear loading……………………………………………………………..7 2.3 Non-dimensional log creep rate vs. log stress for the complete exploratory data set. Tensile and shear data are shown as indicated. Fiber angles θ = 60 o and 45 o and shear data for θ = 90 o shifted to form a master curve…………………………….……………………………...…..8 . 4.1 UMAT Loop ………………………………………………………………….….16 4.2 Plane Stress element…………………………………………………………...…17 4.3 Plane stress transformation from local coordinate (x-y) to global coordinate (1-2)………………………………………………………….……......18 5.1 One element model problem set up with fiber orientation (a) creep loading and (b) constant displacement loading…...……………………24 5.2 Comparison between power law and Robinson creep model for isotropic case under a constant load of 45 MPa when the scalar damage variable is not included in the Robinson creep model. (a) Time evolution of creep strain in Y direction. (b) Time evolution of creep strain in Y direction. (c) Time evolution of creep strain in XY direction……………………….….…..29 5.3 Time evolution of creep strain in the direction of loading (without damage) for orientations of the fiber 90 deg at 45MPa load…………….…….....30

viii

5.4 Time evolution of creep strain in the direction of loading (without damage). (a) for different orientations of the fiber (0, 45, 90 deg) at 45MPa load. (b) for 45 deg fiber orientation at varying loading 20, 45 and 100PMa…………………………………………………….….….…..31 5.5 Time evolution of creep strain with damage effect.(a) comparison when damage evolution is included under tensile load 60 MPa. (b) comparison between tensile loads 70, 75, 80 MPa with fiber orientation 45 deg. (c) comparison between fiber orientation 0, 45, 90 deg under shear loads 40 Mpa………………………………………….32 5.6 Time evolution of stress relaxation in the direction of Loading (with damage)………………………………………………………..….33 5.7 Geometry and boundary condition of square plate problem…………………..….35 5.8 Time evolution of creep strain with different orientations of the fiber (0, 45, 90 deg) at 46MPa load in the direction of loading (without damage). (a) time evolution of strain at 1 direction (b) time evolution of shear strain ……………………………………………….36 5.9 Left, geometry and right, 2D quarter symmetry model for the square plate with a hole problem…………………………………………….…...37 5.10 Contour plot of (a) stress distribution of the power law for isotropic material after elastic step, (b) stress distribution of the Robinsons’ model law for isotropic material after elastic step, (c) stress distribution of the power law for isotropic material after creep step (5 hours), (d) stress distribution of the Robinsons’ model law for isotropic material after creep step (5 hours)………………………..…...…..40 5.11 Stress distribution along the hole (quarter circle) going in a counter clockwise direction for the proposed Creep Damage Model as and for isotropic Power Law Creep model. (a) stress distribution after elastic step (b) stress distribution after 5 hours…….............…………………..…...….…….41 5.12 Contour plot of Creep strain distribution after 5 hours creep response. (a) Robinson Damage Model and (b) Isotropic Power Law……..………....…..42 5.13 Creep strain along the hole (quarter circle) going in a counter clockwise direction for the Robinson Damage Model and for isotropic Power Law Creep model (after 5 hours)……….………………….…..43

ix

5.14 Comparison between Analytic solution and Robinson model (FEA) in elastic step. (a) R=20mm and R/L=13%. (b) R=50mm and R/L=33%…........44 5.15 Creep / time response along the hole (quarter circle) going in a counter clockwise with damage effect under axial stress 10MPa with R=20mm, 50mm, 80mm, L=152.4mm, fiber orientation 45 deg……….....45 5.16 Creep / time response along the hole (quarter circle) going in a counter clockwise with damage effect under axial stress 10MPa with R=50mm, L=152.4mm, fiber orientation 45 deg and 90 deg………….…..46 5.17 Stress concentration zone with damage effect under axial stress 10Mpa fiber orientation 45 deg. (a) stress concentration zone with R/L=0.1. (b) stress concentration zone with R/L=0.3. (c) stress concentration zone with R/L=0.5…………………………..……...….47 5.18 Stress compression zone with damage effect under axial stress 10Mpa fiber orientation 45 deg. (a) stress compression zone with R/L=0.1. (b) stress compression zone with R/L=0.3. (c) stress compression zone with R/L=0.5……………………………...…….....48 6.1 Thin wall tube under tension and torsion with fiber orientation,…………....……53 6.2 Time evolution of creep strain under a constant tensile load 45 MPa of thin-walled tube for different fiber orientations (0, 45, 90 deg), (a) time evolution of maximum principle strain. (b) time evolution of shear strain……………………………………………........54 6.3 Pressure vessel, geometry and mesh and path……………………………….…...57 6.4 Path plot of time evolution of Maximum Strain along path1 of the vessel under inside pressure 0.5Mpa with damage evolution (a) fiber orientation 45 deg. (b) fiber orientation 60 deg. (c) fiber orientation 90 deg ………………………………………………….…...58 6.5 Maximum Strain along path1 of the vessel with fiber orientation 0, 45, 90 deg under inside pressure 0.5Mpa with damage evolution after 10 hours…………………………………………………………………......59 6.6 Time evolution of Maximum Strain along path2 of the vessel with fiber orientation 0, 45, 90 deg under inside pressure 0.5Mpa with damage evolution………………………………….…..…..60 x

6.7 Time evolution of Maximum Principle Strain of the pressure vessel with fiber orientation 60, 90 deg under pressure 1Mpa with damage evolution. (a) fiber orientation 60 deg, (b) fiber orientation 90 deg……………....61 -

xi

CHAPTER I INTRODUCTION

The present work details the development of a computational material model for polymer matrix composites.

These materials are increasingly found in structural

applications for their light weight and durability. However, the numerical modeling of these materials poses several challenges. This is primarily due to the highly anisotropic nature of the creep exhibited by these materials. Also, the damage and failure of the material is of particular interest to designers using the polymer matrix composites. Over recent decades, a sustained effort has been to provide the designer with large computer codes containing comprehensive constitutive equations, sometimes embodying several state variables and the most general mathematical forms, for use in structural analysis in support of design. Although this approach may be appropriate in the final stages of design and where complex histories of stress and temperature are involved, it can be too complicated for many design applications, particularly in the early stages of design, and is seldom the path taken by the designer. The constitutive material law based on a Norton /Bailey type of creep law by Robinson and Binienda (2001) [1,2] is utilized for this computational model. The model is successful in describing polymer matrix composite (PMC) materials having long or continuous reinforcement fibers embedded in a polymer matrix. The objective 1

in choosing a numerical model based on this type of material law is to allow simplicity and utility (on behalf of the designer), while still retaining the essence of the actual material behavior. When strongly reinforced PMC materials and structures operate in the creep range of their polymer matrix (at or near Tg) they undergo time-dependent deformation and eventually fail. Substantial resistance to creep deformation and damage is achieved in materials where long or continuous fibers are embedded in the polymer matrix. By design, much of the load in such materials is carried by the strong fibers that creep very little, if at all. Evidently, design engineers need quantitative tools for predicting the response of PMC materials and structures as a basis of achieving optimal structural designs, e.g., optimal fiber configurations. Although it is often observed that highly anisotropic creeping materials exhibit inelastic compressibility, cf., Robinson and Binienda (2001) [1], creep deformation models do not commonly include dependence on hydrostatic stress; this dependence is a principal feature of the material law that is implemented in this work. The model incorporates a dissipation potential function that is taken to depend on appropriate invariants of stress and material orientation consistent with transverse isotropy. Failure and damage are introduced in the model with a Monkman / Grant type of relationship [2,3,4,5]. The particular viscoelastic model used here is as described by following chapter.

2

CHAPTER II THE ANISTROTROPIC VISCOELASTICITY MODEL

The transversely isotropic viscoelasticity model of concern derives from the more general anisotropic deformation /damage model proposed in Robinson et al.(1992). The present model is an extension of that earlier work in the sense that hydrostatic stress is taken into account in the deformation response. The extension follows Robinson et al. (1994) and Robinson and Binienda (2001)b. The viscoelasticity model has the form e& ij = ε& ije + ε& ij

ε& ij ε& o

=

& =− ψ

where

ε eij

(2.1)

3 n −1 Γij 1 Φ 2 σo ψ n

(2.2)

1 1 ∆ν m ( m + 1) t o ψ

(2.3)

denotes the components of elastic strain, e& ij and ε& ij denote components of the

total and creep (viscous) deformation rates, respectively. σ ij are the components of Cauchy stress; σ o is a reference stress; ε& o , n, m, t o and ν are material parameters; Φ is a dissipation potential function; ∆ is an isochronous damage function and the scalar ψ is the Kachanov continuity. ψ = 1 corresponds to an undamaged material element; ψ = 0 indicates its total loss of load carrying capacity.

3

For transverse isotropy, the functions Φ and ∆ are taken to depend on the following invariants of stress σ ij , deviatoric stress sij and a material orientation tensor D ij

. J2 =

1 s ij s ji 2

J o = D ij s ji

I = σ ii

I o = D ij σ ji

J = D ij s jk s ki

(2.4)

The dissipation potential function is Φ=

and

1 σo

Γij = ∂Φ

3[J 2 − ξ(J − J 2o ) − (ζ − η)J 2o +

∂σ ij

1 (ζ − 4 η)I 2 ] 9

= sij − ξ ( Dik s kj + sik Dkj − 2 J o Dij ) −

1 2 2(ζ − η ) J o ( Dij − δ ij ) + (ς − 4η )I∂ ij 3 9 Φ

(2.5)

(2.6)

is a positive, homogeneous function of degree unity in stress assuring that the

representation (1)-(6) is dissipative. The anisotropy parameters ξ, η and ζ are subject to various inequalities based on physical limitations that are specified in Robinson and Binienda (2001) [2]. Also the isochronous damage function is specified as: ∆ = ∆ (σ ij , D ij ) = ∆ ( N , S) p ⎡⎛ N ⎞ p ⎛ ⎞ ⎤ ⎟ + ⎜ αS ⎟ ⎥ ∆ ( N, S) = ⎢⎜ ⎜σ ⎟ ⎥ ⎢⎜⎝ σ o ⎟⎠ ⎝ o⎠ ⎦ ⎣

(2.7) 1/ p

(2.8)

in which, the invariant N specifies the maximum tensile stress normal to the local fiber direction, and the invariant S denotes the local maximum longitudinal shear stress as: N=

1 1 (I − I o ) + J 2 + J 2o − J 2 4

(2.9)

S = J − J o2

(2.10) 4

The angular brackets in (9) are the MacCauley brackets Fig 2.1 illustrates the isochronous damage function in N- S space for a general power law form as well as the special case when the function is linear in N and S (P=1).

Fig 2.1 Illustration of the isochronous damage function in the normal stress (y axis). P=1 is the linear form.

5

The material parameters are obtained from experiments of thin-walled tubular specimens reinforced with unidirectional, helical fibers at an angle θ = 90 o , 60 o and 45 o under tensile and shear loading (Fig 2.2). The complete exploratory data set is plotted in Fig 2.3. The correlation of the theoretical model and experimental data (solid lines) and of creep response under two of the untested natural stress states (TS)-(dotted line) and (LN)-(dashed line) validates the physics behind the proposed constitutive equation. As a definitive measure of the strength of anisotropy, the creep rate under tensile stress along the fiber direction (LN) is predicted as being less than one thirtieth (1/30) of that under transverse stress (TN).

6

2

τ

T F

F θ

σ

σ

1 T

τ

3

(a)

axial and shear creep strain (%)

0.0040 0.0035 0.0030 axial strain -shear strain

0.0025 0.0020 0.0015 0.0010 0.0005 0

1

2

3

4

5

6

7

time (hr)

(b)

Fig 2.2 Thin wall tube under tension and torsion loading. (a) thin wall tube with fiber orientation. (b) typical experimental creep data under axial and shear loading.

7

. . . _ . log(ε /ε ΤΝ ) or log(γ / V 3 ε TN )

1

0

-1 90-ten 60-ten 45-ten 90-shear

-2

Correlations

(LN) - prediction (TS) - prediction

-3 -0.4

-0.2

0.0

0.2

0.4

log(σ/σ0) or log( 3 τ/σ0)

Fig 2.3 Non-dimensional log creep rate vs. log stress for the complete exploratory data set. Tensile and shear data are shown as indicated. Fiber angles θ = 60 o and 45 o and shear data for θ = 90 o shifted to form a master curve.

8

CHAPTER III THE ISOTROPIC POWER LAW MODEL

The power-law model can be used to do the creep calculation for isotropic material. The constitutive equation of isotropic power law is as following:

ε& cr = Aσ t m

(3.1)

where

ε& cr

is the uniaxial equivalent creep strain rate,

σ

is the uniaxial equivalent deviatoric stress,

t

is the total time, and

A, n, and m are defined by the user as functions of temperature. σ is Mises equivalent stress or Hill's anisotropic equivalent deviatoric stress according to whether isotropic or anisotropic creep behavior is defined (discussed below). For physically reasonable behavior

and n must be positive and

. Since total time is used in the

expression, such reasonable behavior also typically requires that small step times compared to the creep time be used for any steps for which creep is not active in an analysis; this is necessary to avoid changes in hardening behavior in subsequent steps. Based on Robinson creep model we can derive A, n, m. The Robinson model to calculate creep rate in loading direction is:

9

n +1 σ n ε& = ε&0 ( ) (1 − ξ ) 2 σ0

(3.2)

Because it is for isotropy we define

ξ = ζ =η = 0 Now we obtained

ε& = 3.11× 10−15 σ 6.5t 0 A = 3.11× 10−15 , n = 6.5 , m = 0 Depending on the choice of units for either form of the power law, the value of A may be very small for typical creep strain rates. If A is less than 10 −17 , numerical difficulties can cause errors in the material calculations; therefore, use another system of units to avoid such difficulties in the calculation of creep strain increments.

10

CHAPTER IV USER SUBROUTINE UMAT

The material law is implemented computationally as a user defined subroutine (UMAT) in a commercially available FEA code (ABAQUS). User subroutines provide an extremely powerful and flexible tool for analysis. This chapter defines the interfaces for the user subroutines that are available in ABAQUS.

4.1 Introduction of User subroutine (UMAT)

User subroutine UMAT can be used to define the mechanical constitutive behavior of a material; it must update the stresses and solution-dependent state variables to their values at the end of the increment for which it is called it must provide the material Jacobian matrix, for the mechanical constitutive model; it can be used in conjunction with user subroutine USDFLD to redefine any field variables before they are passed in. It is sometimes desirable to set up the FORTRAN environment and manage interactions with external data files that are used in conjunction with user subroutines. For example, there may be history-dependent quantities to be computed externally, once per increment, for use during the analysis; or output quantities that are accumulated over multiple elements in COMMON block variables within user subroutines may need 11

to be written to external files at the end of a converged increment for postprocessing. Such operations can be performed with user subroutine UEXTERNALDB. This user interface can potentially be used to exchange data with another code, allowing for “stagger” between ABAQUS and another code. User subroutines should be written with great care. To ensure their successful implementation, the rules and guidelines below should be followed. Every user subroutine must include the statement. As the first statement after the argument list. The file ABA_PARAM.INC is installed on the system by the ABAQUS installation procedure. It specifies either IMPLICIT REAL*8 (A-H, O-Z) for double precision machines or IMPLICIT REAL (A-H,O-Z) for single precision machines. The ABAQUS execution procedure, which compiles and links the user subroutine with the rest of ABAQUS, will include the ABA_PARAM.INC file automatically. It is not necessary to find this file and copy it to any particular directory; ABAQUS will know where to find it. 1) User subroutines must perform their intended function without overwriting other parts of ABAQUS. In particular, the user should redefine only those variables identified in this chapter as “variables to be defined.” Redefining “variables passed in for information” will have unpredictable effects. 2) When developing user subroutines, test them thoroughly on smaller examples in which the user subroutine is the only complicated aspect of the model before attempting to use them in production analysis work. If needed, debug output can be written to FORTRAN unit 7 to appear in the message (.msg) file or to FORTRAN unit 6 to appear in the data (.dat) file; these units should not be 12

opened by the user's routines since they are already opened by ABAQUS. FORTRAN units 15 through 18 or units greater than 100 can be used to read or write other user-specified information. The use of other FORTRAN units may interfere with ABAQUS file operations. These FORTRAN units must be opened by the user; and because of the use of scratch directories, the full pathname for the file must be used in the OPEN statement. 3) Solution-dependent state variables are values that can be defined to evolve with the solution of an analysis.

4.2 Implementation of UMAT

The equation of motion together with the constitutive law form system consisting of an initial – boundary problem and an ordinary differential equation. The equation of motion is solved with the help of a finite-element package (ABAQUS), and the constitutive law by a solver for ordinary differential equations. The relevant constitutive information is passed to ABAQUS by a subroutine UMAT which has to be supplied by the user. Starting from an equilibrium at time t n , ABAQUS performs an (incremental) loading as well as with the time increment ∆t , and an initial guess ∆ε n , for the strain increment. The user subroutine UMAT has to supply ABAQUS with new Cauchy stress tensor σ (t n + ∆t ), updated according to the constitutive law as well as with the derivative of stress with respect to the strain increment. With this information, a new guess for the strain increment is calculated and the whole procedure is iterated until convergence. The precise information on the Jacobian is essential to achieve fast 13

convergence in Newton-type iteration performed by ABAQUS. The following is the detail introduction of definition of Jacobian Matrix in UMAT and the Newton Ralphson method incorporated in ABAQUS. The following diagram (Fig 4.3) indicates a certain step of UMAT working with ABAQUS. STEP 1: At an equilibrium time t n , ABAQUS will supply time increment ∆t , and

total strain increment ∆ε total (t n ) and total strain ε total (t n ) for the strain increment. And also

σ (t n ) is calculated by previous increment. All these four values will pass to UMAT to calculate new Cauchy stress tensor σ (t n + ∆t ) . STEP 2: Stress Update. According to Robinson creep damage Model, creep strain increment 3 Γ ∆ε = Φ n −1 ∆t 2 σ0 Also, the relationship between total strain, creep strain and elastic strain is ∆ε total (t n ) = ∆ε e (t n ) + ∆ε (t n ) so

∆ε e (t n ) = ∆ε total (t n ) − ∆ε (t n ) In this way we can update Cauchy Stress tensor. The stress increment is ∆σ (t n ) = J ⋅ ∆ε e (t n ) The new Cauchy Stress tensor is

σ (t n + ∆t ) = σ (t n ) + ∆σ (t n ) 14

STEP 3: Strain Update ABAQUS will update strain tensor

ε total (t n + ∆t ) = ε total (t n ) + ∆ε total (t n ) STEP 4: ABAQUS will generate new total strain increment ∆ε total (t n +1 ) STEP 5: ABAQUS equilibrium iterations at new time t n+1 , the maximum iteration number is set to 9 in ABAQUS and error tolerance is set to TOL 5e-3. If less than 9 times iterations the error is less than TOL, it calls convergence, n=n+1 and move to another increment. If the error is larger than TOL, ABAQUS will reduce the ∆t go to step 1 until convegnece. Fig 4.1 briefly illustrates one increment in UMAT calculation.

15

step1 step n ABAQUS supply information at time tn Time increment ∆t total Total strain increment ∆ε (tn ) ε ( t ) Total strain n

∆t

ε ntotal σ n

∆ε total (tn )

UMAT

step2

Γ 3 ∆ε ij = Φ n−1 ij ∆t 2 σo

Dr. Robinson’s creep law

∆ε total (tn ) = ∆ε e (tn ) + ∆ε (tn )

∆εije (tn ) = ∆εijtotal(tn ) − ∆εij (tn ) Stress update

∆σ (tn ) = J ⋅ ∆ε e (tn )

σ n+1 = σ (tn + ∆t ) = σ (tn ) + ∆σ (tn ) step3 Strain update

ABAQUS update

ε n+1 = ε total

total

(tn + ∆t ) = ε total (tn ) + ∆ε total (tn )

step4 ABAQUS generate new ∆ε total (tn+1 )

step5 ABAQUS equilibrium iterations at time tn+1

convergence

Not convergence Reduce ∆t

n=n+1, go to step 1

Fig 4.1

UMAT Loop 16

4.2.1 Jacobian Matrix for plane stress element and shell element

A class of common engineering problems involving stresses in a thin plate or on the free surface of a structural element, such as the surfaces of thin-walled pressure vessels under external or internal pressure, the free surfaces of shafts in torsion and beams under transverse load has one principal stress that is much smaller than the other two. By assuming that this small principal stress is zero, the three-dimensional stress state can be reduced to two dimensions. Since the remaining two principal stresses lie in a plane, these simplified 2D problems are called plane stress problems.

Fig 4.2 Plane Stress element

Assume that the negligible principal stress is oriented in the z-direction. To reduce the 3D stress matrix to the 2D plane stress matrix, remove all components with z subscripts to get,

17

where τ xy

⎡σ x ⎢τ ⎣ yx

τ xy ⎤ σ y ⎥⎦

= τ yx for static equilibrium. The sign convention for positive stress

components in plane stress is illustrated in the above Fig 4.1 on the 2D element.

Jacobian matrix of the constitutive model for plane stress and shell element,

∂∆σ / ∂∆ε , where ∆σ are the stress increments and ∆ ε are the strain increments can be defined in following steps.

STEP 1 .Definition of Transformation Matrix [T]

Fig 4.3 Plane Stress Transformation from local coordinate (x-y) to global coordinate (12)

18

According to the Fig 4.3 we can calculate the resultant force in 1-2 coordinate: ∑ F1 = σ 1 A − σ y sin θA cos θ − τ s cos θA sin θ − σ x cos θA cos θ − τ s sin θA cos θ = 0

(4.1)

∑ F2 = τ 12 A − σ y cos θA sin θ + τ s sin θA cosθ + σ x sin θA cos θ − τ s cos θA cos θ = 0

(4.2)

We can obtain stress components in 1-2 coordinate: σ 1 = m 2σ x + n 2σ y + 2mnτ s

(4.3)

τ 6 = τ 12 = − mnσ x + mnσ y + (m 2 − n 2 )τ s

(4.4)

where

m = cos θ

n = sin θ

(4.5)

We conclude that: 2 n2 2mn ⎤ ⎡σ x ⎤ ⎡σ 1 ⎤ ⎡ m 2 ⎢σ ⎥ = ⎢ n 2 m − 2mn ⎥ ⎢σ y ⎥ ⎥⎢ ⎥ ⎢ 2⎥ ⎢ ⎢⎣τ 6 ⎥⎦ ⎢⎣− mn mn ( m 2 − n 2 )⎥⎦ ⎢⎣ τ s ⎥⎦

⎤ ⎡ ⎤ ⎡ ⎢ εx ⎥ ⎢ ε1 ⎥ ⎢ ε ⎥ = [T ]⎢ ε ⎥ ⎢ y ⎥ ⎢1 2 ⎥ ⎢1 γ ⎥ ⎢ γ6⎥ ⎢⎣ 2 s ⎥⎦ ⎣2 ⎦

→ [σ ]1, 2 = [T ][σ ]x , y

→ [ε ]1, 2 = [T ][ε ]x , y

(4.6)

(4.7)

STEP 2. Transformation of the reduced stiffness matrix According to following procedure we can calculate the transformed reduced stiffness matrix

19

⎡σ 1 ⎤ ⎡Q11 ⎢σ ⎥ = ⎢Q ⎢ 2 ⎥ ⎢ 12 ⎢⎣τ 6 ⎥⎦ ⎢⎣ 0

⎡σ x ⎤ ⎢σ ⎥ = T −1 ⎢ y⎥ ⎢⎣ τ s ⎥⎦

⎡ ⎤ 0 ⎤ ⎢ ε1 ⎥ ⎡σ 1 ⎤ ⎡Q11 ⎥ ⎢ ⎥ 0 ⎥ ε 2 ⇒ [T ]⎢⎢σ 2 ⎥⎥ = ⎢⎢Q12 ⎢ ⎥ ⎢⎣ τ s ⎥⎦ ⎢⎣ 0 2Q66 ⎥⎦ ⎢ 1 γ ⎥ 6 ⎣⎢ 2 ⎦⎥

Q12 Q22 0

⎡Q11 ⎢Q ⎢ 12 ⎢⎣ 0

[ ]

Q12 Q22 0

⎤ ⎡ 0 ⎤ ⎢ εx ⎥ 0 ⎥⎥[T ]⎢ ε y ⎥ ⇒ ⎥ ⎢ 2Q66 ⎥⎦ ⎢ 1 γ ⎥ ⎢⎣ 2 s ⎥⎦

⎡σ x ⎤ ⎡Qxx ⎢σ ⎥ = ⎢Q ⎢ y ⎥ ⎢ yx ⎢⎣ τ s ⎥⎦ ⎢⎣Qsx

Q12 Q22 0

Qxy Qyy Qsy

⎤ ⎡ 0 ⎤ ⎢ ε1 ⎥ 0 ⎥⎥[T ]⎢ ε 2 ⎥ ⇒ ⎥ ⎢ 2Q66 ⎥⎦ ⎢ 1 γ ⎥ s ⎣⎢ 2 ⎦⎥

⎡ ⎤ 2Qxs ⎤ ⎢ ε x ⎥ ⎥ 2Qys ⎥ ⎢ ε y ⎥ ⎢ ⎥ 2Qss ⎥⎦ ⎢ 1 γ ⎥ ⎢⎣ 2 s ⎥⎦

In this way we can obtain each component in Transformed Stiffness Matrix Qxx = m 4 Q11 + n 4 Q22 + 2m 2 n 2 Q12 + 4m 2 n 2 Q66

(4.8)

Eq. (4.8) is the formula to calculate the components in Jacobian Matrix for plane stress element with different orientation.

4.2.2 Newton-Ralphson Method in ABAQUS

Newton's method, also called the Newton-Ralphson method, is a root-finding algorithm that uses the first few terms of the Taylor series of a function f(x) in the vicinity of a suspected root. Newton's method is sometimes also known as Newton's iteration, although in this work the latter term is reserved to the application of Newton's method for computing square roots. For f(x) a polynomial, Newton's method is essentially the same as Horner's method. The Taylor series of f(x) about the point x = x0 + ε is given by 20

f ( x0 + ε ) = f ( x0 ) + f ' ( x0 )ε +

1 f " ( x0 )ε 2 + ... 2

(4.9)

Keeping terms only to first order, f ( x0 + ε ) = f ( x0 ) + f ' ( x0 )ε

(4.10)

This expression can be used to estimate the amount of offset

ε

needed to land closer to

the root starting from an initial guess x0 . Setting f ( x0 + ε ) = 0 and solving (4.10) for

ε = ε 0 gives ε0 = −

f (x0 ) f ' (x0 )

(4.11)

which is the first-order adjustment to the root's position. By letting x1 = x0 + ε 0 , calculating a new ε 1 , and so on, the process can be repeated until it converges to a root using ε0 = −

f ( xn ) f ' ( xn )

(4.12)

Unfortunately, this procedure can be unstable near a horizontal asymptote or a local extremum. However, with a good initial choice of the root's position, the algorithm can by applied iteratively to obtain xn+1 = xn −

f ( xn ) f ' ( xn )

(4.13)

For n=1, 2, 3, .... An initial point

that provides safe convergence of Newton's method

is called an approximate zero. The error ε n+1 after the (n+1) st iteration is given by ε n+1

=

ε n + ( xn+1 − xn )

=

εn −

(4.14)

f ( xn ) f ' ( xn )

(4.15) 21

But =

f ( x0 )

f ( xn−1 ) + f ' ( xn−1 )ε n +

=

f ' ( xn−1 )ε n +

1 f " ( xn−1 )ε n2 + ... 2

f ' ( xn−1 ) + f " ( x)ε n + ...

=

f ' ( xn )

1 f " ( xn−1 )ε n2 + ... 2

(4.16) (4.17) (4.18)

so f ( xn ) f ' ( xn )

1 f " ( xn−1 )ε n2 + ... 2 f ' ( xn−1 ) + f " ( xn−1 )ε n + ...

f ' ( xn−1 )ε n +

=

1 f " ( xn−1 )ε n2 . 2 f ' ( xn−1 )

f ' ( xn−1 )ε n +

≈ =

εn +

f " ( xn−1 ) 2 εn 2 f ' ( xn−1 )

(4.19)

(4.20)

(4.21)

and (4.14) becomes ε n+1

=

ε n − [ε n +

f " ( xn −1 ) 2 εn ] 2 f ' ( xn−1 )

(4.22)

ABAQUS sets the error tolerance TOL 5e-3 and maximum iteration times 9. If after less than 9 times iteration

ε n +1 < TOL

(4.23)

we can call convergence.

22

CHAPTER V PLANE STRESS ELEMENT IMPLEMENTATION

In the previous chapter we introduce the theory of user subroutine UMAT and how to implement user material law through UMAT. When developing user subroutines, test them thoroughly on smaller examples in which the user subroutine is the only complicated aspect of the model before attempting to use them in production analysis work. First we do single element test.

5.1 Single element test

We first investigate the one-element tension test, given in Fig 5.1. The dimension of this element is 10mm x10mm.

23

Y

Y

u

σ

θ

θ

(a)

(b)

Fig 5.1 One element model problem set up with fiber orientation. (a) creep loading and (b) constant displacement loading.

Fig 5.1 shows the schematic with boundary conditions for a single element model. This particular test was done for the reduced integration plane stress elements with linear and quadratic interpolation schemes ( ABAQUS element type CPS4R and CPS8R). The elements were tested for response under constant load (creep) as well as constant displacement (relaxation) boundary conditions. The material parameters used

24

for this purpose are the ones obtained in the thin cylinder experiment described above and are as follows:

α = 0.35 , ν = 10.6 , m = 7.125 , n = 6.5 , σ 0 = 6671.74 psi , ξ = 0.57 , η = 0.1 , ζ = 0.64 , ε&0 = 0.01% perhour ,

t 0 = 12.0hours This particular set of parameters assumes that the isochronous damage function described in equation (3.8) is of the order 1 (P =1). Fig 5.2 shows comparison between power law and Robinson creep model for isotropic case under a constant load of 45 MPa when the scalar damage variable is not included in the Robinson creep model. The Fig 5.2.(a) indicates that at creep strain evolution after 10 hours at 1 direction calculated by power law model and Robinson creep model are agree with each other, the same to 2 direction (Fig. 5.2.(b)). But the creep strains evolution after 10 hours for 1-2 direction are different for each material model (Fig.5.2.(c)) Because Robinson Model is anisotropic constitutive equation while Power Law is isotropic constitutive equation. The Power Law equation to calculate creep strain rate:

ε&11 = ε&22 = γ&12 = ε&0 (

n +1 σ n ) (1 − ξ ) 2 σ0

The Robinson Model to calculate creep strain rate:

ε&11 = ε&0 (

σ n ) (1 − ξ ) σ0

n +1 2

n +1 σ n ε&22 = ε&0 ( ) (1 − η ) 2 σ0

25

γ&12 / 3 = ε&0 (

3τ

σ0

) n (1 − ξ )

n +1 2

According to isotropy we define

ξ = ζ =η = 0 In this case the Power Law equation reduced to

ε&11 = ε&22 = γ& = ε&0 (

σ n ) σ0

Robinson Model reduced to

ε&11 = ε&22 = ε&0 (

γ&12 / 3 = ε&0 (

σ n ) σ0 3τ

σ0

)n

Fig 5.3 shows the response of a single element creep loading test under a constant load of 46 MPa with fiber orientation 90 deg, the scalar damage variable is not included. After compare FEA result with experimental date we observed that both experiment and FEA calculation is with creep rate relatively constant at 0.013%/hr. We also compare UMAT results with experiment data under different tensile stress with fiber orientation 45 and 90 deg. The UMAT results well coincide with experiment data (Table 1 and Table 2). Fig 5.4 shows the response of a single element creep loading test. The Fig 5.4.a shows the time evolution of creep strain under a constant load of 45 MPa for different fiber orientations (0 deg, 45 deg, 90 deg) when the scalar damage variable is not included in the material model. It is observed that when the creep loading is along fiber directions, the PMC has the strongest behavior. The creep strains are progressively 26

higher when the creep loading is at progressively higher angle to the fiber orientation. For the 90 deg orientation, failure can be achieved in roughly 10 hours. Similarly, for a given orientation, the creep strains increase progressively, as the creep load is increased (Fig. 5.4.(b)). Fig 5.5.(a) shows the time evolution of creep with fiber orientations 45 deg under loads equal to 60MPa when the scalar damage variable is included. We find that when include this damage factor the plot is getting curved and the value of the creep strain is higher than that without this damage factor. Fig 5.5.(b) shows the time evolution of creep with fiber orientations 45 deg under different loads (when the scalar da70Mpa, 75Mpa, 80Mpa) with damage variable is included. For creep load less than 70 MPa, the evolution of creep strain with time is linear. But above a critical load, the increased damage results in a larger creep strain. We also plot the time evolution of shear strain with different fiber orientation 0 deg, 45 deg and 90 deg under shear loading 40 MPa (Fig 5.5.(c)). It can be observed that 0 deg and 90 deg have the same strain evolution behavior and 45 deg has more strain evolution than 0 deg and 90 deg.

27

Table 1 fiber orientation 90 deg under tension

σ /σ 0

Experiment data

0.94 1 1.04

0.64 1 1.46

ε& / ε&TN

UMAT Results

ε& / ε&TN

0.67 1 1.34

Table 2 fiber orientation 45 deg under tension

σ /σ 0

Experiment

0.54 0.6 0.65

0.005 0.0092 0.012

data

ε& / ε&TN

UMAT Results Experiment

ε& / ε&TN

data

0.0044 0.0094 0.014

-0.0034

28

γ& / ε&TN

UMAT Results

γ& / ε&TN

-0.0021

.0025 Robinson model power law

.0020

Robinson model power law

0.0000

-.0002 .0015

ε11

ε11

-.0004 .0010

-.0006 .0005

-.0008

0.0000

-.0010

-.0012 0

2

4

6

8

10

0

12

2

4

6

8

10

12

Time(hr)

Time(hr)

(a)

(b)

4e-17 Robinson model power law

3e-17

ε11

2e-17

1e-17

0

0

2

4

6

8

10

12

Time(hr)

(c)

Fig 5.2 Comparison between power law and Robinson creep model for isotropic case under a constant load of 45 MPa when the scalar damage variable is not included in the Robinson creep model. (a) Time evolution of creep strain in Y direction. (b) Time evolution of creep strain in Y direction. (c) Time evolution of creep strain in XY direction

29

0.10 UMAT Experiment

axial strain (%)

0.08

0.06

0.04

0.02

0.00 0

2

4

6

8

time (hr)

Fig 5.3 Time evolution of creep strain in the direction of loading (without damage) for fiber orientation of 90 deg at 45MPa load.

30

.0010 00 452 900

.0008

ε11

.0006

.0004

.0002

0.0000

0

1

2

3

4

5

6

4

5

6

Time(hr)

(a)

.0010 20Mpa 46Mpa 60Mpa

.0008

ε11

.0006

.0004

.0002

0.0000

0

1

2

3

Time(hr)

(b)

Fig5.4 Time evolution of creep strain in the direction of loading (without damage). (a) for different orientations of the fiber (0, 45, 90 deg) at 45MPa load. (b) for 45 deg fiber orientation at varying loading 20, 45 and 100PMa.

31

.0020 creep with damage creep without damage .0015

ε11

.0010

.0005

0.0000

0

2

4

6

8

10

12

Time(hr)

(a) .005 70Mpa 75MPa 80MPa

.004

ε11

.003

.002

.001

0.000

0.0

.5

1.0

1.5

2.0

2.5

Time(hr)

(b) 6e-4 0 deg 45 deg 90 deg

5e-4

4e-4

ε

12

3e-4

2e-4

1e-4

0

0

1

2

3

4

5

6

Time (hr)

(c) Fig 5.5 Time evolution of creep strain with damage effect. (a) comparison when damage evolution is included under tensile load 60 MPa. (b) comparison between tensile loads 70, 75, 80 MPa with fiber orientation 45 deg. (c) comparison between fiber orientation 0, 45, 90 deg under shear loads 40 Mpa. 32

We also apply relaxation test on single element. We apply fixed displacement at the right boundary (u=0.1%). The dimension is still 10mm x 10mm. We compare stress calculated by Robinson’s creep law with that by power law in isotropic case. We find out that after 10 hours test the stress along loading direction calculated by Robinson Creep model is lower than the result calculated by Power law mode, which is due to the damage effect of Robinson creep model.

Power Law

Robinson Model

Fig 5.6: Time evolution of stress relaxation in the direction of loading (with damage).

33

Results for both the creep and relaxation tests were able to match the results from an isotropic Power Law Creep Model when damage was not included. However, the proposed model deviates from the Power Law model when creep anisotropy due to orientation and damage is included. These element tests are achieved for a wide range of loads and orientations, demonstrating the robustness of the numerical scheme. Also, the tests show the usefulness of the constitutive model despite the use of a single scalar parameter for damage variable.

5.2 Square plate

One element tests are demonstrating the robustness of the numerical scheme of the UMAT. In this part we extend our model to relatively complex model – a square plate with multiple elements. The dimension of the plate is 20mm x 20 mm. This particular test was done for the reduced integration plane stress element with linear and quadratic interpolation schemes (ABAQUS element type CPS4R and CPS8R). Geometry and boundary condition of the plate under tension is shown on Fig 5.7. We apply tensile stress 46 Mpa stress along 1 direction with different fiber orientation (0 deg, 45 deg and 90 deg). Fig 5.8 shows the time evolution of creep strain with fiber orientations 0, 45 and 90 deg under tensile loads equal to 46MPa. Fig 5.8.(a) shows the 1 direction strain distribution and Fig 5.8.(b) shows the shear strain distribution. All these results coincide with the results obtained in one element tests.

34

σ

σ

20mm

Fig 5.7 Geometry and boundary condition of square plate problem.

35

.0008 00 450 0 90

.0006

ε11

.0004

.0002

0.0000

0

1

2

3

4

5

6

4

5

6

Time (hr)

(a)

4e-5 00 450 900

3e-5

ε12

2e-5

1e-5

0

-1e-5 0

1

2

3

Time (hr)

(b)

Fig5.8 Time evolution of creep strain with different orientations of the fiber (0, 45, 90 deg) at 46MPa load in the direction of loading (without damage). (a) time evolution of strain at 1 direction. (b) time evolution of shear strain.

36

5.3 Plate with a hole in the middle

This particular test was done for the reduced integration plane stress element with linear and quadratic interpolation schemes (ABAQUS element type CPS4R and CPS8R). Geometry and boundary condition of the plate under tension and torsion is shown on Fig 5.8. We apply tensile stress 10 MPa along 1 direction with different fiber orientation (0 deg, 45 deg and 90 deg).

L=152.4mm

σ R=80mm

Fig 5.9 Left, geometry and Right, 2D quarter symmetry model for the square plate with a hole problem.

37

First of the study we choose isotropic material to compare Robinson creep Model with Power law creep model. The Power law model is one of the creep model ABAQUS use to do creep calculation for isotropic material. Fig 5.10 shows the stress distribution of the Robinsons’ model law for isotropic material and comparative stress distribution for a Power Law type model. It can be seen that stress distribution after instantaneous elastic response for the proposed Creep Damage model (Fig 5.10.(b)) is identical to a Power Law model (Fig 5.10.(a)). This is because in elastic step both power law and Robinson model has the same constitutive equation. But due to the damage effect we take into consideration in Robinson model we can observe some difference existing after 5 hours creep response in Fig 5.10.(c) and Fig 5.10.(d). To further compare these two models we plot stress distribution along the hole (counter clockwise) because this area has the most complex behavior of the whole plate. Fig 5.11 shows the distribution of Stress distribution along the hole, after elastic step (a) and after 5 hours of creep load (b). Fig 5.119.(a) shows the plots along the hole for both creep law are identical to each other after elastic step. Fig 5.119.(b) shows lower stresses in the Robinson Model at the second half of the plot is a result of higher stress relaxation from higher stress concentration area due to the damage evolution in Robinson Model. Fig 5.12 shows the creep strain distribution for Power law and Robinson creep Model after 5 hours creep response. Robinson Model has more creep strain existing in the high stress concentration zone. Fig 5.13 shows the creep strain distribution along the hole (counter clockwise), after 5 hours of creep load. We can observe that at the first part of the path Robinson model is agree with Power law model. At the second half of the path Robinson model undergoes more creep strain than power 38

law that is due to the damage factor we include in the calculation of creep strain in Robinson Model. We also find the highest value located in different area in these two models. The highest value for Power law located in the very end of the path while for Robinson model highest value of stress redistribution is in the area close to the end of the path which is due to the anisotropic constitutive equations of Robinson Model even for isotropic material. To verify the accuracy of the numerical scheme we compare the numerical results with analytical solutions. We choose the radius of the plate 20mm and 50mm with fiber orientation 90 deg (Fig 5.14). We compare each case with analytical solution. Both two plots indicated the accuracy of the numerical result in positive value area. We can find that there is big difference in negative area. Analytic solution has much lower value than numerical solution in negative zone. The discrepancy is due to the reason that analytic solution is for infinite plate and free boundary calculation. Fig 5.15 shows the stress distribution along the hole (counter clockwise), after 5 hours of creep load 10 Mpa with R=20mm, 50mm, 80mm, L=152.4mm and fiber orientation 45 deg. Fig 5.16 shows the stress distribution along the hole (counter clockwise), after 5 hours of creep load 10 Mpa with R=50mm, L=152.4mm.Fiber orientation 45deg and 90 deg. Fig 5.17 shows Stress concentration zone with damage effect under axial stress 10Mpa fiber orientation 45 deg (a)R/L=0.1 (b)R/L=0.3 (c) R/L=0.5. Fig 5.18 shows Stress compression zone with damage effect under axial stress 10Mpa fiber orientation 45 deg (a) R/L=0.1 (b) R/L=0.3 (c) R/L=0.5. All these results indicates that the computational routine (UMAT) successfully describes the rather complex creep/damage phenomenon observed in PMCs. 39

78MPa

78MPa

50MPa

50MPa

38MP

38MPa

20MPa

20MPa

(a)

(b)

40MPa

38MPa 20MPa

20MPa

(c)

(d)

Fig 5.10 Contour plot of (a) stress distribution of the power law for isotropic material after elastic step, (b) stress distribution of the Robinsons’ model law for isotropic material after elastic step, (c) stress distribution of the power law for isotropic material after creep step (5 hours), (d) stress distribution of the Robinsons’ model law for isotropic material after creep step (5 hours).

40

100 Robinson model powerlaw

80

σ11(MPa)

60

40

20

0

-20 0

20

40

60

80

100

120

140

True distance along the hole (counter clockwise)

(a)

50 Robinson model power law

40

σ11(MPa)

30

20

10

0

-10 0

20

40

60

80

100

120

140

True distance along the hole (counter clockwise)

(b)

Fig 5.11 Stress distribution along the hole (quarter circle) going in a counter clockwise direction for the proposed Creep Damage Model as and for isotropic Power Law Creep model. (a) stress distribution after elastic step. (b) stress distribution after 5 hours.

41

0.003% 0.014% 0.033% 0.051% 0.073%

(a)

0.003% 0.014% 0.033% 0.051% 0.073%

(b) Fig 5.12 FEA plot of Creep strain distribution after 5 hours creep response. (a) Robinson Damage Model and (b) Isotropic Power Law.

42

.0012 Robinson model power law

.0010

.0008

ε11

.0006

.0004

.0002

0.0000

-.0002 0

20

40

60

80

100

120

140

True distance along the hole (counter clockwise)

Fig 5.13 Creep strain along the hole (quarter circle) going in a counter clockwise direction for the Robinson Damage Model and for isotropic Power Law Creep model (after 5 hours).

43

40

30

Robinson model Analytic solution

σ11(MPa)

20

10

0 θ -10

-20

-30 0

20

40

60

80

100

60

80

100

θ

(a) 40

30

Robinson model Analytic solution

σ11(MPa)

20

10

0

-10

-20

-30 0

20

40

θ

(b)

Fig 5.14 Comparison between Analytic solution and Robinson model (FEA) in elastic step. (a) R=20mm and R/L=13%. (b) R=50mm and R/L=33%.

44

60 R=80 R=50 R=20

50

σ11(ΜPa)

40

30

20

10

0

-10 0

20

40

60

80

100

120

140

True distance along the hole (counter clockwise)

Fig 5.15 Creep / time response along the hole (quarter circle) going in a counter clockwise with damage effect under axial stress 10MPa with R=20mm,50mm,80mm, L=152.4mm, fiber orientate 45 deg.

45

40 450 900

σ11(MPa)

30

20

10

0

-10 0

20

40

60

80

100

θ (deg)

Fig 5.16 Creep / time response along the hole (quarter circle) going in a counter clockwise with damage effect under axial stress 10MPa with R=50mm, L=152.4mm, fiber orientate 45 deg and 90 deg.

46

20MPa 25 MPa

(a)

20MPa 40Mpa

20MPa

(b)

25 MPa 30Mpa 40Mpa 45Mpa

56MPa

(c) Fig 5.17 Stress concentration zone with damage effect under axial stress 10Mpa fiber orientation 45 deg. (a) stress concentration zone with R/L=0.1. (b) stress concentration zone with R/L=0.3. (c) stress concentration zone with R/L=0.5. 47

-0.5MPa

(a)

-2Mpa

-0.5MPa

-6MPa

(b) -4Mpa -3Mpa

-2Mpa

-0.5MPa

(c) Fig 5.18 Stress compression zone with damage effect under axial stress 10Mpa fiber orientation 45 deg. (a) stress compression zone with R/L=0.1. (b) stress compression zone with R/L=0.3. (c) stress compression zone with R/L=0.5. 48

CHAPTER VI SHELL ELEMENT IMPLEMENTATION

6.1 theory of shell element

Shell elements are surface representations of structures that are much thinner in one direction than the other two (thin-walled structures). The elements are geometrically defined by three or four sided surfaces, and are located in space at the mid-plane of the solid they are representing. The user specifies the thickness of the elements as an input to the software. These elements are used in modeling all types of thin-walled structures, such as airplane and automotive bodies, pressure vessels, sheet metal, and many plastic molded parts. Shell elements have six active degrees of freedom per node, much like beam elements. Because of commonality in degrees of freedom, beam and shell elements are often joined together in mixed models.

49

Shell elements are good for modeling structures that are thin. These elements are usually formulated under the assumptions governing thin plate theory. If a structure is too thick, the behavior of thin plates is no longer seen (shear stresses become large, etc.), and shell elements should not be used. This limit is usually seen at a thickness to width or length (whichever is smaller) of 1/10 and larger. Shell elements generally have a lower limit on this ratio as well. At thickness-towidth ratios between 1/100 and 1/1000, thin plates begin to behave like membranes, with no bending stiffness (like a string in tension, subjected to transverse load). Because of this shell elements cannot be used to model very thin, flexible structures such as fabric or thin membranes.

6.2 Shell element in ABAQUS

ABAQUS includes general-purpose shell elements as well as elements that are valid for thick and thin shell problems. See below for a discussion of what constitutes a “thick” or “thin” shell problem. This concept is relevant only for elements with displacement degrees of freedom. The general-purpose shell elements provide robust and accurate solutions to most applications and will be used for most applications. However, in certain cases, for specific applications, enhanced performance may be obtained with the thin or thick shell elements; for example, if only small strains occur and five degrees of freedom per node are desired. Element type S4R is general-purpose shell. These elements allow transverse shear deformation. They use thick shell theory as the shell thickness increases and become 50

discrete Kirchhoff thin shell elements as the thickness decreases; the transverse shear deformation becomes very small as the shell thickness decreases. Element type S8R should be used only in thick shell problems. Thick shells are needed in cases where transverse shear flexibility is important and second-order interpolation is desired. When a shell is made of the same material throughout its thickness, this occurs when the thickness is more than about 1/15 of a characteristic length on the surface of the shell, such as the distance between supports for a static case or the wavelength of a significant natural mode in dynamic analysis.

6.3 Implementation of shell element

First, the computational model was applied to the thin-walled tube (7.5 in dia) under a tensile stress of 0.5 MPa. This particular test was done for the reduced integration shell elements with linear and quadratic interpolation schemes (ABAQUS element type S4R and S8R). Geometry and boundary condition of thin wall tube under tension and torsion is shown on Fig 6.1. We apply tensile stress along 1 direction with different fiber orientation (0 deg, 45 deg and 90 deg). Because the stress and strain distributions are evenly along the tube, we take one element to analyze the strain development. Fig 6.2 shows the time evolution of creep strain under a constant tensile load of 45 MPa for different fiber orientations (0 deg, 45 deg, 90 deg) when the scalar damage variable is included in the material model. It is observed that when the creep loading is along fiber directions, the PMC has the strongest behavior. The creep strains are progressively higher when the creep loading is at progressively higher angle to the 51

fiber orientation. For the 90 deg orientation, failure can be achieved in roughly 6 hours. Fig 6.3 shows the time evolution of creep strain under a constant shear load of 45 MPa for different fiber orientations (0 deg, 45 deg, 90 deg) when the scalar damage variable is included in the material model.

52

2

T F

F θ 1 T

3

Fig 6.1 Thin wall tube under tension and torsion with fiber orientation

53

(a)

.0052 0 deg 45 deg 90 deg

.0050 .0048 .0046

ε12

.0044 .0042 .0040 .0038 .0036 .0034 .0032 0

2

4

6

8

10

12

Time (hr)

(b)

Fig 6.2 Time evolution of creep strain under a constant tensile load (F) 45 MPa of thinwalled tube for different fiber orientations (0 deg, 45 deg, 90 deg). (a) time evolution of maximum principle strain. (b) time evolution of shear strain. 54

Further the computational model was applied to the pressure vessel (20 cm dia) under a pressure of 0.5 MPa. This particular test was done for the reduced integration shell elements with linear and quadratic interpolation schemes (ABAQUS element type S4R and S8R). Fig 6.3 shows the geometry, mesh and path of the pressure vessel, the thickness is 0.4. First we plot the path along path 1 (Fig 6.3) with different fiber orientation. 0 deg is the direction along the path 1, 90 deg is the direction perpendicular to the path1. It can be observed that 0 deg has the smallest strain deformation along the path1, the higher angle of fiber, the higher strain deformation occurs. We can also find the movement of the strain distribution along path 1 with fiber orientation of 45 deg (Fig. 6.4) after 10 hours evolution. This phenomena is also obvious in other fiber orientation. Fig 6.4 shows path plot of time evolution of Maximum Strain along path1 of the vessel with fiber orientation 45 deg under inside pressure 0.5Mpa with damage evolution. Fig 6.5 shows Maximum Strain along path1 of the vessel with fiber orientation 0, 45, 90 deg under inside pressure 0.5Mpa with damage evolution after 10 hours. Fig 6.6 shows path plot of time evolution of Maximum Strain along path2 of the vessel with fiber orientation 45, 60, 90 deg under inside pressure 0.5Mpa with damage evolution. It can be observed that the obvious difference in maximum principle strain between 45, 60 and 90 deg is at the center part of the pressure vessel and the 90 deg has the highest strain development. Fig 6.7 shows the time evolution of Maximum Principle Strain of the pressure vessel with fiber orientation 60, 90 deg under pressure 1Mpa with damage evolution. And we set the failure maximum strain is 1%. We can observe that the 60 deg vessel 55

reached failure at 4 hours and 90 deg vessel reaches failure about 3 hours. The strain distributions at the dome area for both vessels are close to each other. The obvious discrepancy occurs in the straight part of the vessel and in this area the fiber orientation is important for the strain evolution and 90 deg has the weakest behavior. So it takes less time for 90 deg vessel to first get failure.

56

Path 2

Path 1

Fig 6.3 Pressure vessel, geometry and mesh and path 57

.0036

.0032 45 deg - 10 hours 45 deg - elastic step

.0035

60 deg - 10 hours 60 deg - elastic step

Maximum Principle Strain

Maximun Principle Strain

.0031

.0030

.0029

.0034

.0033

.0032

.0031

.0028 .0030

.0029

.0027 0

20

40

60

0

80

20

40

60

80

Path 1

Path 1

(a)

(b)

.0037

Maximum Principle Strain

.0036

.0035

.0034

.0033

90 deg - 10 hours 90 deg - elastic step

.0032

.0031 0

20

40

60

80

Path 1

(c)

Fig 6.4 Path plot of time evolution of Maximum Strain along path1 of the vessel under inside pressure 0.5Mpa with damage evolution (a) fiber orientation 45 deg. (b) fiber orientation 60 deg. (c) fiber orientation 90 deg.

58

.0038

Maximum Principle strain

.0036

.0034

.0032

.0030

.0028

45 deg 60 deg 90 deg

.0026 0

10

20

30

40

50

60

70

Path 1

Fig 6.5 Maximum Strain along path1 of the vessel with fiber orientation 0, 45, 90 deg under inside pressure 0.5Mpa with damage evolution after 10 hours

59

.0040 45 deg - 10 hours 45 deg - elastic step 60 deg - 10 hours 60 deg - elastic step 90 deg - 10 hours 90 deg - elastic step

Maximum Principle Strain

.0035

.0030

.0025

.0020

.0015 0

10

20

30

40

50

60

Path 2

Fig 6.6 Time evolution of Maximum Principle Strain along path2 of the pressure vessel with fiber orientation 45, 60, 90 deg under pressure 0.5Mpa with damage evolution

60

0.3%

0.2%

0.4%

0.4%

0.6%

0.6%

0.8%

0.8%

1%

1%

Failure area

Failure area

Failure area

Fig 6.7 Time evolution of Maximum Principle Strain of the pressure vessel with fiber orientation 60, 90 deg under pressure 1Mpa with damage evolution. (a) fiber orientation 60 deg. (b) fiber orientation 90 deg.

61

CHAPTER VII CONCLUSION

This paper demonstrates the computational utility of the anisotropic, creep damage model presented in a paper by Robinson, Binienda and Ruggles (2002). The references 1 and 2 describe supporting exploratory experiments are conducted on thin-walled tubular specimens fabricated from a model PMC. Thin-walled tubes are used not for their interest as structural components but because they are convenient specimens for generating multiaxial stress and deformation.

The computational routine (UMAT)

successfully describes the rather complex creep/damage phenomenon observed in PMCs. The routine is used with a commercially available code. The present work demonstrates the utility of the creep damage law in describing the essential physics behind creep damage using a single scalar parameter. The model is successfully applied to a benchmark problem (circular hole a square plate). A primary assumption in the damage model, cf., Robinson et al. (2002), is that the stress dependence of damage evolution is on the transverse tensile and longitudinal shear traction acting at the fiber/matrix interface. Accordingly, the isochronous damage function is taken to depend on the appropriate invariants N and S, i.e., ∆( N, S) . Exploratory data are generated to partially define the isochronous damage curve ∆( N , S) = 1

for the model PMC; evidently, a more extensive data base is required to fully 62

define ∆( N , S) = 1 and to verify that this stress dependence correlates directly with creep failure. This maybe of general interest and the present code needs to be extended to describe the general power law form of the damage curve. Also, the code may easily be extended to calculate creep rupture life based on the deformation rate and the damage variable calculations. These are left for future work.

63

REFERENCES

1. ROBINSON, D.N., BINIENDA, W.K. AND RUGGLES, M.B. (2002). “CREEP OF POLYMER MATRIX COMPOSITES: PART 1- A NORTON/BAILEY CREEP LAW FOR TRANSVERSE ISOTROPY.” JOURNAL OF ENGINEERING MECHANICS, Vol. 129, No. 3, March 2003, pp. 310-317 2. BINIENDA, W.K., ROBINSON, D.N. AND RUGGLES, M.B. (2002). “CREEP FAILURE OF POLYMER MATRIX COMPOSITES (PMC): A MONKMANGRANT RELATIONSHIP FOR TRANSVERSE ISOTROPY”. JOURNAL OF ENGINEERING MECHANICS, Vol. 129, No. 3, March 2003, pp. 318-323 3. LECKIE, F.A. (1986). “THE MICRO- AND MACRO-MECHANICS OF CREEP RUPTURE.” ENGRG. FRACTURE MECH., 25,5, 505-521. 4. LECKIE, F.A., AND HAYHURST, D.R. (1974)“CREEP RUPTURE IN STRUCTURES.” PROC. ROYAL SOCIETY OF LONDON, A340, 323-347. 5. MONKMAN, F.C. AND GRANT, N.J., (1956). “AN EMPIRICAL RELATIONSHIP BETWEEN RUPTURE LIFE AND MINIMUM CREEP RATE IN CREEP-RUPTURE TESTS.” ASTM, 56, 593 – 620. 6. LISSENDEN, C.J., LERCH, B.A., ELLIS, J.R. AND ROBINSON, D.N. (1997). “EXPERIMENTAL DETERMINATION OF YIELD AND FLOW SURFACES UNDER AXIAL-TORSIONAL LOADING.” STP 1280, A.S.T.M, 92-112. 7. ROBINSON, D.N. AND DUFFY, S.F. (1990). “CONTINUUM DEFORMATION THEORY FOR HIGH TEMPERATURE METALLIC COMPOSITES.” J. ENGRG. MECH., ASCE, 116(4), 832-844. 8. ROBINSON, D.N., BINIENDA, W.K. AND MITI-KAVUMA, M. (1992). “CREEP AND CREEP RUPTURE OF METALLIC COMPOSITES.” J. ENGRG. MECH.,. ASCE, 118(8), 1646-1660. 9. ROBINSON, D.N. AND PASTOR, M.S. (1993). “LIMIT PRESSURE OF A CIRCUMFERENTIALLY REINFORCED SIC/TI RING.” COMPOSITES ENGINEERING, 2, (4), 229-238. 64

10. ROBINSON, D.N., TAO, Q. AND VERRILLI, M.J. (1994). “A HYDROSTATIC STRESS-DEPENDENT ANISOTROPIC MODEL OF VISCOPLASTICITY.” NASA TM 106525. 11. ROBINSON, D.N. AND WEI, WEI (1996). “FIBER ORIENTATION IN COMPOSITE STRUCTURES FOR OPTIMAL RESISTANCE TO CREEP FAILURE.” J. ENGRG. MECH., ASCE, 122,(9), 855-860. 12. ROBINSON, D.N. AND BINIENDA, W.K., (2001)A. “OPTIMAL FIBER ORIENTATION IN CREEPING COMPOSITE STRUCTURES.” J. APPL. MECHANICS, 68,(2), 213-217. 13. ROBINSON, D.N. AND BINIENDA, W.K., (2001)B. “MODEL OF VISCOPLASTICITY FOR TRANSVERSELY ISOTROPIC INELASTICALLY COMPRESSIBLE SOLIDS.” J. ENGRG. MECH., ASCE, 127,(6). 14. ROBINSON, D.N., KIM, K.J AND WHITE, J.L. (2002). “CONSTITUTIVE MODEL OF A TRANSVERSELY ISOTROPIC BINGHAM FLUID.” J. APPL. MECHANICS, 69,(1), 1-8. 15. ABAQUS THEORY AND VERIFICATION MANUALS, VERSION 6.2, HKS INC.

65

APPENDICES

66

APPENDIX A UMAT FILE

SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD, * RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN, * TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MATERL,NDI,NSHR,NTENS, * NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT, * DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC) C INCLUDE 'ABA_PARAM.INC' C CHARACTER*80 MATERL DIMENSION STRESS(NTENS),STATEV(NSTATV), *

DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),

*

STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),

*

PROPS(NPROPS),COORDS(3),DROT(3,3),

*

DFGRD0(3,3),DFGRD1(3,3)

DOUBLE PRECISION J0,J,J2,KSI,NN,N,V,M C

ELASTIC PROPERTIES EMOD1=PROPS(1) 67

EMOD2=PROPS(2) ENU=PROPS(3) EG=PROPS(4) THETA=PROPS(5) KSI=PROPS(6) ESI=PROPS(7) ENTA=PROPS(8) N=PROPS(9) E0=PROPS(10) SIG0=PROPS(11) V=PROPS(12) M=PROPS(13) T0=PROPS(14) EBULK=1/(1.0-ENU**2*EMOD2/EMOD1) C C

ELASTIC STIFFNESS

C IF (KSTEP.EQ.1) THEN DO K1=1,NTENS DO K2=1,NTENS DDSDDE(K2,K1)=0 END DO ENDDO 68

C Q11=EMOD1*EBULK Q12=ENU*EMOD2*EBULK Q21=ENU*EMOD2*EBULK Q22=EMOD2*EBULK Q66=EG DDSDDE(1,1)=Q11*COS(THETA)**4.0 *

+2.0*(Q12+2*Q66)*SIN(THETA)**2.0*COS(THETA)**2.0

*

+Q22*SIN(THETA)**4.0

DDSDDE(1,2)=(Q11+Q22-4*Q66)*SIN(THETA)**2.0*COS(THETA)**2.0 *

+Q12*(SIN(THETA)**4.0+COS(THETA)**4.0) DDSDDE(2,1)=DDSDDE(1,2) DDSDDE(2,2)=Q11*SIN(THETA)**4.0

*

+2.0*(Q12+2*Q66)*SIN(THETA)**2.0*COS(THETA)**2.0

*

+Q22*COS(THETA)**4.0 DDSDDE(NTENS,1)=(Q11-Q12-2*Q66)*SIN(THETA)*COS(THETA)**3.0+

*

(Q12-Q22+2*Q66)*SIN(THETA)**3.0*COS(THETA)

DDSDDE(1,NTENS)=DDSDDE(NTENS,1) DDSDDE(NTENS,2)=(Q11-Q12-2*Q66)*SIN(THETA)**3.0*COS(THETA)+

*

(Q12-Q22+2*Q66)*SIN(THETA)*COS(THETA)**3.0

DDSDDE(2,NTENS)=DDSDDE(NTENS,2) 69

DDSDDE(NTENS,NTENS)=(Q11+Q22-2*Q12-2*Q66)*SIN(THETA)**2 * *

*COS(THETA)**2.0+ Q66*(SIN(THETA)**4+COS(THETA)**4.0)

C C1=DDSDDE(1,1) C2=DDSDDE(1,2) C3=DDSDDE(1,NTENS) C4=DDSDDE(2,1) C5=DDSDDE(2,2) C6=DDSDDE(2,NTENS) C7=DDSDDE(NTENS,1) C8=DDSDDE(NTENS,2) C9=DDSDDE(NTENS,NTENS) C C CALCULATE STRESS FROM ELASTIC STRAINS C DO K1=1,NTENS DO K2=1,NTENS STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1) END DO END DO C STATEV(1)=STATEV(1)+DSTRAN(1) 70

STATEV(2)=STATEV(2)+DSTRAN(2) STATEV(3)=STATEV(3)+DSTRAN(3) A1=STATEV(1) A2=STATEV(2) A3=STATEV(3) D1=C1 D2=C5 D3=C2 C END IF C C CALCULATE CREEP STRESS IF (KSTEP.EQ.2) THEN S11=STRESS(1)-(STRESS(1)+STRESS(2)+STRESS(3))/3.D0 S22=STRESS(2)-(STRESS(1)+STRESS(2)+STRESS(3))/3.D0 S33=STRESS(3)-(STRESS(1)+STRESS(2)+STRESS(3))/3.D0 S12=STRESS(NTENS) D11=COS(THETA)**2.D0 D22=SIN(THETA)**2.D0 D12=COS(THETA)*SIN(THETA) J0=D11*S11+D22*S22+2.D0*D12*S12 J=D11*(S11**2.D0+S12**2.D0)+D22*(S12**2.D0+S22**2.D0) * +D12*S12*(S11+S22) 71

J2=.5*(S11**2.D0+S22**2.D0+S33**2.D0)+S12**2 Q1=J-J0**2 Q2=J0**2 TEMP1=(ESI-4.D0*ENTA)*(STRESS(1)+STRESS(2))**2/9.D0 TEMP2=(KSI)*Q1 TEMP3=(ESI-ENTA)*Q2 TEMP4=3.0*(J2-TEMP2-TEMP3+TEMP1) PHI=SQRT(TEMP4)/(E0) TA11=2.D0*S11*D11+2.D0*D12*S12 TA22=2.D0*S22*D22+2.D0*D12*S12 TA12=S12*(D11+D22)+D12*(S11+S22) TB11=2.D0*J0*D11 TB22=2.D0*J0*D22 TB12=2.D0*J0*D12 TONE11=TA11-TB11 TONE22=TA22-TB22 TONE12=TA12-TB12 TTWO11=2.D0*J0*(D11-1.D0/3.D0) TTWO22=2.D0*J0*(D22-1.D0/3.D0) TTWO12=2.D0*J0*D12 TAO11=S11-KSI*TONE11-(ESI-ENTA)*TTWO11+ *

2.D0*(ESI-4.D0*ENTA)*(STRESS(1)+STRESS(2))/9.D0

TAO22=S22-KSI*TONE22-(ESI-ENTA)*TTWO22+ 72

*

2.D0*(ESI-4.D0*ENTA)*(STRESS(1)+STRESS(2))/9.D0 TAO12=S12-KSI*TONE12-(ESI-ENTA)*TTWO12 RSTRAN11=3.D0*PHI**(N-1.D0)*TAO11/(2.D0*E0*SIG0*100) RSTRAN22=3.D0*PHI**(N-1.D0)*TAO22/(2.D0*E0*SIG0*100) RSTRAN12=3.D0*PHI**(N-1.D0)*TAO12/(2.D0*E0*SIG0*100) DCRSTRAN11=RSTRAN11*DTIME DCRSTRAN22=RSTRAN22*DTIME DCRSTRAN12=RSTRAN12*DTIME

STATEV(4)=DCRSTRAN11+STATEV(4) STATEV(5)=DCRSTRAN22+STATEV(5) STATEV(6)=DCRSTRAN12+STATEV(6) c c DAMAGE calculation c I=STRESS(1)+STRESS(2)+STRESS(3) I0=D11*STRESS(1)+2*D12*STRESS(NTENS)+D22*STRESS(2) RE1=J2+.25*J0**2.0-J IF (RE1.LE.0) THEN RE1=0-RE1 END IF NN=.5*(I-I0)+SQRT(RE1) IF (NN.LE.0) THEN NN=0.0 73

END IF IF (Q1.LE.0) THEN Q1=J0**2-J END IF S=SQRT(Q1) DELTA=(NN+0.35*S)/(E0) STATEV(11)=(1-DELTA**V*TIME(2)/T0)**(1/(1+M)) RDASTRAN11=3.0*PHI**(N-1.0)*TAO11/(2.0*E0*SIG0*100) *

/(STATEV(11)**N) RDASTRAN22=3.0*PHI**(N-1.0)*TAO22/(2.0*E0*SIG0*100)

*

/(STATEV(11)**N) RDASTRAN12=3.0*PHI**(N-1.0)*TAO12/(2.0*E0*SIG0*100)

*

/(STATEV(11)**N) DDAMSTRAN11=RDASTRAN11*DTIME DDAMSTRAN22=RDASTRAN22*DTIME DDAMSTRAN12=RDASTRAN12*DTIME STATEV(12)=DDAMSTRAN11+STATEV(12)

STATEV(13)=DDAMSTRAN22+STATEV(13) STATEV(14)=DDAMSTRAN12+STATEV(14) c C CALCUALTE UPDATED STRESS C DS1=C1*(DSTRAN(1)-DCRSTRAN11) 74

DSTRESS11=C1*(DSTRAN(1)-DCRSTRAN11)+C2*(DSTRAN(2)DCRSTRAN22)+ *

C3*(DSTRAN(NTENS)-DCRSTRAN12) DSTRESS22=C4*(DSTRAN(1)-DCRSTRAN11)+C5*(DSTRAN(2)-

DCRSTRAN22)+ *

C6*(DSTRAN(NTENS)-DCRSTRAN12)

DSTRESS12=C7*(DSTRAN(1)-DCRSTRAN11)+C8*(DSTRAN(2)DCRSTRAN22)+ *

C9*(DSTRAN(NTENS)-DCRSTRAN12)

STRESS(1)=DSTRESS11+STRESS(1) STRESS(2)=DSTRESS22+STRESS(2) STRESS(NTENS)=DSTRESS12+STRESS(NTENS) C C CALCULATE UPDATED JACOBIAN DO K1=1,NTENS DO K2=1,NTENS DDSDDE(K2,K1)=0 END DO ENDDO Q11=EMOD1*EBULK Q12=ENU*EMOD2*EBULK Q21=ENU*EMOD2*EBULK Q22=EMOD2*EBULK 75

Q66=EG DDSDDE(1,1)=Q11*COS(THETA)**4.0 *

+2.0*(Q12+2*Q66)*SIN(THETA)**2.0*COS(THETA)**2.0

*

+Q22*SIN(THETA)**4.0

DDSDDE(1,2)=(Q11+Q22-4*Q66)*SIN(THETA)**2.0*COS(THETA)**2.0 *

+Q12*(SIN(THETA)**4.0+COS(THETA)**4.0) DDSDDE(2,1)=DDSDDE(1,2) DDSDDE(2,2)=Q11*SIN(THETA)**4.0

*

+2.0*(Q12+2*Q66)*SIN(THETA)**2.0*COS(THETA)**2.0

*

+Q22*COS(THETA)**4.0 DDSDDE(NTENS,1)=(Q11-Q12-2*Q66)*SIN(THETA)*COS(THETA)**3.0+

*

(Q12-Q22+2*Q66)*SIN(THETA)**3.0*COS(THETA)

DDSDDE(1,NTENS)=DDSDDE(NTENS,1) DDSDDE(NTENS,2)=(Q11-Q12-2*Q66)*SIN(THETA)**3.0*COS(THETA) *

(Q12-Q22+2*Q66)*SIN(THETA)*COS(THETA)**3.0

DDSDDE(2,NTENS)=DDSDDE(NTENS,2) DDSDDE(NTENS,NTENS)=(Q11+Q22-2*Q12-2*Q66)*SIN(THETA)**2 * *

*COS(THETA)**2.0+ Q66*(SIN(THETA)**4+COS(THETA)**4.0) RETURN END

76

APPENDIX B INPUT FILE FOR PRESSURE VESSEL

*Heading ** Job name: Job-1 Model name: Model-1 *Node 1,

0.,

20.,

0.

…… 12445, 19.9938431,

-20., 0.496222973

*Element, type=S8R,ELSET=CREEP 1, 1, 291, 399, 4, 4171, 4172, 4173, 4174 …… 4106, 4170, 164,

3, 163, 12425, 12445, 9401, 12444

*Nset, nset=B1 3, 164, ……., 12445 *Surface, type=ELEMENT, name=vessel creep, SNEG *shell SECTION,ELSET=CREEP,MATERIAL=com1 0.4 *TRANSVERSE SHEAR STIFFNESS 77

2230, 2230, 0 ** MATERIALS-1 *MATERIAL,NAME=COM1 *USER MATERIAL,CONSTANTS=16,TYPE=MECHANICAL 44e4, 7.3e3, .284, 2.7e3,1.047, 0.57, 0.64, 0.1, 6.5, 46., 50., 10.6, 7.125, 12,2.7e3,2.7e3 *DEPVAR 14 *INITIAL CONDITIONS,TYPE=SOLUTION CREEP,0.0,0.0,0.0,0.0,0.0,0.0,0.0 0.0,0.0,0.0,0.0,0.0,0.0,0.0 *Boundary B1, 2, 2 *STEP,INC=30 PRESCRIBED TENSILE STRESS *STATIC 1.E-7,1.E-6 *Dsload vessel, P, 0.5 *EL PRINT,FREQ=1 S, SDV1, SDV4 78

*END STEP *STEP,INC=100000 CREEP STEP *VISCO,CETOL=1E-4 1.E-6,10 ** OUTPUT REQUESTS *PRINT,FREQ=1,RESIDUAL=YES *EL PRINT,FREQ=1 SDV *OUTPUT,FIELD,FREQ =1 *ELEMENT OUTPUT S,E,CE, SDV *NODE OUTPUT U *END STEP

79