Cyclic plasticity modeling and finite element analyzes of a circumferentially notched round bar under combined axial and torsion loadings

Materials and Design 34 (2012) 842–852

Contents lists available at ScienceDirect

Materials and Design journal homepage: www.elsevier.com/locate/matdes

Cyclic plasticity modeling and finite element analyzes of a circumferentially notched round bar under combined axial and torsion loadings Mehmet Firat ⇑ The University of Sakarya, Dept. of Mech. Engineering, 54187 Adapazari, Turkey

a r t i c l e

i n f o

Article history: Received 1 June 2011 Accepted 7 July 2011 Available online 2 August 2011 Keywords: A. Ferrous metals and alloys F. Plastic behavior H. Failure analysis

a b s t r a c t A cyclic plasticity model suitable for fatigue damage modeling of metallic structures is presented and its finite element (FE) implementation is described within the small strain plasticity framework. The model uses the von Mises yield surface in stress space whose evaluation follows an Armstrong–Frederick type of nonlinear kinematic hardening rule, and the normality hypothesis in conjunction with the associative flow rule is assumed. An incremental implicit-iterative algorithm was employed for the numerical solution of resulting stress–strain equations, and the continuum tangent obtained from plasticity model was used in FE implementation. The developed FE computational model is applied in the cyclic deformation analysis of a circumferentially notched specimen in combined axial force-torsion loading tests. The computed notch root deformations were compared with measured notch root strain histories. An assessment of model predictions showed that non-proportional loading tests have been simulated with a good accuracy. The computed strain loops were in accord with experimental data and matched qualitatively with measured shear – axial strain histories irrespective of loading path of the test. In proportional balanced torsion-axial loading, the nonlinear shear strain – axial strain loops were also simulated properly. The errors in notch root strains were dependent on the loading path shape, and compared to axial strains, the shear strain errors were relatively greater. The computer solution times were also acceptable. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction The fatigue-safe design of load-bearing structural parts is a critical ingredient of manufacture of durable products in automotive and related industries. While component testing under accelerated service programs is still the principal means of product sign-off for ground vehicles, the need for efficient computational methodologies has been growing to ensure durability early in the computer engineering phase [1]. In the last decades, local stress–strain approaches have been established as practical engineering tool for fatigue life predictions in several industrial design applications in this context [2]. In local approach, it is assumed that the smooth specimens tested under strain control can simulate fatigue damage in the fatigue-critical location of a component if both the highlystressed region and the specimen experience identical stress and strain histories. Furthermore, a multiaxial damage criterion is postulated on the basis of stress and/or strain components in order to describe a functional relation for the equivalent fatigue damage that material is subjected in both cases [3]. The local stress–strain evaluations on the basis of fatigue loading cycles constitutes basic input of the fatigue damage criterion, and therefore, an accurate calculation of material deformation response is essential in ⇑ Tel.: +90 264 2955449; fax: +90 264 2955450. E-mail address: fi[email protected] 0261-3069/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.matdes.2011.07.022

assessing cyclic damage at the critical design features such as notches, fillets, welds, and shoulders [4]. The calculation of local stress–strain response at notches is a standard boundary value problem that can be solved by satisfying the equilibrium equations, the material constitutive equation, and compatibility equations together with the boundary conditions. Early research on the notch stress–strain analysis was focused on the determination of theoretical stress concentration factors that were computed with elasticity theory or determined from photoelasticity analysis [5]. The use of elastic stress concentration factors in connection with superposition technique may yield accurate stresses and strains under multiaxial loading conditions as long as the maximum stress on the net section remains less than material yield stress [6]. Beyond elasticity, however, analytical solutions were limited to simple monotonic loadings and notch geometries leading uniaxial stress or strain state under plane stress and plane strain deformation conditions [7,8], and consequently, is not suitable for multiaxial elasto-plastic stress–strain analysis [9]. On the other hand, numerical techniques, in particular finite element (FE) methods, provide reliable solutions that can be applied in the variety of component geometry under general fatigue loadings [10]. With the remarkable increase in the computer power and computational capabilities of FE software, incremental FE solutions for fairly complex loading histories are now feasible in the stress and strain analysis for the purpose of cyclic damage

843

M. Firat / Materials and Design 34 (2012) 842–852

assessment in an industrial environment. An accurate fatigue life analysis requires, however, the use of plasticity constitutive models suitable for simulating complex multiaxial cyclic deformations [11]. Experimental investigations with smooth bar and thinwalled tube metallic specimens have demonstrated that the actual material stress–strain response exhibits fairly complex features such as Baushinger effect, cyclic hardening and cyclic softening, ratcheting or stress shakedown and nonproportional hardening [11–15]. Consequently, it is the plasticity model that determines the efficiency of FE simulation, and in order to achieve a high accuracy for a particular application, a proper plasticity modeling as well as its FE implementation is a compulsory task in this computational fatigue approach. The cyclic plasticity modeling is not new in the literature particularly for the low-cycle fatigue research [16,17]. Restricting to isothermal deformation at room temperature, the cyclic stress–strain response of metallic materials depends primarily on the history of deformation and strain-rate effects are usually neglected for many engineering metals and alloys [18–20]. At the macroscopic scale, the material is uniform and homogenous and can be assumed initial isotropic. Furthermore, there exists an initial yield stress that distinguishes the limit of elastic and elasto-plastic deformation in simple tension test. In literature, constitutive equations describing material stress–strain response under aforementioned conditions are known as rate-independent plasticity models [18–20]. In most of these models, the notion of the kinematic hardening in conjunction with a backstress variable is well-established in the description of the experimentally observed anisotropy of deformation [18–20]. Prager–Ziegler linear kinematic hardening [18], Besseling sublayer [17] and Mroz’s multi-surface theories [12], the two-surface models of Dafalias and Popov [21] and Krieg [22] and Armstrong–Frederick nonlinear kinematic hardening model [23] can be considered as the fundamental theories in this context. There have been many enhancements in recent decades including kinematic models developed by Chaboche [15], Onno and Wang [24] and Jiang and Sehitoglu [25], and an evaluation of metal plasticity models in literature indicates that various possibilities exist for both qualitatively and quantitatively improved description of the cyclic deformations material experience under fatigue loadings. Furthermore, some of the models have been also devised for fatigue life analyses and good results were also reported in several engineering applications [26–30]. In this study, a rate-independent plasticity model suitable for cyclic deformation modeling of metallic structures is presented and its FE-implementation is described within the small strain plasticity framework. The model is based on the multi-component additive decomposition of the total backstress and uses the von Mises yield surface in stress space whose evaluation follows a Armstrong–Frederick type of nonlinear kinematic hardening rule. An incremental implicit-iterative algorithm is employed for the numerical solution of resulting stress–strain plasticity equations and implemented into Ansys FE program [31]. The presented computational model is applied in the FE simulation of a circumferentially notched specimen loaded under proportional and nonproportional axial force-torsion cycles. The computed notch root deformations are compared with experimental data in terms of measured notch root strain components. An assessment of FE predicted strain histories indicates a good agreement with measurements and proved the performance of the developed computational model both qualitatively and quantitatively.

2. Cyclic deformation modeling Among many constitutive models, the nonlinear kinematic hardening model of Armstrong and Frederick [23] became one

of mostly respected plasticity framework in the cyclic deformation modeling of engineering metals and alloys. With the primary aim of describing the nonlinear stress–strain hysteresis loop under cyclic straining, Armstrong and Frederick modified the linear kinematic hardening model [23] by introducing a recovery term related to the strain memory effect into the backstress evaluation.

da ¼ adep  badp

ð1Þ p

In this equation, de and dp denote the increment of plastic strain tensor and equivalent plastic strain increment respectively, and ais a material constant and b is a scalar function of the plastic strain path. Due to the explicit dependence of yield surface translation, da, to the current yield surface center, a, the hardening rule is called as the nonlinear kinematic hardening rule [18–20]. The Armstrong and Frederick (A–F) model can simulate both the Bauschinger effect and the nonlinear stress–strain curve successfully under cyclic deformations. However, stress–strain hysteresis loops under tension– compression cycles were not accurate [19,20]. Chaboche [15] expressed the backstress tensor as a series expansion of several back stress components each having the same structure of the A–F model.

a¼

m X

aðiÞ

ð2Þ

i¼1

daðiÞ ¼

2 ðiÞ p C de  cðiÞ adp 3

ð3Þ

According to this concept, a represents the total backstress composed of m individual parts and a(i) is a single backstress component having its own set of model parameters, C(i) and c(i). The model parameters are calculated using the material stabilized cyclic stress–strain curve determined by means of cyclic strain tests at different strain amplitudes. In this multi-component form of A–F model, the mean stress relaxation and cyclic strain accumulation (ratcheting) phenomena observed under asymmetric stress- and strain-controlled loading conditions were described properly. The decision on the number of backstress components is related to the shape of the cyclic curve of a given material, and the use of three to five back stress components are considered to be accurate for most of the metals [15]. This model has been further integrated with isotropic hardening rule in order to simulate the cyclic hardening and other modeling enhancements have been introduced by Chaboche and his coworkers [18]. Onno and Wang [24] introduced a critical state concept in the modeling of the high-cycle ratcheting phenomenon under load-controlled, asymmetric cyclic loadings. Based on the Chaboche’s additive decomposition of total backstress, Onno and Wang modified the dynamic recovery term in backstress evaluation equation with a critical state parameter, dp(i)

daðiÞ ¼

2 ðiÞ p ðiÞ h de  fðiÞ adp 3

ð4Þ

As in the Chaboche model, h(i) and f(i) are material parameters for each backstress part. dp(i) stands for the part of equivalent plastic strain increment causing dynamic recovery during cyclic loading and defined with the following inequality [24]:

0 6 dp

ðiÞ

6

 12 2 p de : dep 3

ð5Þ

Jiang-Sehitoglu [25] introduced the concept of limiting surfaces for each backstress components, and unified the A–F type models previously developed by Chaboche and Onno–Wang on this basis. According to the Jiang–Sehitoglu kinematic model [25], the total backstress is composed of m number of backstress components whose evaluation is given with the following expression

844

M. Firat / Materials and Design 34 (2012) 842–852

0 ðiÞ

da

1

!vðiÞ þ1 kaðiÞ k n LðiÞ Adp ði ¼ 1; 2; . . . ; mÞ r ðiÞ

ðiÞ ðiÞ @

¼c r

ð6Þ

where n is normal tensor of the yield surface at the current stress, and the quantity L(i) is the unit tensor of the ith backstress component defined by

LðiÞ ¼

aðiÞ ði ¼ 1; 2; . . . ; mÞ kaðiÞ k

ð7Þ

c(i) and r(i) are the model parameters for each backstress part and the exponent v(i) is introduced in order to control the ratcheting rate. Using the limiting surface concept, the determination of material parameters in the backstress evaluation equation becomes relatively simple compared to the graphical or numerical approaches previously required for Chaboche and Onno–Wang models. An investigation of recent fatigue studies in literature shows the relative success of A–F type nonlinear kinematic hardening models in cyclic deformation analysis of several engineering metals and alloys at room temperature applications [11,15,17, 24,27,29]. Recognizing the efficiency in additive backstress evaluation and simplicity in calculation of kinematic model parameters, Jiang–Sehitoglu model is chosen for the multiaxial stress–strain analysis based on the FE method and implemented into Ansys FE program in this study. 2.1. A rate-independent cyclic plasticity model In small-strain plasticity, the increment of linearized Green strain tensor de is decomposed additively into elastic and plastic strain increments

de ¼ dee þ dep

ð8Þ

The same strain increment may also be divided into deviatoric and hydrostatic parts as follows:

de ¼ de þ dem I

ð9Þ

where I is the second-order identity tensor. de denotes the incremental deviatoric strain and dem is mean strain increment calculated with

dem ¼

1 1 ðde : IÞ ¼ trðdeÞ 3 3

ð10Þ

The plastic deformation of metals is incompressible, depm ¼ 0. As a result, the plastic strains are deviatoric (distortional) and the mean strain is purely elastic and recoverable, dem ¼ deem . Then, Eq. (9) may also be written as:

de ¼ dee þ dep þ dem I

ð11Þ

The increment of Cauchy stress tensor dr is also divided into deviatoric and hydrostatic parts

dr ¼ dS þ drm I

ð12Þ

where drm is mean stress increment computed by

drm ¼

1 1 ðdr : IÞ ¼ trðdrÞ 3 3

ð13Þ

The stress and elastic strain relations are described with the Hooke’s law and expressed with the following equations in terms of hydrostatic and deviatoric components respectively

drm ¼ 3Kdem

ð14Þ

and

dS ¼ 2Gdee

ð15Þ

In above expressions, K is the material bulk modulus, and G denotes the elastic shear modulus which is related to Young’s modulus of elasticity, E, through the Poisson’s constant, m. Furthermore, the yield surface f defining the boundary of elastic and plastic deformations in stress space is described by the von Mises yield function

rffiffiffi 2 f ¼ kS  ak  r0 ¼ 0 3

ð16Þ

where r0 is the material yield stress. The yield surface is assumed to be a non-deforming hypersphere, which is allowed to translate and expand during plastic deformation [18–20]. The increment of plastic strain tensor is collinear with the exterior normal to the yield surface at the current stress and remains orthogonal during plastic deformation [18–20]

dep ¼

1 ðdS : nÞn h

ð17Þ

In above expression, h is the plastic hardening modulus function and n is defined as the unit normal tensor denoting the yield surface normal at the current stress point

n¼

Sa kS  ak

ð18Þ

The consistency condition ensures yield condition hold during plastic deformation and is expressed mathematically with the following equation

df ¼ 0

ð19Þ

The size of yield surface will be constant for a cyclically stabilized material stress–strain response, and using the von Mises yield function, the enforcement of consistency condition leads to the following equality during an increment of loading

dS : n  da : n ¼ 0

ð20Þ

The plastic hardening modulus function h may be further expressed with

h¼

dS : n dp

ð21Þ

or using Eq. (20)

h¼

da : n dp

ð22Þ

The total backstress, a, will be decomposed into m number of backstress parts a(i) whose evaluation follows the nonlinear kinematic hardening rule developed by Jiang and Sehitoglu [25]. By substituting Eq. (6) into Eq. (22) and rearranging, the plastic modulus function is obtained in the following explicit form

h¼

m X i¼1

0

1 !vðiÞ þ1 ðiÞ k a k cðiÞ r ðiÞ @1  LðiÞ : nA ði ¼ 1; 2; . . . ; mÞ r ðiÞ

ð23Þ

In Eq. (23), c(i), r(i) and v(i) are three scalar parameters for each backstress part, and L(i) is the unit tensor for the ith backstress component defined previously in Eq. (7). The three plasticity parameters are divided into two groups due to their different functions in modeling the cyclic material response. Considering a balanced proportional loading tests in which the mean strain over a single loading cycle is zero, the material stress response will stabilize eventually after a number of loading cycles. Since no ratcheting is expected for balanced loading conditions experimentally, the mean stress will also be zero over a stabilized deformation cycle. Therefore, the stabilized material deformations is described solely with the plasticity model parameters c(i) and r(i). Consequently, the

845

M. Firat / Materials and Design 34 (2012) 842–852

material cyclic stress–strain curves determined with uniaxial tension–compression tests can be employed directly by integrating the backstress evaluation Eq. (6) along the uniaxial test loading path [25]. The pre-selection of ‘‘m’’ number of stress and strain  pairs raðiÞ ; epaðiÞ from the stabilized cyclic stress–strain curve is the essential step in the application of this parameter identification technique [27]

cðiÞ ¼

rffiffiffi 2 1 3 epaðiÞ

r ðiÞ ¼

2 HðiÞ  Hðiþ1Þ 3 cðiÞ

ði ¼ 1; 2; . . . ; mÞ ði ¼ 1; 2; . . . ; mÞ

ð24Þ ð25Þ

In the above equations, epaðiÞ is the plastic strain amplitude at point (i) on the cyclic stress–strain curve and H(i) designates the slope between points (i  1) and (i) calculated with

HðiÞ ¼

raðiÞ  raði1Þ ði ¼ 1; 2; . . . ; mÞ epaðiÞ  epaði1Þ

ð26Þ

The stress and plastic strain amplitudes take the following initial values and the hardening slope at the last point is zero

rað0Þ ¼ r0 ¼; eað0Þ ¼ 0; Hðmþ1Þ ¼ 0

ð27Þ

There is, however, no direct method for the determination of plasticity model parameters v(i) under a general non-proportional loading since the three model parameters are coupled in the plasticity calculations. Instead, a trial-and-error technique controlling the rate of ratcheting was recommended in the determination of ratcheting exponents v(i) for a given cyclic stress and strain history. For this identification technique, Jiang and Sehitoglu [25] noted the need for additional experimental data under unbalanced loading tests in which the mean stress over a loading cycle is non-zero and ratcheting (cyclic strain accumulation) occurs in one of the strain component. The calculated stress–strain values can be fitted to the experimental response considering a pre-determined plasticity parameter set and their initial values may be adjusted to minimize the difference in the calculation-experiment cyclic material response. While ratcheting under load-controlled conditions are important in stress–strain response of metallic components, the rate of ratcheting can be assumed to have a minimal influence on the fatigue damage prediction [11,32], and therefore v(i) parameters for all backstress components may be taken as constants. Moreover, they may be assumed to be zero for cases, where no experimental data for unbalanced loading of smooth specimens is available and a fast ratcheting may be expected in the plasticity calculation. 3. Stress computation The cyclic plasticity model presented in the previous section is implemented into Ansys FE program through the user material programmable features. Ansys program uses a FE formulation based on the principle of virtual displacement that leads to a set of nonlinear equations in displacement field [31]. Accordingly the linearization of nonlinear FE equations is essential for the numerical solution with the Newton–Raphson algorithm and this process leads to both stresses and stress differentials in the resulting FE algebraic equations [31]. As a result, a FE implementation task consists of calculation of both stress tensor for a given increment of total strain tensor and its derivative with respect to the same incremental strain tensor. In order to compute the stress tensor for a given strain increment determined for a single FE integration point from the global equilibrium iteration, the ordinary differential equations describing the plasticity stress–strain relations are discretized in time using Backward Euler method [31,33]. Next, following the solution method proposed by Chaboche and Cailletaud [34], a nonlinear scalar function is derived to describe the

convergence condition for increment of accumulated plastic strain during the time stepping. The nonlinear scalar equation is solved iteratively by a successive substitution and updating the total backstress and yield function during iterations. Once the convergence to increment of accumulated plastic strain for a given time step is found, the stress and strain tensors at the end of time step are updated. 3.1. Backward Euler discretization and implicit time integration A time interval, Dt from state n to n + 1 will be considered, and it is assumed that all state variables at state n are known and the increment of total strain tensor at time point n + 1 is given. Then, the discretized set of equations employed in the local integration can be summarized as follows.

enþ1 ¼ en þ Denþ1 enþ1 ¼ eenþ1 þ epnþ1 þ ðem Þnþ1 I

ð29Þ

epnþ1 ¼ epn þ Depnþ1

ð30Þ

ðrm Þnþ1 ¼ 3Kðem Þnþ1 rnþ1 ¼ Snþ1 þ ðrm Þnþ1 I

ð32Þ

ð28Þ

ð31Þ

Snþ1 ¼ Sn þ DSnþ1   DSnþ1 ¼ 2G Denþ1  Depnþ1

ð34Þ

Depnþ1

¼ Dpnþ1 nnþ1

ð35Þ

Snþ1  anþ1 ¼ qffiffi

ð36Þ

nnþ1

ro

2 3

anþ1 ¼

m X

aðiÞ nþ1 ¼ an þ Danþ1 ¼ an þ

i¼1

ðiÞ anþ1

ð33Þ

0

m X

ðiÞ Danþ1

ð37Þ

i¼1

1  ðiÞþ1 0  ðiÞ 1v a   nþ1 B ðiÞ C ¼ anðiÞ þ cðiÞ r ðiÞ @nnþ1  @ ðiÞ A Lnþ1 ADpnþ1 r

aðiÞ

ðiÞ nþ1  Lnþ1 ¼   ðiÞ  anþ1 

rffiffiffi 2 ro ¼ 0 fnþ1 ¼ kSnþ1  anþ1 k  3 pnþ1 ¼ pn þ Dpnþ1

ð38Þ

ð39Þ

ð40Þ ð41Þ

In Eqs. (28)–(41), the subscripts n and n + 1 indicates the values of the corresponding quantity at the respective time points, and D represents an increment between time points n and n + 1. In the first step, an elastic predictor step is held, and then a plastic corrector step is carried out if the trial elastic stress state is outside the yield surface [33]. A trial stress, Strial nþ1 is calculated assuming an elastic response for step n + 1, and therefore,

Strial nþ1 ¼ Sn þ 2GDenþ1

ð42Þ

and using the trial stress, the trial yield function is computed, trial fnþ1

  rffiffiffi 2  trial  ¼ Snþ1  an   r0 3

ð43Þ

trial If the condition fnþ1 < 0 holds, the material response is elastic and the trial stress is equal to the stress tensor at step n + 1. Consequently, the increments of equivalent plastic strain, plastic strain and backstress tensors at step n + 1 are set zero, and elastic strain trial and stress tensors are updated. If fnþ1 > 0 inequality holds, the state at step n + 1 is outside of the yield surface at step n, and the deformation is elasto-plastic. Then a plastic corrector step is required to determine the state variables at step n + 1, and the unknown deviatoric stress tensor at step (n + 1) is given by

846

M. Firat / Materials and Design 34 (2012) 842–852

p Snþ1 ¼ Strial nþ1  2GDenþ1

ð44Þ

and the relative stress at step n + 1, Bn+1, becomes

Bnþ1 ¼ Snþ1  anþ1 ¼ Strial nþ1 

2GDepnþ1

 anþ1

ð45Þ

Using Eqs. (37)–(39), the increment of backstress component at step n + 1 is rewritten as follows

ðiÞ anþ1 ¼ anðiÞ þ cðiÞ r ðiÞ nnþ1 Dpnþ1  cðiÞ

  ðiÞ  ðiÞ v anþ1  ðr ðiÞ Þv

ðiÞ

ðiÞ anþ1 Dpnþ1

ffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  Pm ðiÞ ðiÞ   trial Pm ðiÞ ðiÞ  Strial  P a S  P a : i¼1 nþ1 n i¼1 nþ1 n nþ1 nþ1 qffiffi 1¼0   P ðiÞ m 2 r þ 2G þ i¼1 Pnþ1 cðiÞ rðiÞ Dpnþ1 3 o

ð55Þ

or rearranging

Dpnþ1

  qffiffi  trial Pm ðiÞ ðiÞ  Snþ1  i¼1 Pnþ1 an   23ro ¼ P ðiÞ ðiÞ ðiÞ 2G þ m i¼1 Pnþ1 c r

ð56Þ

ð46Þ 3.2. Iterative solution scheme

ðiÞ nþ1

and rearranging for a

0

1  ðiÞ 0  ðiÞ 1v anþ1  ðiÞ B ðiÞ ðiÞ anþ1 @1 þ cðiÞ @ ðiÞ A Dpnþ1 C A ¼ aðiÞ n þ c r nnþ1 Dpnþ1 r

ð47Þ

ðiÞ

At this point, a non-dimensional factor, Pnþ1 , for the ith backstress component is defined as below

The expression derived for the increment of accumulated plastic strain at step n + 1, Dpn+1, is a non-linear equation in that all the ðiÞ non-dimensional factors, Pnþ1 , for the m number of backstress ðiÞ components, Pnþ1 , are dependent on the increment of accumulated plastic strain at step n + 1 itself. Therefore an iterative solution approach is followed as described in the following steps, Initialization phase

1 1  ðiÞ 0  ðiÞ 1v anþ1  B C @1 þ cðiÞ @ ðiÞ A Dpnþ1 A

ðiÞ Pnþ1 ¼0

ð48Þ Step 0: Assume an initial, 0th, iteration value for Dpn+1 as zero, ðiÞ and therefore 1 for all Pnþ1 using Eq. (48).

r

Iteration and substitution phase ðiÞ

and the backstress component, anþ1 , at step n + 1 is rewritten



ðiÞ ðiÞ anþ1 ¼ Pnþ1 anðiÞ þ cðiÞ r ðiÞ nnþ1 Dpnþ1



ð49Þ

The relative stress may also be rewritten as p Bnþ1 ¼ Snþ1  anþ1 ¼ Strial nþ1  2GDenþ1 

m X

ðiÞ anþ1

ð50Þ

i¼1

or expanding for the ith backstress component p Snþ1  anþ1 ¼ Strial nþ1  2GDenþ1 m   X ðiÞ ðiÞ ðiÞ  Pnþ1 aðiÞ n þ c r nnþ1 Dpnþ1

ð51Þ

i¼1

The above equation may be rewritten first by inserting the expression for Depnþ1

Snþ1  anþ1 ¼ Strial nþ1  2GDpnþ1 nnþ1 

m X

ðiÞ PðiÞ nþ1 an

i¼1



m X

ðiÞ Pnþ1 cðiÞ r ðiÞ nnþ1 Dpnþ1

ð52Þ

Step 1: Then calculate the first iteration value for Dpn+1 using Eq. (56). Step 2: Substitute the first iteration value for Dpn+1 into the expressions for relative stress at step n + 1, Sn+1  an+1 given in Eq. (54) and for increment of deviatoric plastic strain at step n + 1, Depnþ1 in Eq. (35) using nn+1 in Eq. (36). Step 3: Use the first iteration values for Depnþ1 and nn+1 to update ðiÞ all Pnþ1 in Eq. (48). Step 4: Finally the second iteration value for Dpn+1 is calculated again using defined scalar equation given in (35). Step 5: Repeat the calculations over the steps 1–4 by the direct substitution of former iteration values to calculate the current iteration value for Dpn+1, and continue this direct substitution and iteration loop until the convergence condition to the solution of increment of accumulated plastic strain at step n + 1 is achieved (Eq. (57)). Once iterations stopped, all state variables (stress and strain tensors together with the backstress components) are updated for increment n + 1.

   count1  Dpiteration   nþ1  6 convergence tolerance 1  count   Dpiteration nþ1

ð57Þ

i¼1

and also for nn+1

3.3. Implementation into Ansys

Snþ1  anþ1 ¼ Strial nþ1  2GDpnþ1 

m X

Snþ1  anþ1  kSnþ1  anþ1 k

ðiÞ Pnþ1 cðiÞ r ðiÞ Dpnþ1

i¼1

m X

ðiÞ PðiÞ nþ1 an

i¼1

Snþ1  anþ1 kSnþ1  anþ1 k

ð53Þ

The above expression may be grouped for relative stress, Sn+1  an+1, as follows

qffiffi

Snþ1  anþ1

 Pm ðiÞ ðiÞ  ro Strial nþ1  i¼1 Pnþ1 an ¼ qffiffi   P 2 ðiÞ ðiÞ Dp r þ 2G þ mi¼1 PðiÞ nþ1 nþ1 c r 3 o 2 3

ð54Þ

The above form of relative stress can be substituted into the yield condition at step n + 1, fn+1, in Eq. (40) in order to define a scalar function for the increment of accumulated plastic strain at step n + 1, Dpn+1, as below

The second essential task in the FE implementation of the presented plasticity model is the calculation of the stress differential, C Lep , which is the derivative of stress tensor with respect to the increment of total strain tensor

C Lep ¼

@ðrnþ1 Þ @ðDenþ1 Þ

ð58Þ

The preferred approach in the determination of C Lep is the derivation of the consistent tangent in order to preserve the quadratic convergence property of Newton–Raphson algorithm in the numerical solution of FE equations [31,33]. The consistent tangent operators for the original A–F model and as well as the models of Chaboche and Onno–Wang have been obtained in the previous studies in the literature and also extended for multi-component models including isotopic hardening [34–38]. A non-symmetric consistent tangent is a common feature for A–F type plasticity

M. Firat / Materials and Design 34 (2012) 842–852

models due to the backstress term in the nonlinear kinematic hardening rules. Although the stability of the numerical solution is granted due to the implicit solution [33], the accuracy of stress computation is influenced significant by the magnitude of loading increment [34,35,37], and it is recommended to reduce the strain increment per step for accuracy contemplation especially for non-proportional loads [36,38]. This situation, however, reduces the benefits of using consistent tangent operator in the elastoplastic FE analysis for the purpose of fatigue damage modeling. This is because the non-proportionality of the loading histories is a typical characteristic of fatigue analysis with nonlinear FE analyses, and consequently, a high number loading increments is usually essential to describe the loading shape fittingly. This condition effectively leads to the use of small incremental strains per load step and makes the use of consistent tangent unfavorable due additional computational costs associated with the calculation of an inverse of a fourth order tensor and the use of non-symmetric matrix in each iteration. Therefore, the classical approach based on the continuum tangent Cep obtained directly from the presented cyclic plasticity model is employed as a substitute for C Lep in Ansys implementation

4

1 4G2 C ep ¼ 2G I  ðI  IÞ þ KðI  IÞ  ðnnþ1  nnþ1 Þ 3 Z nþ1

ð59Þ

4

where I is the fourth order identity tensor and Z is given by

Z nþ1 ¼ 2G þ

m X i

8 9  ðiÞ 0  ðiÞ 1v > > < = anþ1  ðiÞ cðiÞ r ðiÞ  @ ðiÞ A anþ1 : nnþ1 > > r : ;

ð60Þ

847

Both the stress integration algorithm and the continuum tangent matrix are coded into the Ansys Userpl.f subroutine using Fortran programming language and a custom executable is generated using the Intel Fortran Compiler 9.1 [39] in the windows environment.

4. Cyclic test simulations of a notched shaft The presented plasticity model is employed in the FE analysis of a circumferentially notched round bar under combined axial and torsion loading cycles, and the computed notch root strain histories are compared with experimental data collected in biaxial tests performed by Barkey [40]. The notched specimen is a shaft – like component with a net section diameter 25.4 mm and a circumferential notch of radius 12.7 mm (Fig. 1a). Specimens made of normalized 1070 steel stock were machined to the geometry and heat treated to alleviate the manufacturing stresses. The Young’s modulus and Poisson’s constant were also determined to be 210 GPa and 0.3. The stabilized cyclic stress– strain curve of the material was obtained by curve fitting the true stress–strain data from an incremental step testing of a single smooth specimen. Cyclic strength coefficient and cyclic hardening exponent were given as 1736 MPa and 0.199 respectively, and the initial yield stress is determined as 250 MPa. During cyclic loading tests, notch root strain data were measured with a three-element strain gage rosette connected to a digital acquisition unit and converted into Cartesian tensor components (Fig. 1b). The details of testing procedures and experimental work may be found in [40].

(a)

(b)

Fig. 1. (a) Notched round shaft geometry and dimensions (mm) and (b) schematics showing round bar test loads and notch root strain state.

848

M. Firat / Materials and Design 34 (2012) 842–852

Fig. 2. Loading paths showing nominal shear stress in vertical axis and nominal axial stress in horizontal axis in Barkey’s tests [40].

From Barkey’s experiments, six cyclic loading tests are chosen for FE analysis in this study, and the corresponding biaxial loading paths in terms of nominal axial and shear stresses are shown schematically in Fig. 2. The amplitudes of nominal stresses controlled in the test rig are also listed in Table 1. In FE analysis of each test conditions, the time-histories of axial force and torsion moments are calculated by generating identical loading paths, and the corresponding amplitudes, Fa and Ta, are determined using the definitions of nominal axial and shear stresses, Sa and St.

F a ¼ Sa  Anet J T a ¼ St  net r net

Eight-noded solid hexahedral elements (Ansys solid185 element) are used in the 3-D mesh that is generated by sweeping the 2-D pattern-mesh about the axis of symmetry through a uniform sweep angle of 6°. In order to achieve a computationally-efficient FE mesh for the cyclic loading simulation, a mesh convergence study is conducted in that the number of elements over the notch root section is decreased while the smallest element is always located at the notch root. The 2-D pattern-meshes in convergence study

ð61Þ ð62Þ

Anet and rnet are the cross section area and radius of the notch root respectively, and Jnet denotes the polar moment of inertia of the same cross section. The calculated time histories are used as the time-dependent force boundary condition applied to the nodes over the area at the specimen smooth end (Fig. 1a). The half of the specimen is included in the FE model due to the geometric symmetry as well as mechanical properties and boundary conditions.

Table 1 The nominal stress amplitudes and corresponding loading paths in FE simulations. Test case

Nominal axial stress Sa (MPa)

Nominal shear stress St (MPa)

Loading path

1 2 3 4 5 6

139 296 296 296 296 296

87.5 193 193 193 193 193

Path Path Path Path Path Path

1 1 2 3 4 5

849

M. Firat / Materials and Design 34 (2012) 842–852

Mesh 1

Mesh 2

Mesh 3

Fig. 3. Notch zone mesh layouts in mesh-convergence study.

Table 2 The computed stress concentration factors with three FE models. Pattern mesh

Number of elements

Number of nodes

Axial stress concentration factor, Ka

Shear stress concentration factor, Kt

1 2 3

16,140 10,245 4020

17,418 11,670 4571

1.438 1.421 1.413

1.132 1.137 1.145

are shown in Fig. 3. Linear elastic stress analyses with unit axial force and torsion moments were conducted to calculate the axial and shear stress concentration factors. The calculated concentration factors with all FE meshes are also provided in Table 2. Barkey [40] determined axial and shear stress concentration factors with strain gage measurements and reported as 1.45 and 1.17 respectively. An investigation of the calculated concentration factors indicates that all FE meshes for the specimen are acceptable and only slight changes on the order of 1% or 2% are determined as the number of elements decreases. Considering the computational gain for a nonlinear time-history analysis, the specimen FE model generated from the pattern mesh 3 is used in the cyclic deformation analyses (Fig. 4). The total backstress is composed of 10 parts (m = 10) and the plasticity model parameters, c(i) and r(i), for each backstress component are calculated using the material cyclic stress–strain curve. Accordingly, the cyclic curve is digitized with 10 stress– strain pairs (r(i), e(i)) and a total of 20 plasticity model parameters were calculated with Eqs. (24)–(27). For all backstress parts, a constant ratcheting exponent is taken (v(i) = 1 with i = 1–10). The same set of model parameters are employed in all test simulations, and first ten loading cycles were calculated using an equally-spaced 80 increments per loading cycle. The only exception is the test case 1 for which the FE calculation was done for a single loading cycle due to the elastic deformations. The convergence condition for the global Newton–Raphson iterations during load stepping is based on residual force criteria of maximum value 0.1%[31]. All FE analyses were performed using a four-processor Windows-64XP PC.

4.1. A performance assessment FE-computed shear strain eyz and axial strain ezz histories at the notch root of the specimen are compared with measurements. In all simulated test cases, except for path 5 with a proportional axial force-torsion variation only, a preloading step is held, where the axial force and torsion increase proportionally to the load point, where the equivalent stress becomes maximum, and then the equivalent stress fluctuates between the maximum value and minimum values for four-times in a single cycle. Consequently, the equivalent stress at the end of the preloading determines the notch root deformation to be either elastic or elasto-plastic for continuation of loading cycles. Considering the first loading case along a box-type axial force-torsion loading path, the FE-calculated strain histories are purely elastic and the shear strain – axial strain plot is a perfect box (Fig. 5). The measured axial – shear strain plot is also of the box-like shape and the loops follow the FE-calculated (elastic) strain plot very closely. The relatively small rounding at each corner is related to the controller error in the actual test [40]. The high correlation between experimental and computed strain plots proved the accuracy of material parameters and FE modeling used in numerical analysis. The notch root deformation response changed from elastic to elasto-plastic regime with increased nominal stress values in test case 2, and the experimental strain plot becomes distorted leading a clockwise rotation of the box-shape. This shape change indicates the elasto-plastic coupling of shear and axial strain components due to the plastic deformation occurred at the notch root. An exploration of computed strain histories for this test case also shows the same coupling effect fourtimes in each loading cycle. Fig. 6 compares graphically measured and computed strain variations. There is also a steadily shift in the experimental shear strain – axial strain loop as the number of loading cycles increases, and an asymmetry developed gradually (Fig. 6). This experimental trend has been also simulated successfully, and the computed strain plot shows a fairly similar asymmetry with box-type loading up to 10 cycles. But the rate of

Fig. 4. The three dimensional FE model of the notched specimen.

M. Firat / Materials and Design 34 (2012) 842–852 0.1500

0.48

0.0750

0.32

0.0000

0.16

shear strain (%)

shear strain (%)

850

-0.0750

-0.1500 -0.1500

-0.1000

-0.0500

0.0000

0.0500

0.1000

-0.16

0.1500

axial strain (%) experiment

0.00

-0.32

model prediction

Fig. 5. Measured and calculated notch root strains for test program 1.

-0.48 -0.36

-0.18

0.00

0.18

0.36

axial strain (%) experiment

model prediction

0.6

Fig. 8. Measured and calculated notch root strains for test 4.

Cyclic Shift 0.4

0.32 0

shear strain (%)

shear strain (%)

0.48 0.2

-0.2

Coupling in shear and axial strains

-0.4

-0.6 -0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.16

0.00

-0.16

axial strain (%) experiment

model prediction

-0.32

Fig. 6. The shift in the shear strain – axial strain response with increasing number of cycles under box type of loading in test 2.

-0.36

-0.18

0.00

0.18

0.36

axial strain (%) experiment

model prediction

Fig. 9. Measured and calculated notch root strains for test 5. 0.56

0.40

Forward strain path

shear strain (%)

0.24

0.08

-0.08

-0.24

Reverse strain path

-0.40

-0.56 -0.36

-0.18

0.00

0.18

axial strain (%) experiment

model prediction

Fig. 7. Measured and calculated notch root strains for test 3.

0.36

shift in first loading cycle is relatively high compared to the experimentally observed case. The appearance of this gradual shift in the notch root deformation response is related to the unbalanced cyclic loading applied to the specimen along a box-type axial-torsion loading history. The kinematic hardening plasticity model implemented in the presented simulations has the proper attribute to simulate cyclic deformations under unbalanced loading conditions [25]. From the FE analysis of zigzag-type loading path (path 3), elasto-plastic cyclic deformations were predicted at the notch root due to the deviation of forward and reverse strain loops (Fig. 7). In the case of elastic notch deformations, the shape of shear strain – axial strain plot during a cycle is still form identical to that of nominal shear – nominal stresses since the time-histories of both nominal stresses were prescribed in the FE modeling for each loading cycle. Experimental data confirmed the FE results and a similar difference along forward and reverse strain paths is observed in the measured shear strain and axial strain plot. Compared with measurements, the computed strain histories are matched qualitatively with the

851

M. Firat / Materials and Design 34 (2012) 842–852 0.06

0.5

shear strain axial strain

0.375

0.05

error in strain (%)

shear strain (%)

0.25 0.125 0 -0.125

0.04

0.03

0.02

-0.25 0.01

-0.375 0

-0.5 -0.3

1

-0.15

0

0.15

0.3

model prediction

3

4

5

6

Loading case

axial strain (%) experiment

2

Fig. 11. The strain calculation error (using the SMSD parameter) as a function of test loading case.

Fig. 10. Measured and calculated notch root strains for test 6.

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1 XN  i i SMSD ¼ e  e computed measured i¼1 N measured strains, and the experimental trends in shear strain – axial strain histories are tracked properly. The measured and computed strain variations along v-type paths (path 4 and path 5) are shown graphically in Figs. 8 and 9. From the measured strain plot, it can be observed that the lines in nominal stress cycles are now converted into two different branches for forward and return paths in measured strain histories. The computed strain histories are in accord with experimental results and the strain ranges in each cycle are computed with a good accuracy. Furthermore, a strain differences with each loading cycle is also simulated qualitatively with presented cyclic plasticity model, and this prediction is confirmed from the measured shear and axial strain histories irrespective of loading path of the test. However, cyclic shift in measured shear strain history is not computed and symmetric strain loops with a smaller axial strain range are computed for both test conditions. Considering the FE-computed strain histories for proportional loading path (path 5), and a nonlinear strain loop is predicted (Fig. 10). The redistribution of stresses is responsible for different forward and reverse branches during cyclic loading and the strain localization at the notch root appeared in the form of a closed, symmetric hysteresis loop in shear strain – axial strain plot. A fairly similar loop is observed in the measured strain plot, but it is not symmetric with respect to origin. The measured axial strain range is, moreover, larger than the FE prediction and the experimental loop was moved toward the compressive side. Furthermore, similar to the previous experimental results for paths 3 and 4, the slight shift of shear strain history was not simulated. This is consistent with the theory of implemented plasticity model and no shift in notch strain history is expected with each deformation cycle under a balanced loading path. In a recent study, Ye et al. [41] simulated a proportional axial force – torsion loading of a notched shaft using the same plasticity model. Their analysis with a different material also produced a nonlinear strain loop, and similarly, no cyclic shift along shear strain direction was calculated. The performance of the presented computational model is evaluated in view of the accuracy of the calculated strain histories and the computer solution times. Considering the measured and FEcomputed strain histories as separate time series, an average error of FE strains is calculated by means of an arithmetic mean of square root of mean of squared difference (SMSD) parameter. SMSD parameter is defined for both shear and axial strain components and is given by

ð63Þ

where N is the number of strain data points in measurements and calculations. The SMSD parameters calculated for each test case separately are shown graphically in Fig. 11. An examination of the computed SMSD parameters indicates that the error in computed strains is a function of the loading path shape and maximum error occurred in the FE analysis of box-type loading test. For all tests with elastic–plastic notch deformations, the average error in FEcomputed shear strains is consistently higher than the error in axial strain and determined less than 0.06% (Fig. 11). An assessment of the computer computational times showed that the FE analysis of the box-type loading test took approximately 75 min computer clock time on the four-processor Intel (Windows-64XP) machine and is determined to be the maximum among all test cases simulated. The average number of global equilibrium iterations per increment is about 6 for this loading path and higher than that for other test conditions. This is an acceptable cost considering the probable accuracy enhancement of existing damage models or to develop new models for the fatigue life predictions. 5. Conclusions In this paper, a cyclic plasticity model suitable for fatigue damage modeling of metallic structures is presented and its FE-implementation is described within the small strain plasticity framework. The model is based on the multi-component additive decomposition of the total backstress and uses the von Mises yield surface whose evaluation in stress space follows the Jiang–Sehitoglu nonlinear kinematic hardening rule. An incremental implicititerative algorithm was employed for the numerical solution of resulting stress–strain equations, and the continuum tangent obtained directly from plasticity model was used in Ansys FE implementation. The developed FE computational model is applied in the cyclic deformation simulation of a circumferentially notched specimen loaded under proportional and nonproportional axial force-torsion loadings. The computed notch root deformations were compared with measured notch root strain histories. An assessment of model predictions showed that all loading tests have been simulated with a good accuracy for both elastic and elastoplastic notch deformations. The computed strain loops were in accord with experimental data and matched qualitatively with measured shear – axial strain histories irrespective of loading path

852

M. Firat / Materials and Design 34 (2012) 842–852

of the test. In proportional balanced torsion-axial loading, the nonlinear shear strain – axial strain loops were also predicted properly. The average error for both notch strain components is a function of the loading path shape, and compared to shear strains, the axial strain errors were relatively smaller. The overall computational cost of FE analyses was also acceptable. Acknowledgments Author thanks Prof. Dr. E. Barkey for providing the technical material and experimental data of notched specimen test programs. Author also thanks Dr. T. Ögüt for supporting the use of Ansys software in this study. References [1] Lee DC, Han CS. CAE (computer aided engineering) driven durability model verification for the automotive structure development. Finite Eleme Anal Des 2009;45:324–32. [2] Conle FA, Chu C. Fatigue analysis and the local stress–strain approach in complex vehicular structures. Int J Fatigue 1997;19(1):S317–23. [3] Firat M, Kocabıcak U. Analytical durability modeling and evaluation – complementary techniques for physical testing of automotive components. Eng Fail Anal 2004;11(4):655–74. [4] Tipton SM, Nelson DV. Advances in multiaxial fatigue life prediction for components with stress concentrations. Int J Fatigue 1997;19(6):503–15. [5] Peterson RE. Stress Concentration Factors. New York: John Wiley & Sons; 1977. [6] Knop M, Jones R, Molent L, Wang C. On the Glinka and Neuber methods for calculating notch tip strains under cyclic load spectra. Int J Fatigue 2000;22:743–55. [7] Neuber H. Theory of stress concentration for shear strained prismatic bodies with arbitrary stress–strain law. J Appl Mech 1961;28:544–50. [8] Molski K, Glinka G. A method of elastic–plastic stress and strain calculation at a notch root. Mater Sci Eng 1981;50:93–100. [9] Buczynski A, Glinka G. An analysis of elasto-plastic strains and stresses in notched bodies subjected to cyclic non-proportional loading paths. In: Carpinteri A, de Freitas M, Spagnoli A, editors. Biaxial/multiaxial fatigue and fracture. Amsterdam: Elsevier; 2003. p. 265–84. [10] Savaidis A, Savaidis G, Zhang Ch. FE fatigue analysis of notched elasto-plastic shaft under multiaxial consisting of constant and cyclic components. Int J Fatigue 2001;23:303–15. [11] Jiang Y, Hertel O, Vormwald M. An experimental evaluation of three critical plane multiaxial fatigue criteria. Int J Fatigue 2007;29:1490–502. [12] Mroz Z. An attempt to describe the behavior of metals under cyclic loads using a more general workhardening model. Acta Mech, Cilt 1969;7(2–3):199–212. [13] Lamda HS, Sidebottom OM. Cyclic plasticity for nonproportional paths. Part I – Cyclic hardening, erasure of memory, and subsequent strain hardening experiments. J Eng Mater Technol, Cilt 1978;100:96–104. [14] Drucker DC, Palgen L. On the stress–strain relations suitable for cyclic and other loadings. J Appl Mech, Cilt 1981;48:479–85. [15] Chaboche JL. Constitutive equations for cyclic plasticity and cyclic viscoplasticity. Int J Plast 1989;5:247–302.

[16] McDowell DL, Socie DF, Lamda SH. Multiaxial nonproportional cyclic deformation. In: Amzallag C, Leis BN, Rabbe P, editors. Low-cycle fatigue and life prediction. ASTM STP, vol. 770. Philadelphia: American Society for Testing and Materials; 1980. p. 500–18. [17] Besseling JF. A theory of elastic, plastic, and creep deformation of an initially isotropic material showing anisotropic strain-hardening, creep recovery, and secondary creep. J Appl Mech 1958;17:529–36. [18] Lemaitre J, Chaboche J-L. Mechanics of solid materials. Cambridge University Press; 1990. [19] Lubliner J. Plasticity theory. New York: Macmillan Publishing Company; 1990. [20] Khan AS, Huang S. Continuum theory of plasticity. New York: John Wiley & Sons Inc.; 1995. [21] Dafalias YF, Popov EP. A model of nonlinearly hardening materials for complex loading. Acta Mech, Cilt 1975;21(2–3):173–92. [22] Krieg RD. A practical two surface plasticity theory. J Appl Mech 1975:641–6. [23] Armstrong PJ, Frederick CO. A mathematical representation of the multiaxial bauschinger effect. Report RD/B/N 731, Central Electricity Generating Board; 1966. [24] Ohno N, Wang JD. Kinematic hardening rules with critical state of dynamic recovery. Part I – Formulation and basic features for ratcheting behavior. Int J Plast, Cilt 1993;9:375–390–650. [25] Jiang Y, Sehitoglu H. Modeling of cyclic ratcheting plasticity. Part I – Development of constitutive relations. J Appl Mech 1996;63:720–5. [26] Mingzhou S, Qiang G, Bing G. Finite element analysis of steel members under cyclic loading. Finite Elem Anal Des 2002;39:43–54. [27] Kocabıcak U, Firat M. A simple approach for multiaxial fatigue damage prediction based on FEM post-processing. Mater Des 2004;25(1):73–82. [28] Liu Y, Liu L, Stratman B, Mahadevan S. Multiaxial fatigue reliability analysis of railroad wheels. Reliab Eng Syst Safety 2008;93:456–67. [29] Choa SK, Yanga YS, Sona KJ, Kimb JY. Fatigue strength in laser welding of the lap joint. Finite Elem Anal Des 2004;40:1059–70. [30] Firat M, Kozan R, Ozsoy M, Mete OH. Numerical modeling and simulation of wheel radial fatigue tests. Eng Failure Anal 2009;16(5):1533–41. [31] Ansys Theory Manual. Canonsburg, PA: Ansys Inc.; 2003. [32] Jiang Y, Sehitoglu H. Cyclic ratcheting of 1070 steel under multiaxial stress states. Int J Plast 1994;10(5):579–608. [33] Simo JC, Hughes TJR. Computational inelasticity. New York: Springer-Verlag Inc.; 1998. [34] Chaboche JL, Cailletaud G. Integration methods for complex plastic constitutive equations. Comput Methods Appl Mech Eng 1996;133:125–55. [35] Doghri I. Fully implicit integration and consistent tangent modulus in elastoplasticity. Int J Numer Methods Eng 1993;36:3915–32. [36] Hartmann S, Haupt P. Stress computation and consistent tangent operator using non-linear kinematic hardening models. Int J Numer Methods Eng 1993;36:3801–14. [37] Kobayashi M, Ohno N. Implementation of cyclic plasticity models based on a general form of kinematic hardening. Int J Numer Methods Eng 2002;53:2217–38. [38] Hartmann S, Luhrs G, Haupth P. An efficient stress algorithm with applications in viscoplasticity and plasticity. Int J Num Meth Eng 1997;40:991–1013. [39] Intel Fortran Compiler 9.1, Programmer’s Manual and Documentation. Intel Corporation; 2006. [40] Barkey ME. Calculation of notch strains under multiaxial nominal loading. Ph.D. Dissertation, University of Illinois at Urbana-Champaign; 1993. [41] Ye D, Hertel O, Vormwald M. A unified expression of elastic–plastic notch stress–strain calculation in bodies subjected to multiaxial cyclic loading. Int J Solids Struct 2008;45:6177–89.