BEM implementation of Robinson's viscoplastic model in creep analysis of metals using a complex variable numerical technique

Advances in Engineering Software 33 (2002) 805–816 www.elsevier.com/locate/advengsoft

D/BEM implementation of Robinson’s viscoplastic model in creep analysis of metals using a complex variable numerical technique C.P. Providakis* Department of Applied Sciences, Technical University of Crete, GR-73100 Chania, Greece Accepted 17 July 2002

Abstract A new boundary element approach for the solution of time-dependent inelastic problems arising in creeping metallic structures subjected to the combined action of high temperature gradient and quasi-static mechanical loading conditions is investigated. The new approach allows the use of complex variable techniques in the boundary element procedure for the evaluation of stress components as derivatives of the displacement integral equations. This methodology makes faster and more accurate the conventional boundary element method. To validate the efficiency of the proposed method in the implementation of Robinson’s viscoplastic model, the results obtained using the present methodology are compared with those obtained by using known analytical and finite element solution for the analysis of a thick-walled internally pressurized cylinder and an experimental cylindrical thrust chamber in plane strain under general thermomechanical loading histories. q 2002 Elsevier Science Ltd. All rights reserved. Keywords: Complex variables; Viscoplasticity; Robinson’s model; Boundary element method; Thermomechanical loading

1. Introduction Recently, remarkable progress has been made in the technological applications of metals and alloys at elevated temperatures and thus, in order to improve the high temperature inelastic behavior of structures, especially those used in nuclear and aerospace industries, some new and more realistic constitutive models (called viscoplastic models) have been proposed. A very comprehensive review of numerous viscoplastic models was developed by Walker [1] and Lindholm et al. [2], Lemaitre and Chaboche [3], Freed and Walker [4] and Saleeb et al. [5]. The main mathematical feature of these newer viscoplastic models is that the non-elastic strain rates can be expressed as functions of the current values of stress, temperature and some certain well defined state variables. The mathematical formulation of these viscoplastic models became very complex since they incorporate constitutive equations, which were highly non-linear and mathematically stiff. Thus, it is concluded that many practical problems, which involve complex geometries as well as complex thermomechanical loading histories such as those arising in hot gas-path components of gas turbine engines * Tel./fax: þ30-821-37437. E-mail address: [email protected] (C.P. Providakis).

and rocket engines, have to be solved in the environment of such numerical solution methodologies which could make these viscoplastic models adaptable for realistic structural and life analyses of this kind of structural components. Robinson’s model [6,7] is one of these viscoplastic models using state variables that aim to represent inelastic behavior of metallic structures more faithfully than is possible with the traditional inelastic models in a variety of applications at elevated temperature thermomechanical loading. Arya [8,9], Arya and Arnold [10] and Arya and Kaufman [11] proposed a finite element methodology (FEM) for the inelastic analysis of structures with material behavior governed by the Robinson’s constitutive viscoplastic model. They used a rate formulation of the governing differential equation of the problem. This implementation was exercised first on several uniaxial problems involving isothermal and non-isothermal cyclic loading and later more complex structural problems were treated. This FEM was found to work efficiently for the selected numerical examples. In Providakis and Kourtakis [12] demonstrated that the domain/boundary element method (D/BEM) is a very powerful method with several potential advantages over the FEM for solving non-linear time-dependent inelastic deformation problems in the presence of high temperature gradients. One of the main advantages is that the number of

0965-9978/02/$ - see front matter q 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 9 6 5 - 9 9 7 8 ( 0 2 ) 0 0 0 4 1 - 8

806

C.P. Providakis / Advances in Engineering Software 33 (2002) 805–816

unknowns in the final algebraic system of equations is proportional to the number of boundary elements in D/BEM as opposed to the total number of nodes in FEM. However, one of the most difficult problem in BEM-related problems was the numerical simulation of domain-based effects (e.g. inertial and interior loading effects, inelastic and thermal strain effects). This is caused by the existence of domain integrals in the formulation, which can only seldom be directly transferred into boundary integrals forms. In high temperature and time-dependent creeping problems, domain integrals arise due to thermal and inelastic term effects. However, some inherent difficulties, such as the strong singularities arising in the derivation of the domain integrals, have stymied the further development of this method. Several approaches have been developed to deal with the removal of the strong singularities arising in the domain integrals such as those referenced in Gao and Davies [13]. The present paper could be considered as an extension of the previous work presented in Providakis and Kourtakis [12] where a boundary element methodology was proposed for the implementation of Robinson’s viscoplastic model to the analysis of time-dependent non-linear inelastic deformations of metallic structural components subjected to high temperature thermomechanical loading histories. The present formulation is based in the same as the above work boundary element implementation of Robinson’s viscoplastic model but now adopts, actually for the first time, the use of a complex variable numerical technique to eliminate the strong singularity involving in the evaluation of internal stresses in this kind of complicated inelastic problems in the existence of high temperature gradients. This is achieved by an appropriate introduction of complex variables into the displacement integral equations. This numerical approach which was proposed by Lyness and Moler [14] was applied to the computation of the derivatives of any complicated function by using a step size in the imaginary part of the variable. The proposed new approach of D/BEM method can be considered numerically as an ‘exact derivative method’ since its numerical accuracy can be guaranteed by taking a sufficiently small step size. Following this concept, the stresses at any internal point could be accurately obtained by using only the numerical derivative of the displacement integral equations which, however, involve only weakly singular integrals. This technique makes faster and more accurate the conventional boundary element algorithm since the complicated and error sensitive methodologies proposed by other researchers for evaluating the internal stresses could be now eliminated. In the present boundary element formulation the time rates of change of displacements and stresses are generated at any time step. The time histories of the quantities of interest are then obtained by using a time integration scheme with automatic time-step control. For the implementation of the proposed boundary element formulation isoparametric linear boundary elements were used. Due to the lack of existence of other available solutions, to compare with the

present D/BEM approach results an analytical solution is developed and discussed herein which is, however, valid only for the case of the axisymmetric analysis of a thickwalled thermoviscoplastic cylinder. Then, a comparison of the two solutions is achieved by developing two different programs, the first one for the present viscoplastic D/BEM and the second one for the analytical axisymmetric approach. In order to avoid programming bias and make the comparison as meaningful as possible both of the programs are created by the same research group. Finally, numerical examples are presented and compared with analytical and FEM solutions in the last section of this paper for plane strain deformation problems of cylinders subjected to mechanical and high temperature loading histories.

2. Robinson’s viscoplastic model derivation A small displacement and a small strain formulation are employed and the rate of the total strain is assumed to consist of 1_eij ; 1_nij the elastic and non-elastic (including plastic, creep, relaxation, etc.) strain rate components and 1_Tij the thermal strain rate component given by the temperature rate field through the expression 1_Tij ¼ aT_ dij with a being the thermal expansion coefficient. The mathematical structure of many of the state variable models of viscoplastic deformation can be summarized by the following equations 1_ij ¼ 1_eij þ 1_nij þ 1_Tij q_ ðkÞ iij

¼

1_nij ¼ hðsij ; qðkÞ ij ; TÞ

1_kk ¼ 0 ð1Þ

f ðsiij ; qðkÞ ij ; TÞ

where sij is the stress tensor, T is the temperature and qðkÞ ij are the state variables. Robinson’s model, used in this paper, was found to be very well suited to non-isothermal structural analysis. The inelastic strain rate in this model is governed by the following flow law ( f ðFÞSij F . 0 and Sij Sij . 0 n 1_ij ¼ 0 F # 0 or F . 0 and Sij Sij # 0 ð2Þ Fn f ðFÞ ¼ 2m where the effective stress Sij is the difference between the deviatoric stress Sij and the thermal deviatoric stress aij ; m is a material constant and the variable F n is given later. The state variable in the Robinson’s model [6] is the deviatoric stress or backstress aij which accounts for the kinematic hardening. According to the Robinson’s model, the growth law governing the evolution of internal state variable aij is given by

a_ ij ¼ hðakl Þ1_ij 2 rðakl Þaij

ð3Þ

The first term denotes the hardening process with

C.P. Providakis / Advances in Engineering Software 33 (2002) 805–816

accumulated deformation and the second term, a recovery or softening process proceeding with time. Under the steadystate conditions, these two mechanisms balance each other and consequently, a_ ij ¼ 0: The expressions for the functions hðakl Þ and rðakl Þ are given by h¼

h¼

2m H ; Gb 2m H Gb0

r¼

;

r¼

RGm2b pffiffiffi ; G . 0 and Sij aij . 0 I2 b RGm2 0

1 2

ð4Þ

ð

þ

½Uij ðxP ; xQ Þt_j ðxQ Þ 2 Tij ðxP ; xQ Þ_uj ðxQ ÞdGðxQ Þ

ð V

Sij Sij ;

F¼

J2 21 k2

I2 ¼

1 2

Sij ¼ sij 2

1 3

skk dij

ð5Þ

ð6Þ

and the parameters m, H, m, b, and G0 can be considered as temperature independent constants. The temperature-dependent constants m and R can be obtained by using the equations

m ¼ m expð2u1 Þ

ð7Þ

and R ¼ 9:0 £ 1028 expð2u2 Þ   1 1 2 u1 ¼ ð23:8u 2 2635Þ 811 u   1 1 2 u2 ¼ 40; 000 811 u

ð8Þ

where u is the absolute temperature expressed in Kelvin (K). The scalar k accounts for the isotropic hardening (or softening) and can be considered as the threshold (Bingham) stress. In the present derivation, the scalar k is assumed to be temperature independent.

3. Boundary element formulation The Navier equation for the displacement rates for plane stress, in the presence of non-elastic strain are 1þn u_ 1 2 n k;ki 2n n 2ð1 þ nÞ _ F_ i þ 21_nij;j þ 1_ þ ðaTÞ;i 12n G 1 2 n kk;i

_ ½Xjki ðxP ; xq Þ1_njk ðqÞ þ X^ jki ðxP ; xq Þdjk aTðqÞ

dVðxq Þ

aij aij ;

where the deviatoric stress, Sij ; is defined as

¼2

¼

G

pffiffiffi ; G # G0 or Sij aij # 0 I2

pffiffiffi I2 G¼ k

u_ i;kk þ

rate boundary conditions must be prescribed. For the plane strain case ð1zz ¼ 0Þ the integral representation of the solution for a source point P with coordinates xP ; and a field point Q with coordinates xQ on the boundary of the body (with F_ i ¼ 0) has the form ðdij 2 Cij ðxp ÞÞ_uj ðxP Þ

where J2 ¼

807

ð9Þ

where Fi is the prescribed body force per unit volume, G, n and a are the shear modulus, Poisson’s ratio and coefficient of linear thermal expansion, respectively, ui is the displacement vector. Suitable traction and displacement

ð10Þ

where dij is the Kronecker delta, P, Q are boundary points, q is interior point, G and V are the boundary and the surface of the body, respectively. The kernels Uij ; Tij ; Xjki and X^ jki are known singular solutions due to a point load in an infinite elastic solid in plane strain [15] and defined as functions of the coordinates xP and xQ at the source point P and field point Q, respectively. The traction and displacement rates are denoted by t_ and u_ ; respectively. The coefficients of Cij ðxp Þ are known functions of the included angle at the boundary corner at P, the angle between the bisector of the corner angle and the x-axis. Eq. (10) is a system of integral equations for the unknown traction and displacements rates in terms of their prescribed values on the boundary, the nonelastic strain rates and the temperature rates profile. The displacement rates for a point p in the interior of the body can be obtained as a special case of Eq. (10) by taking the term dij 2 Cij ðxp Þ being equal to 1. Now, the coordinate of the source point p should be converted to a complex variable by adding a small imaginary part ibm to the mth coordinate with a step size bm become sufficiently small in the order of machine zero (up to 10235). Then, Eq. (10) can be transformed to the same equation but with xp replaced by xp þ ibm and by simply declaring the corresponded variables in the kernels functions as complex. This easily suggests that the conventional computer algorithm [12] of the thermoviscoplastic analysis of creeping metallic structures remains exactly the same as the present one except complex variable operations have to be involved in its formulation. The stress rates can be obtained by using Hooke’s law in terms of the derivatives of the displacements rates at the source point P through the equations ! ›u_ j ›u_ i Gn ›u_ k s_ ij ðxp Þ ¼ G ðx Þ þ ðx Þ þ 2 ðx Þd ›xj p ›x i p 1 2 2n ›xk p ij     1þn _ p Þdij 22 G1_nij ðxp Þ þ G aTðx ð11Þ 1 2 2n where G is the shear modulus. Now, once the evaluation of the rates of the displacement ui ðxp þ ibÞ is done, the derivatives of the displacement rates

808

C.P. Providakis / Advances in Engineering Software 33 (2002) 805–816

can be obtained by using the following equation given from the procedure of Ref. [14] for the numerical evaluation of the derivative of a function of complex variables

›u_ i 1 ðxp Þ ¼ Im½_uðxpm þ ibm Þ bm ›x m

ð12Þ

It is clear that by employing the above conversion into Eq. (11), the stress rates will deal with displacement kernels Uij ðxpm þ ibm ; xQ Þ and traction kernels Tij ðxpm þ ibm ; xQ Þ which display a weak singularity of order log r and 1/r, respectively. Thus, the resulted stress integral equations is no more singular than that of displacement integral equations and, actually, they present exactly the same order of singularity. The present procedure is seems to be considered numerically as an exact derivative method since its numerical accuracy can be guaranteed by taking a sufficiently small step size b in the order of machine zero (up to 10235) and, in contrast to other numerical methods, it seems that by using the present D/BEM methodology no subtractive cancellation errors occur and the results tend to be independent of the step size bm as the step size becomes sufficiently small (up to 10235).

The non-elastic strain rates terms at an arbitrary point within the element Vi can be defined by using the equations b

1_njk ¼ N b ðz1 ; z2 ÞE_ njk

X~ j ¼ N a ðzÞX~ aj

u_ j ¼ N a ðzÞU_ aj

t_j ¼ N a ðzÞP_ aj

where E_ njk and T_ b represent the vectors of the increments of the non-elastic strain and temperature rates terms at an arbitrary point j inside the interior element Vi : The isoparametric boundary element representation of the integral Eq. (10) utilizing the function expansions (13) – (15) can be written as ðdij 2 Cij ÞU_ j ðxj Þ ¼

L ð X

where N a ðzÞ is a set of polynomial shape functions defined on element Gl; z is an intrinsic coordinate on Vl which varies between 2 1 and þ 1 and the superscript a is summed from _ aj P_ aj are 1 to n which is the number of nodes on Gl. X~ aj ; U vectors containing the nodal values of coordinates, displacement rates and the boundary traction rates, respectively. For the discretization of the non-elastic strain and temperature rates integrals, it is considered that the coordinates of an arbitrary interior point j within the interior element Vi can be calculated by the equation x~ j ¼ N b ðz1 ; z2 Þ~xbj

Gl

l¼1

2

L ð X

N

ð14Þ

where x~ j is the vector that contains the cartesian coordinates of an arbitrary interior point j within the element Vi ; N b ðz1 ; z2 Þ is a set of polynomial shape functions and X~ bj is the vector of the cartesian coordinates related to the nodal point of the element Vi ; z1 and z2 are intrinsic coordinates on any interior element Vi and the superscript b is summed from 1 to m  which is the number of nodes on element Vi :

a

Gl

l¼1

þ

 ~ zÞdGðXÞÞ ~ P_ aj N a ðzÞUijp ðxj ; Xð

N ð X



~ zÞÞdGðXÞ ~ ðzÞTijp ðxj ; Xð

_ aj U b

Vn

n¼1

p N b ðz1 ; z2 ÞXjki ðxj ; x~ ðz1 ; z2 ÞÞdVð~xÞE_ njk

N  ð X þ a

N Vn

n¼1

ð13Þ

ð15Þ

b

4. Matrix formulation and solution strategy Eqs. (10) and (11) can be expressed in numerical form by discretizing the boundary and the interior into a number of two noded linear elements and four noded linear isoparametric interior elements, respectively. The discretization of boundary integrals is performed by using the coordinates, the displacement and traction rates of an arbitrary point within the element Gl which can be calculated by the following equations

T_ ¼ N b ðz1 ; z2 ÞT_ b

b

ðz1 ; z2 ÞX^pjki ðxj ; x~ ðz1 ; z2 ÞÞdVð~xÞ



 T_ b ð16Þ

The stress rates can be evaluated in a similar way by using Eq. (11) while the derivatives can be obtained by using the equations

›u_ i ðx Þ ›x j j ¼

L ð X l¼1

2

L ð X l¼1

þ

N X

ð

n¼1

þ

N X n¼1

a

! 1 p ~ zÞÞdGðXÞ ~ P_ aj N ðzÞ Im½Uij ðxjj þ ibj ; Xð bj Gl a

! 1 p _ aj ~ ~ N ðzÞ Im½Tij ðxjj þ ibj ; XðzÞÞdGðXÞ U bj Gl a

b 1 p N ðz1 ; z2 Þ Im½Xjki ðxjj þ ibj ; x~ ðz1 ; z2 ÞÞdVð~xÞE_ njk bj Vn

!

b

ð Vn

N b ðz1 ; z2 Þ

! 1 Im½X pjki ðxjj þ ibj ; x~ ðz1 ; z2 ÞÞdVð~xÞ T_ b bj

ð17Þ where L is the number of boundary elements and N the number of interior surface elements; Gl is the lth boundary elements ðG ¼ SGl Þ and Vn is the nth surface element ðV ¼ p SVn Þ; Uijp ; Tijp ; Xjki and X pjki are the corresponding tensors for the boundary integrals, inelastic and thermal effect tensors for surface integrals of Eqs. (10) and (17). The evaluation of the coefficients of the matrices in Eqs. (16) and (17) needs a number of complicated integration procedures. Since analytical integration of the integrals in these equations is not possible in general, the Gaussian quadrature technique was used. For the boundary singular cases which occur

C.P. Providakis / Advances in Engineering Software 33 (2002) 805–816

when the field and the source point are situated over the same boundary element special approaches were employed as in Providakis approach [16] in plate bending problems. The evaluation of the singular surface inelastic and temperature rates profile integrals was done by subdividing each fournoded linear isoparametric interior element into a number of triangular subelements with their appexes located at the origin point. Then, by applying a boundary nodal point collocation procedure one can obtain the following system of equations in matrix form n

_T

½A{_u} ¼ ½B{t_} þ ½E{1_ } þ ½T{b } ~ t_} þ {s_ } ¼ ½B{

~ 1_nn } ½E{

_T

~ b } þ ½T{

DSN

derived in this section for the case of a thermoviscoplastic long thick-walled cylinder problem. Those formulas have been also derived earlier by Mukherjee [17], Cordts and Kollman [18] and Arya [8]. The stress rate equation of equilibrium for the case of a thick-walled cylinder with internal and external radius a and b, respectively, which is subjected to time-dependent thermal loading is given in terms of cylindrical coordinates r and u by

› ðr s_ r Þ ¼ s_ u ›r

ð21Þ

ð18Þ

and the compatibility equation is

ð19Þ

› ðr 1_ Þ ¼ 1_r ›r u

where the coefficients of matrices [A ] and [B ] contain integrals with complex variables of the type ð ð Tijp ðPM ; QÞdSq ; Uijp ðPM ; QÞdSq ð20Þ

ð22Þ

By using Eq. (1) the total radial strain rate 1_r can be given by 1_r ¼ 1_e þ 1_n þ 1_T ¼

DSN

p and the matrices [E ] and [T ] involve Xjki and X^ pjki real ~ ~ variables integral terms while ½E; and ½T involve their complex conversions, respectively. However, the vector {1_n } is known at any time through the constitutive Eq. (2) and the stress rates of Eq. (11) while the vector {b_ T } is known through the assumed solution of the appropriate steady state diffusion equation subject to slowly varying surface temperatures. Half of the total number of components of {_u} and {t_} are prescribed through the boundary conditions while the other half are unknowns.

4.1. Time integration algorithm The initial distribution of the state variables have to be prescribed while the initial value of the non-elastic strain is set to zero. Thus, the only existing strains at time step t ¼ 0 are elastic and then, the thermal and initial stresses and displacements can be obtained from the solution of the corresponding elastic problem. By the use of Eqs. (18) and (19) the displacement and stress rates can be obtained at time step t ¼ 0 while the rates of change of the state variables can be computed from Eq. (3). Thus, the initial rates of all the relevant variables are now known and their values at a new time Dt can be obtained by integrating forward in time. The rates are then obtained at time Dt and so on, and finally the time histories of all the variables can be computed. Another important task in this approach is the choice of a suitable time integration scheme. For the purposes of the present paper, an Euler type algorithm with automatic time-step control is employed.

5. Analytical solutions For the purposes of this paper and for comparison studies reasons, the stress and strain rate analytical formulations are

809

1 ½s_ 2 nsu  þ 1_nr þ aT_ E r

ð23Þ

With similar expressions for the total circumferential strain rate 1_u : The inelastic strain rate can by found by using the correspondent Eq. (2) from Robinson’s viscoplastic model. From equations of equilibrium (21) and the expression (23) for strain rates one can obtain the equation ðr aE E _ dr þ s_ r ¼ 2 2 Tr r ð1 2 nÞ a r 2 ð1 2 n2 Þ ðr

ðr  r cðrÞdr 2 rð1nu þ n1nz Þdr a

þ

r 2 2 a2 r2

a

!

A B þ 2 2 r

ð24Þ

where A and B are constants of integration and ðr 1_n 2 1_n r u dx cðrÞ ¼ x a Eq. (24) in conjunction with the appropriate boundary conditions

s_ r ¼ 0 at r ¼ a;

s_ r ¼ 0 at r ¼ b

ð25Þ

yields to the following expressions for the radial stress rates ! E r 2 2 a2 b2 Eð1 2 2nÞ s_ r ¼ I2 þ I1 2 2 2ð1 2 n2 Þ 2ð1 2 n2 Þ b 2 a2 r 2 ! 1 r 2 2 a2 b2 2 r 2 a2 2 p _ I  2 I3 2 2 4 r b 2 a2 b2 2 a2 r 2 ! aE r 2 2 a2 ð26Þ 2 I5 2 2 I6 ð1 2 nÞr 2 b 2 a2 where ðr 1_n 2 1_n r u I1 ¼ dx x a

810

C.P. Providakis / Advances in Engineering Software 33 (2002) 805–816

Fig. 2. Geometry of the model of an experimental cylindrical thrust chamber.

6. Numerical analyses Fig. 1. Geometry of a thick-walled cylinder in plane strain subjected to inside and outside time-dependent thermal loading Ta ðtÞ and Tb ðtÞ:

ð b 1n 2 1n r u dx ¼ ðI1 Þr¼b x a ðr I3 ¼ x1_nz dx I2 ¼

a

I4 ¼ I5 ¼ I6 ¼

ðb a

ðr

x1_nz dx ¼ ðI3 Þr¼b xT_ dx

6.1. Boundary element model

a

ðb a

xT_ dx ¼ ðI5 Þr¼b

Eq. (26) in association with Eq. (25) yield the expressions for the circumferential stress rates ! E r 2 þ a2 b2 Eð1 2 2nÞ s_ u ¼ I2 2 I1 2 2 2ð1 2 n2 Þ 2ð1 2 n2 Þ b 2 a2 r 2 ! 1 r 2 þ a2 E I4 2 ð1_nu þ n1_nz Þ  2 I3 þ 2 r 1 2 n2 b 2 a2 r 2 þ b2 a2 aE þ b2 2 a2 r 2 ð1 2 n2 Þr 2 ! r 2 þ a2 2 _  I5 þ 2 I6 2 Tr b 2 a2

The details of the boundary element discretization model and thermoviscoplastic stress analysis are provided in this section. Numerical analysis has been performed and presented here for a long thick-walled cylinder of internal and external radii a and b (Fig. 1) and an experimental cylindrical thrust chamber (Fig. 2), under plane strain thermoviscoplastic deformation conditions and subjected to mechanical and thermal loading. The numerical values of the material parameters for the alloy 2.25Cr-1 Mo steel and the copper-based alloy NARloy-Z used in the subsequent analysis were taken from Refs. [8,10], respectively.

The boundary and internal element discretization models used for the viscoplastic stress and thermal analysis performed under the assumption of Robinson’s model consisted of 32 two-noded linear boundary and 35 four noded interior elements. Because the thick-walled cylinder and the cylindrical thrust chamber were symmetrical, only one-quarter of the thick-walled cylinder and one-half of the cylindrical thrust chamber cross section was modeled.

þ p_

ð27Þ

Finally, the above stress and strain rate equations together with the associated constitutive equations of Robinson’s model must be integrated to obtain the stress and strain distribution for the whole cross section of the cylinder. Time integration is performed by a scheme similar to the one employed in the above mentioned boundary element implementation. A computer program has been constructed by the same team to perform the analytical integration procedure and to compare with the results obtained by the proposed in this paper thermoviscoplastic D/BEM approach.

Fig. 3. Boundary element mesh for the thick walled cylinder.

C.P. Providakis / Advances in Engineering Software 33 (2002) 805–816

811

Fig. 4. Redistribution of circumferential stress along the radial dimension of a thick walled cylinder in plane strain at different time steps.

6.2. Thermal analysis

the thermal analysis and to generate the time-dependent cross sectional temperature rates.

The temperature rate field for the case of the thick-walled cylinder at a point with a distance r from the center of the cylinder is assumed to be given by the equation T_ ¼

1 ðT_ lnðr=aÞ þ T_ a lnðb=rÞÞ lnðb=aÞ b

ð28Þ

which actually is the solution of the steady state diffusion equation subject to slowly varying internal and external surface temperatures Ta ðtÞ and Tb ðtÞ; respectively. However, for the case of the thrust chamber the temperature rates profile was evaluated in a more advanced manner using the thermal analyzer program SINDA (see Ref. [18]) to perform

7. Results and discussion 7.1. Example 1 Consider the thick-walled cylinder of Fig. 1 with internal and external radii assumed to be equal to 0.16 and 0.25 in., respectively. The cylinder is assumed to be subjected to an internal pressure which increases linearly with time from 0 to 25.16 MPa in 10 s and it is held constant thereafter. The temperature is assumed to be constant (isothermal

Fig. 5. Error in circumferential stresses as a function of the number of internal surface element subdivisions.

812

C.P. Providakis / Advances in Engineering Software 33 (2002) 805–816

Fig. 6. Redistribution of circumferential strain along the radial dimension of a cylinder at different time steps.

conditions) throughout the cylinder at 550 8C. The boundary element mesh is shown in Fig. 3. Fig. 4 shows the hoop (circumferential) stress distribution across the wall of the cylinder as obtained by the present complex variable BEM procedure at 2.5 and 22 s and the analytical expression (22) at 2.5 s, respectively. The excellent agreement between the analytical and complex variable BEM results gives confidence in the use of the present boundary element implementation of the Robinson’s viscoplastic model. One can also see in Fig. 4 that the stress redistributes rapidly in the beginning but the rate of redistribution reduces with time. To compare the results obtained by the present D/ BEM approach with those obtained by the conventional D/ BEM of Ref. [12] the error of circumferential stresses at a point on the internal surface of the cylinder is computed by

both BEM approaches and is depicted as a function of the number of the internal surface element subdivisions required for the evaluation of internal stresses in Fig. 5. It can be shown that the divergence of the solution of the present D/BEM which employs complex variables techniques is much slower than that of the conventional D/BEM of Ref. [12]. As clearly concluded from Fig. 5 there is a significant error by both the BEM approaches up to number of 100 internal surface element subdivisions with a much less error resulted from the present D/BEM as compared to the error comes from the solution obtained by the conventional BEM. The present D/BEM gives at the point of 100 internal surface subdivisions an error of 4% while the conventional BEM gives an error of 13%. The hoop strain distribution as obtained by using the present D/BEM and the

Fig. 7. Circumferential stresses vs radial distance of a cylinder at different time steps.

C.P. Providakis / Advances in Engineering Software 33 (2002) 805–816

Fig. 8. Circumferential strains vs radial distance of a cylinder at different time steps.

analytical procedure described above procedure has been plotted in Fig. 6 for two different values of time. Fig. 6 shows that the maximum and minimum values of strain occur at the internal and external surface of the cylinder, respectively. 7.2. Example 2 Consider the same problem as the above cylinder but the temperature distribution is held constant at 430 8C (isothermal condition). Figs. 7 and 8 show the hoop stress and strain distributions along the radial axis at different time steps. In order to further verify the BEM results and to ensure the correct implementation of Robinson’s model in the present code, results as obtained by the finite element method of Arya [8] are compared to the results yielded by the present D/BEM procedure which employs complex variables numerical techniques. Both FEM and present D/BEM results are plotted in Figs. 7 and 8 and compare well. 7.3. Example 3 Consider an experimental thrust chamber with a cross section shown in Fig. 2 which was also investigated in Arya and Arnold [10] by a finite element approach. The model of the section of the thrust chamber, analyzed here, is shown by the shaded area in this figure and is discretized by a mesh shown in Fig. 9 to perform a complicated thermomechanical analysis under mechanical pressure and thermal loading histories shown in Fig. 10. Fig. 11 shows the circumferential displacements at C of the model as obtained by the present D/BEM and the finite element solution of Ref. [10] as a function of the number of the loading cycles. The displacement values depicted are the values at the

Fig. 9. Mesh of the experimental thrust chamber model.

813

814

C.P. Providakis / Advances in Engineering Software 33 (2002) 805–816

Fig. 10. Pressure and thermal loading histories.

completion of the corresponding loading cycle. The values of the displacement at point C shown in Fig. 11 also give the change in thickness between points C and D since point D is constrained. A close examination of this figure reveals the experimentally verified in Ref. [19] and so-called ‘doghouse effect’ that greater thinning of the channel wall across section CD is implied as the number of cycles is considered to be increased. Fig. 12 shows the radial displacement at the points A and B with respect of the number of cycles. A doghouse effect is also implied from the inspection of Fig. 12 since the difference in the displacement at points A and B is actually the change in wall thickness (thinning) across the section AB of the model. The very good correlation existed

between the results obtained by the present D/BEM and the FEM of Arya and Arnold [10] also proves the accuracy and the efficiency of the use of the proposed methodology to the implementation of Robinson’s model. One can also observe that FEM requires almost the same d.o.f with the present D/ BEM for the same level of accuracy.

8. Conclusions This paper presents, for the first time, an effective implementation of Robinson’s viscoplastic model to the direct formulation of a D/BEM for the analysis of creeping

Fig. 11. Circumferential displacements at point C.

C.P. Providakis / Advances in Engineering Software 33 (2002) 805–816

815

Fig. 12. Radial displacements at points A and B.

metallic structures subjected to mechanical and thermal loading histories. A new numerical methodology based on complex variable techniques is introduced to the evaluation of internal stresses with integral equations only weakly singular. To validate the results an analytical solution of the problem is also presented. On the basis of the excellent agreement between the analytical and the present boundary element results the proposed approach is proved to be a powerful tool for the solution of complicated inelastic problems of practical importance as in the cases of heat-exchanger metallic tubes and pressure vessels. It can be also concluded that the present boundary element methodology provides results with the same accuracy as the conventional method and has the potential of reduced computation time. The proposed D/BEM algorithm is quite general and can solve in a easy manner the problem of the computation of internal stresses without any analytic differentiation of the displacement integral equations. Of course, the present results cannot be considered as conclusive since the proposed approach needs to be further compared, about its possibilities as a competitive computational method, for more complex engineering problems, with other numerical solutions especially with FEM approaches which is an immediate follow up step for this kind of BEM analysis.

Acknowledgements The author is grateful to Professors D.E. Beskos for encouragement and helpful discussions during the course of this work.

References [1] Walker KP. Research and development program for nonlinear structural modelling with advanced time-temperature dependent constitutive relationships. PWA-5700-50, Pratt and Whitney Aircraft, NASA Contract NAS3-22055, NASA CR165533; 1981. [2] Lindholm US, Chan KS, Bodner SR. Constitutive modelling for isotropic materials. NASA CR-174718; 1984. [3] Lemaitre J, Chaboche JL. Mechanics of solid materials. New York: Cambridge University Press; 1990. [4] Freed AD, Walker KP. Viscoplasticity with creep and plasticity bounds. Int J Plasticity 1993;9:213–42. [5] Saleeb AF, Arnold SM, Castelli MG, Wilt TE, Graf W. A general hereditary multimechanism-based deformation model with application to the viscoelastoplastic response of titanium alloys. Int J Plasticity 2001;17:1305–50. [6] Robinson DN, Swindeman RW. Unified creep-plasticity constitutive equations for 2.25Cr– 1Mo steel at elevated temperature. ORNL TM8444, Oak Ridge National Laboratory TN; 1982. [7] Robinson DN, Bartolotta PA. Viscoplastic constitutive relationships with dependence on thermomechanical history. Cleveland, OH: Lewis Research Center; 1985. NASA CR-174836. [8] Arya VK. Analytical and Finite element solutions of some problems using viscoplastic model. Comput Struct 1989;33(4):957 –67. [9] Arya VK. Application of finite element based solution technologies for viscoplastic structural analysis. NASA Report CR185196; 1990. [10] Arya VK, Arnold SM. Viscoplastic analysis of an experimental cylindrical thrust chamber liner. AIAA J 1992;30(3):781– 9. [11] Arya VK, Kaufman A. Finite element implementation of Robinson’s unified viscoplastic model and its application to some uniaxial and multiaxial problems. Engng Comput 1989;6:237–47. [12] Providakis CP, Kourtakis S. BEM implementation of Robinson’s model to the thermoviscoplastic response of metallic structural components. In press. [13] Gao XW, Davies TG. An efficient boundary element algorithm for 2D and 3D elastoplastic problems. Int J Solids Struct 2000;37: 4987– 5008.

816

C.P. Providakis / Advances in Engineering Software 33 (2002) 805–816

[14] Lyness JN, Moler CB. Numerical differentiation of analytic functions. SIAM J Numer Anal 1967;4:202–10. [15] Mukherjee S. Corrected boundary integral equations in planar thermoelastoplasticity. Int J Solids Struct 1977;13:331 –5. [16] Providakis CP. A general and advanced boundary element transient analysis of elastoplastic plates. Engng Anal Bound Elem 1996;17: 133–43. [17] Mukherjee S, Thermoviscoplastic response of cylindrical structures

using a state variable theory, vol. 2. Cambridge, England: ICM 3; 1979. [18] Smith JP. Systems improved numerical differencing analyzer (SINDA): user’s manual. Redondo Beach, CA: TRW Systems Group; 1971. TRW-14690-H001R0-00. [19] Quentmeyer RJ. Experimental fatigue life investigation of cylindrical thrust chambers. NASA TM X-73665; 1977.