NASA Technical Reports Server (NTRS) 20150000901: Verification and Validation of a Three-Dimensional Generalized Composite Material Model

0.87 MB·English

by  NASA Technical Reports Server (NTRS)

NASA/TM—2015-218446
Veri(cid:191) cation and Validation of a Three-Dimensional
Generalized Composite Material Model
Canio Hoffarth, Joseph Harrington, and Subramaniam D. Rajan
Arizona State University, Tempe, Arizona
Robert K. Goldberg and Kelly S. Carney
Glenn Research Center, Cleveland, Ohio
Paul DuBois
George Mason University, Fairfax, Virginia
Gunther Blankenhorn
Livermore Software Technology Corporation, Livermore, California
January 2015
NASA STI Program . . . in Pro(cid:191) le
Since its founding, NASA has been dedicated to the • CONFERENCE PUBLICATION. Collected
advancement of aeronautics and space science. The papers from scienti(cid:191) c and technical
NASA Scienti(cid:191) c and Technical Information (STI) conferences, symposia, seminars, or other
program plays a key part in helping NASA maintain meetings sponsored or cosponsored by NASA.
this important role.
• SPECIAL PUBLICATION. Scienti(cid:191) c,
The NASA STI Program operates under the auspices technical, or historical information from
of the Agency Chief Information Of(cid:191) cer. It collects, NASA programs, projects, and missions, often
organizes, provides for archiving, and disseminates concerned with subjects having substantial
NASA’s STI. The NASA STI program provides access public interest.
to the NASA Aeronautics and Space Database and
its public interface, the NASA Technical Reports • TECHNICAL TRANSLATION. English-
Server, thus providing one of the largest collections language translations of foreign scienti(cid:191) c and
of aeronautical and space science STI in the world. technical material pertinent to NASA’s mission.
Results are published in both non-NASA channels
and by NASA in the NASA STI Report Series, which Specialized services also include organizing
includes the following report types: and publishing research results, distributing
specialized research announcements and feeds,
• TECHNICAL PUBLICATION. Reports of providing information desk and personal search
completed research or a major signi(cid:191) cant phase support, and enabling data exchange services.
of research that present the results of NASA
programs and include extensive data or theoretical For more information about the NASA STI
analysis. Includes compilations of signi(cid:191) cant program, see the following:
scienti(cid:191) c and technical data and information
deemed to be of continuing reference value. • Access the NASA STI program home page at
NASA counterpart of peer-reviewed formal http://www.sti.nasa.gov
professional papers but has less stringent
limitations on manuscript length and extent of • E-mail your question to [email protected]
graphic presentations.
• Phone the NASA STI Information Desk at
• TECHNICAL MEMORANDUM. Scienti(cid:191) c 757-864-9658
and technical (cid:191) ndings that are preliminary or
of specialized interest, e.g., quick release • Write to:
reports, working papers, and bibliographies that NASA STI Information Desk
contain minimal annotation. Does not contain Mail Stop 148
extensive analysis. NASA Langley Research Center
Hampton, VA 23681-2199
• CONTRACTOR REPORT. Scienti(cid:191) c and
technical (cid:191) ndings by NASA-sponsored
contractors and grantees.
NASA/TM—2015-218446
Veri(cid:191) cation and Validation of a Three-Dimensional
Generalized Composite Material Model
Canio Hoffarth, Joseph Harrington, and Subramaniam D. Rajan
Arizona State University, Tempe, Arizona
Robert K. Goldberg and Kelly S. Carney
Glenn Research Center, Cleveland, Ohio
Paul DuBois
George Mason University, Fairfax, Virginia
Gunther Blankenhorn
Livermore Software Technology Corporation, Livermore, California
Prepared for the
13th US and International LS-DYNA Users Conference
sponsored by the Livermore Software Technology Corporation
Dearborn, Michigan, June 8–10, 2014
National Aeronautics and
Space Administration
Glenn Research Center
Cleveland, Ohio 44135
January 2015
Acknowledgments
Authors Hoffarth, Harrington and Rajan gratefully acknowledge the support of the Federal Aviation Administration through
Grant #12-G-001 entitled “Composite Material Model for Impact Analysis”, William Emmerling, Technical Monitor.
Trade names and trademarks are used in this report for identi(cid:191) cation
only. Their usage does not constitute an of(cid:191) cial endorsement,
either expressed or implied, by the National Aeronautics and
Space Administration.
Level of Review: This material has been technically reviewed by technical management.
Available from
NASA STI Information Desk National Technical Information Service
Mail Stop 148 5301 Shawnee Road
NASA Langley Research Center Alexandria, VA 22312
Hampton, VA 23681-2199
Available electronically at http://www.sti.nasa.gov
Verification and Validation of a Three-Dimensional Generalized
Composite Material Model

Canio Hoffarth, Joseph Harrington, and Subramaniam D. Rajan
Arizona State University
Tempe, Arizona 85004

Robert K. Goldberg and Kelly S. Carney
National Aeronautics and Space Administration
Glenn Research Center
Cleveland, Ohio 44135

Paul DuBois
George Mason University
Fairfax, Virginia 22030

Gunther Blankenhorn
Livermore Software Technology Corporation
Livermore, California 94551
Abstract
A general purpose orthotropic elasto-plastic computational constitutive material model has been
developed to improve predictions of the response of composites subjected to high velocity impact. The
three-dimensional orthotropic elasto-plastic composite material model is being implemented initially for
solid elements in LS-DYNA as MAT213. In order to accurately represent the response of a composite,
experimental stress-strain curves are utilized as input, allowing for a more general material model that can
be used on a variety of composite applications. The theoretical details are discussed in a companion
paper. This paper documents the implementation, verification and qualitative validation of the material
model using the T800-F3900 fiber/resin composite material.
Introduction
Composite materials are now beginning to provide uses hitherto reserved for metals in structural
systems such as airframes and engine containment systems, wraps for repair and rehabilitation and
ballistic/blast mitigation systems. While material models exist that can be used to simulate the response of
materials in these demanding structural systems, they are often designed for use with specific
homogeneous materials such as metals (Refs. 1 and 2), polymers (Ref. 3) and wood (Ref. 4). While
material models exist to simulate the impact response of composites, most are tailored specifically for a
particular application and have limitations - purely elastic, no rate sensitivity, implementation for solid
elements only, and limited damage and failure characterization (Refs. 5 to 10).
A few examples of different approaches to composite damage modeling include Littell et al.
(Ref. 11), Matzenmiller et al. (Ref. 12), and others (Refs. 13 to 18). Yen’s comprehensive model (Ref. 19)
is implemented in LS-DYNA as *MAT_162. Availability of constitutive models for composite materials
in LS-DYNA is discussed in some length in our companion paper (Ref. 20) and that explains the
motivation for this study.

NASA/TM—2015-218446 1
In the next section we briefly discuss the theoretical development of the constitutive model. This is
followed by details of the numerical algorithm, verification tests for a composite used in the aerospace
industry, the use of the composite in a high speed contact example and finally, concluding remarks
including details of our ongoing and future work.
Orthotropic 3D Elasto-Plastic Composite Material Model
The material law that the model is built upon describes the elastic and permanent deformation of the
composite with full three-dimensional implementation suitable for solid elements. Only the relevant
equations used in the numerical algorithm are presented here. The theoretical basis and discussions are
available in a companion paper (Ref. 20). Current development of the model includes the complete rate
independent elasto-plastic deformation model, with damage and failure to be added later.
The commonly used Tsai-Wu composite failure model is generalized as a yield surface for the
plasticity model and is defined as

T
(cid:170)(cid:86) (cid:186) (cid:170)(cid:86) (cid:186) (cid:170)F F F 0 0 0 (cid:186)(cid:170)(cid:86) (cid:186)
11 11 11 12 13 11
(cid:171) (cid:187) (cid:171) (cid:187) (cid:171) (cid:187)(cid:171) (cid:187)
(cid:86) (cid:86) F F F 0 0 0 (cid:86)
(cid:171) 22(cid:187) (cid:171) 22(cid:187) (cid:171) 12 22 23 (cid:187)(cid:171) 22(cid:187)
(cid:171)(cid:86) (cid:187) (cid:171)(cid:86) (cid:187) (cid:171)F F F 0 0 0 (cid:187)(cid:171)(cid:86) (cid:187) (1)
f((cid:86))(cid:32)a(cid:14)(cid:11)F1 F2 F3 0 0 0(cid:12)(cid:171)(cid:86)33(cid:187)(cid:14)(cid:171)(cid:86)33(cid:187) (cid:171) 013 023 033 F 0 0(cid:187)(cid:171)(cid:86)33(cid:187)
(cid:171) 12(cid:187) (cid:171) 12(cid:187) (cid:171) 44 (cid:187)(cid:171) 12(cid:187)
(cid:171)(cid:86) (cid:187) (cid:171)(cid:86) (cid:187) (cid:171) 0 0 0 0 F 0(cid:187)(cid:171)(cid:86) (cid:187)
23 23 55 23
(cid:171) (cid:187) (cid:171) (cid:187) (cid:171) (cid:187)(cid:171) (cid:187)
(cid:172)(cid:86)31(cid:188) (cid:172)(cid:86)31(cid:188) (cid:172) 0 0 0 0 0 F66(cid:188)(cid:172)(cid:86)31(cid:188)

where a = –1 and

1 1 1 1
F1(cid:32)(cid:86)1T1(cid:16)(cid:86)1C1 F11(cid:32) (cid:86)1T1(cid:86)1C1 F44 (cid:32) (cid:86)122 (1a)
1 1 1 1
F2 (cid:32)(cid:86)T22 (cid:16)(cid:86)C22 F22 (cid:32) (cid:86)T22(cid:86)C22 F55 (cid:32) (cid:86)223
1 1 1 1
F (cid:32) (cid:16) F (cid:32) F (cid:32)
3 (cid:86)T (cid:86)C 33 (cid:86)T (cid:86)C 66 (cid:86)2
33 33 33 33 31
2 Fi(cid:14)Fj 1
F (cid:32) (cid:16) (cid:16) (F (cid:14)F (cid:14)F ) i, j(cid:32)1,2,3,k (cid:32)i(cid:14)3
ij (cid:11)(cid:86)i(cid:16)j(cid:12)2 (cid:86)i4(cid:16)5j 2 ii jj kk
45

The stress components of the yield function coefficients correspond to the current yield stresses
associated with a set of experimental stress-strain curves, which vary with the effective plastic strain, thus
allowing the model to describe different hardening properties in each direction. The non-associative flow
surface is defined as

h(cid:32) H (cid:86)2 (cid:14)H (cid:86)2 (cid:14)H (cid:86)2 (cid:14)2H (cid:86) (cid:86) (cid:14)2H (cid:86) (cid:86) (cid:14)2H (cid:86) (cid:86) (cid:14)H (cid:86)2 (cid:14)H (cid:86)2 (cid:14)H (cid:86)2 (2)
11 11 22 22 33 33 12 11 22 23 22 33 31 33 11 44 12 55 23 66 31

where the H terms are the constant flow rule coefficients. The rest of the details can be found in the
ij
companion paper (Ref. 20).
The first six flow rule coefficients are computed directly from the assumed flow rule coefficient
value, H and the plastic Poisson’s ratios.
11

NASA/TM—2015-218446 2
p p
(cid:16)(cid:81) (cid:81)
H (cid:32)(cid:16)(cid:81)p H H (cid:32) 23 12 H
12 12 11 23 p 11
(cid:81)
H (cid:32)(cid:16)(cid:81)pH 21 (3a)
13 13 11 p
(cid:81)
(cid:81)p H (cid:32) 13 H
H (cid:32) 12 H 33 (cid:81)p 11
22 p 11 31
(cid:81)
21

The last three flow rule coefficients (H , H , H ) are calculated iteratively with the user-defined
44 55 66
range to best fit the input experimental shear stress versus shear strain curves (Ref. 20). For example, to
determine H , the following optimization problem is posed using the input curve for the “1-2” shear case.
44
Find H to minimize
44

2
n
(cid:170)(cid:11) (cid:12) (cid:11) (cid:12) (cid:186)
f(H )(cid:32)(cid:166) (cid:86)ˆ (cid:16) (cid:86)ˆ (3b)
44 (cid:171) 22 12 (cid:187)
(cid:172) i i (cid:188)
i(cid:32)1

such that

Hmin (cid:100)H (cid:100) Hmax (3c)
44 44 44

(cid:11) (cid:12)
where n is the number of data points on the master curve, (cid:86)ˆ is the ith effective stress value from the
22
i
(cid:11) (cid:12)
master curve, and (cid:86)ˆ is the effective stress value for the shear curve using an assumed value of H .
12 44
i
Numerical Algorithm
In this initial implementation, we do not account for rate and temperature dependence, damage
accumulation nor failure. The following experimental data are needed as LS-DYNA input.

(1) A total of 12 true stress versus true strain curves under quasi-static test conditions from (a) uniaxial
tension in 1, 2 and 3-directions, (b) uniaxial compression in 1, 2 and 3-directions, (c) shear in 1-2,
2-3 and 3-1 planes and (d) off-axis (e.g., 45(cid:113)) uniaxial tension or compression in 1-2, 2-3 and 3-1
planes. The tabulated x-y data for each curve is read in and used appropriately.
(2) The modulus of elasticity, Poisson’s ratio and average plastic Poisson’s ratio obtained from the
tension and compression tests.

The effective plastic strain can be written in terms of the experimental stress versus total strain data.
In the 1-direction for example, this is written as

(cid:167) (cid:86)c (cid:183)(cid:189)
(cid:86)c (cid:168)(cid:72)p (cid:32)(cid:72) (cid:16) 11(cid:184)(cid:176)
11(cid:168)(cid:169) 11 11 E11(cid:184)(cid:185)(cid:176)(cid:190)(cid:159)(cid:86)1c1(cid:11)(cid:72)ep(cid:12) (4)
(cid:72)p (cid:32)(cid:179)(cid:11)(cid:86) d(cid:72)p /h(cid:12) (cid:176)
e 11 11 (cid:176)(cid:191)

where (cid:86)c is the experimental compressive true stress in the 1-direction, (cid:72) is the total true strain in the
11 11
1-direction, E is the elastic modulus in the 1-direction, (cid:72)p is the true plastic strain in the 1-direction,
11 11
NASA/TM—2015-218446 3
(cid:72)p is the effective plastic strain and h is the current value of the plastic potential function, which is
e
equivalent to the effective stress, shown in Equation (2). The material model uses a typical elastic stress
update as

(cid:305)n(cid:14)1 (cid:32)(cid:305)n (cid:14)C(cid:39)t:(cid:11)(cid:304)(cid:5)(cid:16)(cid:304)(cid:5)p(cid:12) (5)

where C is the orthotropic stiffness matrix, (cid:39)t is the time step, (cid:304)(cid:5) is the total strain rate and (cid:304)(cid:5)p is the
plastic strain rate. The stiffness matrix is written in terms of the compliance matrix as

(cid:16)1
(cid:170) 1 v v (cid:186)
(cid:16) 21 (cid:16) 31 0 0 0
(cid:171) (cid:187)
E E E
(cid:171) 11 22 33 (cid:187)
(cid:171) 1 v (cid:187)
(cid:16) 32 0 0 0
(cid:171) (cid:187)
E E
(cid:171) 22 33 (cid:187)
(cid:171) 1 (cid:187)
0 0 0
(cid:171) (cid:187)
E
C(cid:32)S(cid:16)1(cid:32)(cid:171) 33 (cid:187) (6)
(cid:171) 1 (cid:187)
0 0
(cid:171) (cid:187)
G
(cid:171) 12 (cid:187)
(cid:171) 1 (cid:187)
Sym 0
(cid:171) (cid:187)
G
(cid:171) 23 (cid:187)
(cid:171) 1 (cid:187)
(cid:171) (cid:187)
G
(cid:172) 31(cid:188)

where E are the elastic moduli in the principal material directions, G are the elastic shear moduli and
ij ij
v are the elastic Poisson’s ratios. The yield stresses used to determine the flow rule coefficients are
ij
summarized into a single vector, corresponding to each of the experimental test curves as
qT (cid:32)(cid:170)(cid:86)T (cid:86)T (cid:86)T (cid:86)C (cid:86)C (cid:86)C (cid:86) (cid:86) (cid:86) (cid:86)C (cid:86)C (cid:86)C (cid:186) . Lastly, the consistency
(cid:172) 11 22 33 11 22 33 12 23 31 45(cid:16)12 45(cid:16)23 45(cid:16)31(cid:188)
condition is written in terms of the rate of yield function as follows

(cid:119)f (cid:119)f
f(cid:5) (cid:32) (cid:305)(cid:5) (cid:14) q(cid:5) (cid:32)0 (7)
(cid:119)(cid:305) (cid:119)q

where the terms are as defined earlier. Equation (7) can be expanded as follows by introducing the
elasto-plastic constitutive law

f(cid:5) (cid:32) (cid:119)f (cid:167)(cid:168)C:(cid:304)(cid:5)(cid:16)C:(cid:79)(cid:5) (cid:119)h(cid:183)(cid:184)(cid:14) (cid:119)f (cid:79)(cid:5) dq (cid:32)0 (8)
(cid:119)(cid:305)(cid:169) (cid:119)(cid:305)(cid:185) (cid:119)q d(cid:79)

where (cid:305)(cid:5) is written in terms of the constitutive matrix and strain rates and (cid:79) is the effective plastic strain.
Solving for the effective plastic strain rate produces the following equation.

NASA/TM—2015-218446 4
(cid:119)f
C:(cid:304)(cid:5)
(cid:5) (cid:119)(cid:305)
(cid:79)(cid:32) (9)
(cid:119)f (cid:119)h (cid:119)f dq
C: (cid:14)
(cid:119)(cid:305) (cid:119)(cid:305) (cid:119)q d(cid:79)

To numerically solve for the increment in effective plastic strain, a radial return algorithm is used.
Initially, a perfectly elastic response is assumed. Thus, an elastic predictor is used to compute the stresses
as

(cid:11)(cid:86) (cid:12)e (cid:32)(cid:11)(cid:86) (cid:12)n (cid:14)C (cid:39)t(cid:11)(cid:72)(cid:5) (cid:12)(cid:14)C (cid:39)t(cid:11)(cid:72)(cid:5) (cid:12)(cid:14)C (cid:39)t(cid:11)(cid:72)(cid:5) (cid:12)
11 11 11 11 12 22 13 33
(cid:11)(cid:86) (cid:12)e (cid:32)(cid:11)(cid:86) (cid:12)n (cid:14)C (cid:39)t(cid:11)(cid:72)(cid:5) (cid:12)(cid:14)C (cid:39)t(cid:11)(cid:72)(cid:5) (cid:12)(cid:14)C (cid:39)t(cid:11)(cid:72)(cid:5) (cid:12)
22 22 12 11 22 22 23 33
(cid:11)(cid:86) (cid:12)e (cid:32)(cid:11)(cid:86) (cid:12)n (cid:14)C (cid:39)t(cid:11)(cid:72)(cid:5) (cid:12)(cid:14)C (cid:39)t(cid:11)(cid:72)(cid:5) (cid:12)(cid:14)C (cid:39)t(cid:11)(cid:72)(cid:5) (cid:12)
33 33 13 11 23 22 33 33
(10)
(cid:11)(cid:86) (cid:12)e (cid:32)(cid:11)(cid:86) (cid:12)n (cid:14)C (cid:39)t(cid:11)(cid:72)(cid:5) (cid:12)
12 12 44 12
(cid:11)(cid:86) (cid:12)e (cid:32)(cid:11)(cid:86) (cid:12)n (cid:14)C (cid:39)t(cid:11)(cid:72)(cid:5) (cid:12)
23 23 55 23
(cid:11)(cid:86) (cid:12)e (cid:32)(cid:11)(cid:86) (cid:12)n (cid:14)C (cid:39)t(cid:11)(cid:72)(cid:5) (cid:12)
31 31 66 31

With the elastic trial stresses computed, a trial yield function value can be calculated and used to
determine if the load step is elastic or plastic. It should be noted that the elastic trial stress is used as the
stress in the next step if the trial yield function value is less than or equal to zero. Otherwise, the stress
state is beyond the yield surface and must be brought back using a radial return and plastic corrector, (cid:39)(cid:79).
If the trial yield function happens to be greater than zero, then (cid:39)(cid:79) must be greater than zero. The value of
(cid:39)(cid:79) is determined using a secant iteration, with (cid:39)(cid:79)1(cid:32)0 for the first iteration and the second value for the
increment of the plastic multiplier determined from consistency,

(cid:119)f :(cid:11)(cid:305)e (cid:16)(cid:305)n(cid:12) (cid:119)f :(cid:11)(cid:305)e (cid:16)(cid:305)n(cid:12)
(cid:119)(cid:305) (cid:119)(cid:305)
(cid:39)(cid:79)2 (cid:32) e (cid:124) e (11)
(cid:119)f (cid:119)h (cid:119)f dq (cid:119)f (cid:119)h
C: (cid:14) C:
(cid:119)(cid:305) (cid:119)(cid:305) (cid:119)q d(cid:79) (cid:119)(cid:305) (cid:119)(cid:305)
e e e e

where the derivatives of q are taken as zero, meaning that the response is perfectly plastic. The initial
calculation of the second increment corresponds to perfect plasticity in order to ensure that the stress state
returns to the interior of the yield surface, resulting in a negative value of the yield function and bounding
the solution. From each estimate of (cid:39)(cid:79), the corresponding corrected stresses can be computed as well as
the yield stresses obtained from the input stress versus effective plastic strain ((cid:79)) data. With these terms
calculated for both estimates of(cid:39)(cid:79), the yield function value for each can be computed with the third
estimate evaluated as

(cid:39)(cid:79)2 (cid:16)(cid:39)(cid:79)1
(cid:39)(cid:79)3 (cid:32)(cid:39)(cid:79)1(cid:16) f1 (12)
f2 (cid:16) f1

The corrected stresses and yield stresses associated with (cid:507)(cid:540)3 are then calculated and used to obtain a
new estimate of the yield function value. Convergence of the secant iteration is determined by the
following conditions, where f is the value of the yield function given the effective plastic strain (cid:507)(cid:540)3,
3
NASA/TM—2015-218446 5
(cid:173)(cid:176)(cid:39)(cid:79)1(cid:32)(cid:39)(cid:79)3 (cid:173)(cid:176)(cid:39)(cid:79)1(cid:32)(cid:39)(cid:79)1
f3 (cid:33)0(cid:159)(cid:174) ; f3(cid:31)0(cid:159)(cid:174) ; f3 (cid:124)0(cid:159)(cid:39)(cid:79)(cid:32)(cid:39)(cid:79)3 (13)
(cid:176)(cid:175)(cid:39)(cid:79)2 (cid:32)(cid:39)(cid:79)2 (cid:176)(cid:175)(cid:39)(cid:79)2 (cid:32)(cid:39)(cid:79)3

with the secant iteration continuing if the new yield function value is not less than some defined tolerance.
Once convergence is met, the plastic multiplier increment to return the stress state to the yield surface is
known and the stresses can be updated as

(cid:119)h
(cid:305)n(cid:14)1(cid:32)(cid:305)e(cid:16)C:(cid:39)(cid:79) (14)
(cid:119)(cid:305)
e

Finally, the yield stresses are updated as well, using the new value of the overall effective plastic
strain, (cid:79), in each input curve to determine the corresponding yield stress level as

qn(cid:14)1(cid:32)q(cid:11)(cid:79)n (cid:14)(cid:39)(cid:79)3(cid:12) (15)

This results in anisotropic strain hardening, as a yield stress increase in each direction is initiated with
plasticity in any direction, but at different levels. A detailed algorithm is presented below.

Step 1: Preprocessing
(a) Convert stress/strain input curves to stress versus effective plastic strain using Equation (4).
(b) Store initial yield stresses in q.
(c) Calculate and store the elastic moduli from the initial yield stress and strain values.
(d) Store the three elastic Poisson’s ratio values and flow rule coefficients from input.
(e) Compute optimal values of the flow rule coefficients so as to match the input curves as closely as
possible.

Step 2: Initialization
Parameters available at start of step: q , (cid:305) , ((cid:304)(cid:5) ,(cid:39)t )(cid:111)(cid:39)(cid:304) , (cid:79) .
n n n n n n

Step 3: Elastic predictor
Compute and set the yield function coefficients using Equation (1a).
Construct the constitutive matrix using Equation (6).
Compute elastic trial stresses, (cid:305)e , using Equation (10).
n(cid:14)1
Compute the trial yield function, ftrial, using the elastic trial stresses in Equation (1).
n(cid:14)1
If ftrial (cid:100)0, the elastic solution is correct, set (cid:39)(cid:79) (cid:32)0 and go to stress update. Else go to plastic
n(cid:14)1 n
corrector.

Step 4: Plastic corrector
Set (cid:39)(cid:79)1 (cid:32)0.
Calculate (cid:39)(cid:79)2 from Equation (11).
Loop through secant iteration for specified max number of iterations:
(cid:11) 1 2(cid:12)
Compute the new estimate of the stress for each plastic multiplier increment (cid:39)(cid:79) ,(cid:39)(cid:79) using
Equation (14).
NASA/TM—2015-218446 6