3D modeling of crack growth in electro-magneto-thermo-elastic coupled viscoplastic multiphase composites

Available online at www.sciencedirect.com

Applied Mathematical Modelling 33 (2009) 1014–1041 www.elsevier.com/locate/apm

3D modeling of crack growth in electro-magneto-thermo-elastic coupled viscoplastic multiphase composites Bojing Zhu a,b, Taiyan Qin a,* a

b

College of Science, China Agriculture University, Beijing 100083, PR China College of Earth Sciences, Graduate University of Chinese Academy of Sciences, Beijing 100049, PR China Received 10 March 2007; received in revised form 13 December 2007; accepted 19 December 2007 Available online 1 January 2008

Abstract This work presents a time-domain hypersingular integral equation (TD-HIE) method for modeling 3D crack growth in electro-magneto-thermo-elastic coupled viscoplastic multiphase composites (EMTE-CVP-MCs) under extended incremental loads rate through intricate theoretical analysis and numerical simulations. Using Green’s functions, the extended general incremental displacement rate solutions are obtained by time-domain boundary element method. Three-dimensional arbitrary crack growth problem in EMTE-CVP-MCs is reduced to solving a set of TD-HIEs coupled with boundary integral equations, in which the unknown functions are the extended incremental displacement discontinuities gradient. Then, the behavior of the extended incremental displacement discontinuities gradient around the crack front terminating at the interface is analyzed by the time-domain main-part analysis method of TD-HIE. Also, analytical solutions of the extended singular incremental stresses gradient and extended incremental integral near the crack fronts in EMTE-CVP-MCs are provided. In addition, a numerical method of the TD-HIE for a 3D crack subjected to extended incremental loads rate is put forward with the extended incremental displacement discontinuities gradient approximated by the product of time-domain basic density functions and polynomials. Finally, examples are presented to demonstrate the application of the proposed method. Ó 2007 Elsevier Inc. All rights reserved. Keywords: Electro-magneto-thermo-elastic coupled viscoplastic multiphase composites; Three-dimensional crack; Time-domain hypersingular integral equation; Extended incremental integral; Boundary element method

1. Introduction The development of piezoelectric and piezomagnetic composites has its roots in the early work of [1–5], who proposed that the combination of piezoelectric/piezomagnetic phases may exhibit a new material property-the electromagnetic coupling effect. Since then, the combination coupling effect of two or more phases in *

Corresponding author. Tel.: +86 10 62736992; fax: +86 10 62736777. E-mail addresses: [email protected] (B.J. Zhu), [email protected] (T.Y. Qin).

0307-904X/$ - see front matter Ó 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2007.12.024

B.J. Zhu, T.Y. Qin / Applied Mathematical Modelling 33 (2009) 1014–1041

Nomenclature RiJ F_ J r_ ij D_ i B_ i #_ i f_ i f_ e f_ m f_ # r_ nij e_ piJ e_ ciJ e_ TiJ e_ M iJ aiJ rY k_ b_ w ^ a T_ E CiJKl PiJ cijkl elij dlij eil gil lil iij Bil gil kij Z_ Kl U_ K u_ i /_ u_ !_ S p X U_ IJ T_ IJ e_ IJk T_ J

the second-order extended incremental stress displacement rate matrix tensor the extended dummy incremental body load column vector the second-order incremental stress tensor the incremental electric displacement vector the incremental magnetic flux vector the incremental heat flow vector the incremental body force vector the incremental free charge density vector the incremental magnetic flux density vector the incremental heat flow density vector the second-order nonlinear incremental stress tensor the second-order plastic strain rate tensor the second-order creep strain rate tensor the second-order thermal strain rate tensor the second-order thermo-material strain rate tensor the second-order hardening tensor the yield surface radius the positive linear constant the positive linear coefficient the creep potential the thermal expansion coefficient the thermal increment the electro-magneto-thermo-elastic constant matrix the fourth-order extended elastic stiffness tensor the second-order extended thermal stress constants tensor the fourth-order elastic stiffness tensor the third-order piezoelectric constants tensor the third-order piezomagnetic magnetoelectric coupling tensor the third-order dielectric permittivities tensor the second-order magnetic permeabilities tensor the second-order magnetic constants tensor the second-order thermal stress constants tensor the second-order pyroelectric coefficients tensor the second-order pyromagnetic coefficients tensor the second-order heat conduction coefficients tensor the second-order extended incremental strain field matrix the extend incremental displacement rate vector the incremental displacement rate vector the incremental electric potential rate vector the incremental magnetic potential rate vector the incremental thermal potential rate vector the crack surface interior point in Domain the domain occupied by the EMTE-CVP-MCs the extended displacement fundamental incremental solutions the extended traction fundamental incremental solutions the extended strain fundamental incremental solutions the extended incremental traction

1015

1016

p_ i q_ b_ ^_ # S+ S e_ J U ~u_ j ~_ j / ~_ j u ~_ j !

B.J. Zhu, T.Y. Qin / Applied Mathematical Modelling 33 (2009) 1014–1041

the the the the the the the the the the the

incremental mechanical loads rate incremental electrical loads rate incremental magnetic loads rate incremental thermal loads rate upper side of the crack surface S lower side of the crack surface S extended incremental displacement discontinuities gradient incremental displacement discontinuities gradient incremental electric potential discontinuities gradient incremental magnetic potential discontinuities gradient incremental thermal potential discontinuities gradient

electro-magneto-thermo-elastic coupled viscoplastic multiphase composites (EMTE-CVP-MCs) have been measured by many researches [6]. But most of these theoretical achievements are concerned with material properties and characters, little numerical modeling and result on 3D crack expand and spread problem has been investigated. Nowadays, EMTE-CVP-MCs are being used more and more widely in space planes (sensors), supersonic airplanes (gas turbines), rockets, missiles, nuclear fusion reactor, and submarines (smart skin systems) due to their capabilities to respond in an advantageous manner to changing environment. EMTE-CVP-MCs typically exhibit pronounced nonlinear response under sufficiently higher EMTE incremental loading conditions. These extended incremental loads are causing higher extended incremental stresses than under corresponding static extended incremental loads and may induce crack initiation, crack growth and finally lead to fracture or failure of structures. Due to the importance of the reliability of these structures, viscoplastic fracture and failure analysis of EMTE-CVP-MCs have received considerable interest in the recent decade. For the 2D nonlinear crack problem, using finite difference method, Aboudi and Achenbach [7,8] considered Modes I and III crack in elastic–viscoplastic material, and investigated the dynamic effects due to high crack-tip speeds and compared the dynamic results with those obtained on the basis of quasi-static problem formulation. Riedel [9] evaluated the crack growth in elastic–viscoplastic solids, and obtained that the governing load parameter was the stress intensity factor (SIF) K1 if the bulk of the specimen was predominantly elastic and it was the J-integral in a fully-plastic situation when large creep strains were still confined to a small zone. The time-dependent crack problem in inelastic deformation plates was analyzed by Morjaraia [10]. By applying Green’s functions method, Achenbach [11] considered crack-tip fields in viscoplastic material, and then obtained that the leading term of the near-tip particle velocity was of order r1/4, and the higher-order terms were of the forms rlog r and r, the rlog r term was absent for the Mode III case. After analyzed creep crack initialization and growth problem by a very simple local fracture law, Bui et al. [12] discussed creep zone size, stress redistribution, and influence of local fracture criterion. Based on the finite element method, Krishna et al. [13] investigated steady-state dynamic crack growth problem in viscoplastic material, and analyzed the effect of rate sensitivity and crack velocity on the fracture toughness and the temperature distribution. In the analysis of anti-plane steady-state fast-crack propagation problem in an elastic–viscoplastic material, Wu and Lai [14] found that the stress field near the tip of a semi-infinite crack could be characterized by the sum of an elastic SIF Ka and plastic SIF Kp. Gao [15] presented time-dependent crack tip stress field in elastic–viscoplastic material and obtained that the power law singular behavior prevailed for high viscosity while the logarithm singular behavior applied to low viscosity. Using degenerate transform method, the nonlinear singular part of 1/r3 was developed by Feng [16]. Applying boundary element method, Baveendra [17] considered the SIF of steady planar thermoplasticity crack problem. Based on boundary integral equation method, Sung and Liou [18] observed transient crack growth problem in a material under small-scale yielding conditions, and derived the near-tip time-dependent SIF and discussed the relationship between numerical results and crack tip speeds. Zhu and Cescotto [19] presented an energy-based anisotropic elasto-visco-plastic damage model at finite strain, and demonstrated that the anisotropic damage model was suitable for large-scale computation. In analysis of fast crack propagation under dynamic Mode I load conditions, a finite element model was developed by

B.J. Zhu, T.Y. Qin / Applied Mathematical Modelling 33 (2009) 1014–1041

1017

Leppin and Wriggers [20], and the numerical results demonstrated that these physical effects were of major importance when load speed was increased to a high level and could not be neglected. Dynamic central crack subject to impact tensile loading inisotropically hardening elastic–viscoplastic solid was considered by Siegmund and Needleman [21], and they found that SIFs increased dramatically at a certain value of the crack speed that depended on the cohesive surface strength, the material flow strength, the characterization of strain rate hardening and the impact velocity. By using hybrid displacement discontinuity method, a crack and multiple crack fatigue growth problems have been investigated by Yan [22,23]. To the 3D nonlinear crack problem, using boundary element method, 3D elastoplastic permeable surface crack was analyzed by Wang [24]. Applying dual boundary element method, Cisilino [25] studied the 3D elastoplastic fracture mechanics. Based on elastic–plastic finite element model method, crack growth problem was considered by Zhang [26]. By boundary element method combined with the line-spring model method, Liang [27] considered the elastoplastic and creep surface crack problems. Green’s functions for incremental nonlinear elasticity was obtained by Bigoni and Capuani [28]. Using initial stress boundary element method employing an iteration procedure in each load increment, Sladek [29] considered the strain-softening damage problems. Oleaga [30] introduced a 3D dynamic crack propagation law by Hamilton principle. The application of hypersingular integral equation method to linear elastic fracture mechanics is now well established and used in electromagnetic composites. The accuracy and efficiency of the hypersingular integral equation method suggested that its extension to electromagnetic thermoelastoplastic viscoplastic problems may also be effective. Nonlinear boundary element method first developed by Swedlow and Cruse [31], and they gave the extended form of Somigliana’s identity. Using boundary element method, Riccardella [32] first obtained the numerical result of elastoplasticity planar problems. Boundary integral equation for time-dependent inelastic deformation problems were considered by Kumar, Mukherjee and Mendelson [33–35]. The integral form of the gradient rate for 3D thermoelastoplastic problem was investigated by Bui [36], and he corrected some erroneous expressions reported in the literature. Mukherjee and Kumar [37] obtained the time-dependent inelastic deformation problems. Telles and Brebbia [38] investigated the inelastic boundary integral equations for 2D and 3D problems. The inelastic boundary integral equation in initial strain and stress was obtained by Brebbia and Telles [39,40]. A 3D elastoplasticity and viscoplasticity problem was considered by Banerjee and Davies [41]. Dargush [42] developed 3D thermoplasticity problem in theory. 3D hypersingular boundary integral equation in elastoplasticity solids was analyzed by Huber et al. [43], and numerical method of strongly singular domain integrals was carried out in the paper. Sladek [44] considered displacement gradient within boundary element formulation for small strain plasticity problems. By applying symmetric Galerkin boundary element method, dynamic elastic–plastic problem was investigated by Frangi and Maier [45]. Using time-independent Green’s functions and boundary element method, Liu [46] analyzed nonlinear surface crack problem. Gao and Davies [47] observed removing strong singularities in domain integrals of elastoplasticity problem. Based on 3D boundary element method, Hatzigeorgion and Beskos [48] considered the transient dynamic elastoplastic problem. Brun et al. [49,50] developed the homogeneous incremental nonlinear elasticity strain and obtained some numerical results. Using time-domain boundary element method, Soares et al. [51] investigated nonlinear elastodynamic problem. Recently work by Bertoldi et al. [52] has handled elastoplastic problems by boundary element method without domain integrals. However, relatively little work has been done on problem of 3D crack growth in EMTE-CVP-MCs, because of the current limitations both practical (computing time) and theoretical (accurate 3D formulations of dislocation shielding and image force). There exists a need for a general purpose, accurate theoretical model that can provide efficient numerical and analysis of the problem. The time-domain hypersingular integral equation (TD-HIE) method based on the incremental displacement discontinuities fundamental solutions take the crack’s upper and lower surfaces as one boundary, so it needs no subdivision of the domain. Combined with the time-domain finite-part integral method [53], the method is an efficient method to analyze crack problems. In this paper, the general solutions of extended incremental displacement rate (the creep viscoplastic incremental displacement rate, the electrical incremental potential rate, the magnetic incremental potential rate and the thermal incremental potential rate) are obtained by time-domain boundary element method. Then, applying time-domain main-part method and extended boundary conditions, a 3D crack growth problem in EMTE-CVP-MCs subjected to the extended incremental loads rate (the incremental mechanical load rate,

1018

B.J. Zhu, T.Y. Qin / Applied Mathematical Modelling 33 (2009) 1014–1041

the incremental electrical load rate, the incremental magnetic load rate and the incremental thermal load rate) is reduced to solve a set of TD-HIEs. The unknown functions are the extended incremental displacement discontinuities gradient (the incremental creep viscoplastic displacement discontinuities gradient, the incremental electrical potential gradient, the incremental magnetic potential gradient and the incremental thermal potential gradient) of the crack surface. Also, the singularity of the extended incremental displacement discontinuities gradient is analyzed. The analytical solution of the extended incremental singular stresses gradient (the incremental singular stress gradient, the incremental singular electric displacement gradient, the incremental singular magnetic displacement gradient and the incremental singular thermal stress gradient) and extended incremental J_ gi integral near the crack front in EMTE-CVP-MCs are given. In addition, the numerical method of the TD-HIEs for a 3D crack subjected to extended incremental loads rate is proposed and some numerical solutions are calculated. For the special case, the results are compared with those obtained by Providakis, Li, and Morjaria and Mukherjee [54,55,10]. Finally, the results show that the present method yields smooth variations of extended crack incremental opening displacements (CIOD) along the crack front accurately. 2. Basic equations and mathematical model The nonlinear governing equations and constitutive relationships for EMTE-CVP-MCs can be expressed by incremental tensors forms RiJ ;i þ F_ J ¼ 0:

ð1Þ

In this paper, summation from 1 to 3 over repeated lowercase (1 to 6 uppercase) subscripts is assumed, and a subscript comma denotes the partial differentiation with respect to the coordinates (i.e., x1, x2, x3, or x ,y, z), hi is the McAuley symbol. In addition, the extended incremental stress displacement rate and the extended dummy incremental body loads, RiJ, and F_ J , are defined, respectively, by 8 r_ ij J ¼ j ¼ 1; 2; 3; > > > < D_ J ¼ 4; i ð2Þ RiJ ¼ _ > J ¼ 5; B i > > :_ #i J ¼ 6; 8 _ > f j  r_ nij J ¼ j ¼ 1; 2; 3; > > > < _ J ¼ 4; f e ð3Þ F_ J ¼ > _ > J ¼ 5; > f m > : _ f # J ¼ 6: Using the extend Hart’s elastoplastic creep models, a small displacement and a small strain formulation are employed, the nonlinear incremental stress, r_ nij , can be defined as r_ nij ¼ EiJKl e_ niJ ¼ EiJKl ð_epiJ þ e_ ciJ þ e_ TiJ þ e_ M iJ Þ:

ð4Þ

The mathematical structure of many of the state variable models of extended elastoplastic creep deformation can be summarized by the following equations: _ ðriJ ; aiJ ; rY Þ=orij ; e_ piJ ¼ koF _ e_ c ¼ bow=or iJ ;

iJ e_ TiJ e_ M iJ

¼^ aT_ diJ ; ¼ rKl T_ oEiJKl =oT :

ð5Þ ð6Þ ð7Þ ð8Þ

The elastoplastic creep incremental constitutive equations are written as RiJ ¼ EiJKl Z_ Kl ;

ð9Þ

where the electro-magneto-thermo-elastic (EMTE) constant and the extended incremental strain, EiJKl and Z_ Kl , take the form,

B.J. Zhu, T.Y. Qin / Applied Mathematical Modelling 33 (2009) 1014–1041

 EiJKl ¼

C iJKl

PiJ

0

0

0 T

kiJ

 Z_ Kl ¼ U_ K;l

!_ ;l

1019

 ð10Þ

;

ð11Þ

;

where the extended elastic stiffness, the extended thermal stress constants and the extend incremental displacement rate, CiJKl PiJ and U_ K , are defined, respectively, by 8 cijkl ; J ; K ¼ 1; 2; 3; > > > > > J ¼ 1; 2; 3; K ¼ 4; elij ; > > > > > e ; J ¼ 4; K ¼ 1; 2; 3; ikl > > > d ikl ; J ¼ 5; K ¼ 1; 2; 3; > > > > > > > gil ; J ¼ 4; K ¼ 5 or J ¼ 5; K ¼ 4; > > > eil ; J ; K ¼ 4; > > : lil ; J ; K ¼ 5; 8 > < iij ; J ¼ 1; 2; 3; K ¼ 6; PiJ ¼ 1il ; J ¼ 4; K ¼ 6; ð13Þ > : gil ; J ¼ 5; K ¼ 6; 8 u_ i ; K ¼ 1; 2; 3; > > > < /; _ K ¼ 4; U_ K ¼ ð14Þ > _ u; K ¼ 5; > > :_ !; K ¼ 6: 3. Time-domain boundary integral equations for a crack in EMTE-CVP-MCs 3.1. Boundary conditions of a crack surface The inelastic mechanical boundary conditions of crack are always defined by stress free crack surface. Several EMTE boundary conditions were proposed in the literature. Among these EMTE boundary conditions, two different conditions are applied widely. Those are permeable and impermeable conditions. For the first one, the normal extended incremental displacement rate and extended incremental potential rate should be continuous across the crack surface, _ D_ þ 3 ¼ D3 ;

/_ þ ¼ /_  ;

_ B_ þ 3 ¼ B3 ;

u_ þ ¼ u_  ;

_ #_ þ 3 ¼ #3 ;

!_ þ ¼ !_  ;

ð15Þ

where superscripts + and  denote the upper and lower crack surface, respectively. Proposed impermeable conditions on the crack faces, _ D_ þ 3 ¼ D3 ¼ 0;

/_ þ ¼ /_  ¼ 0;

_ B_ þ 3 ¼ B3 ¼ 0;

u_ þ ¼ u_  ¼ 0;

_ #_ þ 3 ¼ #3 ¼ 0;

!_ þ ¼ !_  ¼ 0:

ð16Þ

This paper presents an analysis for the crack problems based on boundary condition (16). 3.2. Fundamental solutions for a crack in 3D EMTE-CVP-MCs Consider a 3D crack S in finite EMTE-CVP-MCs. A fixed rectangular Cartesian system xi (i = 1, 2, 3) is used. The crack is assumed to be in the x1x2 plane, and normal to the x3-axis. Using the EMTE form of the Somigliana’s identity, the extended incremental displacement rate at interior point p are expressed as Z Z Z e_ IJk r_ n dv; ðT_ IJ U_ J þ U_ IJ T_ J Þds þ ð17Þ U_ IJ f_ J dv þ U_ I ¼  S þ þS 

Jk

X

X

1020

B.J. Zhu, T.Y. Qin / Applied Mathematical Modelling 33 (2009) 1014–1041

where the relationships between the extended fundamental incremental solutions and Green’s functions take form as T_ IJ ¼ EkJMn U_ IM;n nk ; e_ IJk ¼ U_ IJ ;k :

ð18Þ ð19Þ

The extended incremental traction (the elastic incremental traction, the normal incremental charge flux densities, the normal incremental magnetic current densities, and the normal incremental heat flow densities), T_ J , on the boundary and can be defined as 8 p_ j ¼ r_ jl nl  r_ niJ ;i ; J ¼ j ¼ 1; 2; 3; > > > > < q_ ¼ D_ l nl ; J ¼ 4; ð20Þ T_ J ¼ _ _ > b ¼ Bl nl ; J ¼ 5; > > > : ^_ J ¼ 6; # ¼ #_ ;J ni ; where the incremental mechanical loads rate, the incremental electrical loads rate, the incremental magnetic ^_ are on the crack surface due to internal _ and #, _ b, loads rate and the incremental thermal loads rate, p_ i , q, or external incremental loads rate, and can be obtained from the solution for the incremental loads rate of the un-cracked solids. Noting that the extended incremental displacement discontinuities gradient is written as 8 > ~ _ u_ j ¼ u_ þ J ¼ j ¼ 1; 2; 3; > j u j ; > > > > ~ _ _ u ¼ u  u J ¼ 5; > j j j ; > > > > :! ~_ j ¼ !_ þ  !_  ; J ¼ 6: j j ^ ; tÞ ¼ T_ þ Using the relationships T_ IJ ð^ p; ^ qþ ; tÞ ¼ T_ IJ ð^ p; q p; ^q; tÞ, and U_ IJ ðp; qþ ; tÞ ¼ U_ IJ ðp; q ; tÞ, where ^p, and IJ ð^ ^ q are the source point and the field point, respectively. The extended incremental displacements rate (17) can be rewritten as Z Z Z e_ J ds þ _ IJ f_ J dv þ U e_ IJk r_ nJk dv ð22Þ T_ þ U U_ I ¼  IJ Sþ

X

X

using solution (21) and constitutive equation (9), the corresponding second-order incremental stress displacement rate matrix tensor, RiJ, is shown as Z Z Z _ þ e _ _ _ ð23Þ DKiJ f K dv  W_ KiJl e_ nJK dv; S KiJ U K ds þ RiJ ¼  Sþ

X

X

where the integral kernel functions are as follows: S_ KiJ ð^ p; ^ q; tÞ ¼ EiJMn T_ MK;n ð^ p; ^ q; tÞ; _ _DKiJ ð^ p; ^ q; tÞ ¼ EiJMn U MK;n ð^ p; ^ q; tÞ;

ð24Þ

p; ^ q; tÞ ¼ EiJMn ðU_ KJ ;Mn ð^ p; ^ q; tÞ þ U_ KM;Jn ð^p; ^q; tÞÞ: W_ KiJn ð^

ð26Þ

ð25Þ

3.3. Time-domain boundary integral equations For a 3D crack in finite EMTE-CVP-MCs, S± is the crack surface. Using the boundary conditions, the time-domain boundary integral equations and time-domain hypersingular integral equations can be obtained. Let the source point p be taken to the boundary C and represented by P, the boundary integral equation can be derived from Eq. (22), as follows: Z Z Z e_ J ds þ _ IJ f_ J dv þ U e_ IJk r_ nJk dv; ð27Þ C IJ U_ I ¼  T_ þ U IJ Sþ

X

X

B.J. Zhu, T.Y. Qin / Applied Mathematical Modelling 33 (2009) 1014–1041

1021

where CIJ is the constant related to the boundary point P. Using the elastoplastic creep condition, the electric condition, the magnetic condition and the thermal boundary condition of the crack surface, the time-domain hypersingular integral equations can be obtained as Z Z _ _ K dv  _þ U e _ S f ¼ ds þ ð28Þ D W_ KiJl e_ nJK dv; Rþ K KiJ þ KiJ iJ S

X

X

where means that the integral must be interpreted as a time-domain finite-part integral. The first integral in (28) has the order r3, and is a hypersingular one, solving Eq. (28), all the unknown functions can be obtained. 4. A crack in electro-magneto-thermo-elastic coupled viscoplastic multiphase composites For transversely isotropic EMTE-CVP-MCs, the fourth-order EMTE constant tensor, cijkl, and the thirdorder EMTE constant tensor, elij, can be written as follows: cijkl ¼ ðc13  c12 Þðdij d3k d3l þ d3i d3j dkl Þ þ ðc44  c66 Þðdjk d3i d3l þ dik d3j d3l þ dil d3j d3k þ djl d3i d3k Þ þ ðc11 þ c33  2c13  4c44 Þd3i d3j d3k d3l þ c12 dij dkl þ c66 ðdik djl þ dil djk Þ;

ð29Þ

elij ¼ e31 dij d3l þ e15 ðdil d3j þ djl d3i Þ þ ðe33  e31  2e15 Þd3i d3j d3l ;

ð30Þ

eil ¼ e11 dil þ ðe33  e11 Þd3i d3l :

The format of the constant tensor dlij is identity to the constant elij, the format of the constant lil and gil is identity to the constant eil, here c66 = (c11  c12)/2. 4.1. General solutions For transversely isotropic EMTE-VP-MCs, using the method given in Ref. [28,56–60], general solutions of extended incremental displacement rate can be written as an explicit expression. _ incremental monopole J_ , and incremental 4.1.1. Combination of point incremental force P_ z , incremental charge Q, _ thermal force ! _ J_ , and !_ at point ^pðn1 ; n2 ; n3 ; tÞ, the extended incremental displacements rate at For a combination of P_ z , Q, point ^ qðx1 ; x2 ; x3 ; tÞ can be expressed as !, 5 5 X X b 1 ðxa  na ÞdaJ þ j1 vi R1 Aim R jin dðnþ1ÞJ  R2 ðx3  n3 Þji5 v2 d6J 4p; ð31Þ G_ mJ ¼ i

i

i1

i¼1

i

i

n¼2

where a ¼ 1; 2; m ¼ 3–6; b i ¼ Ri þ vi z; R 2

2

2 0:5

Ri ¼ ððx1  n1 Þ þ ðx2  n2 Þ þ ðvi zÞ Þ ; T

Ai ¼ ½Ai1 Ai2 Ai3 Ai4 0 ¼ W 1 ; 0:5 v1 ¼ c0:5 66 c44 ; 0:5 v6 ¼ k0:5 11 k33

and vt(t = 2–5) are the four roots of the algebra equations: n1v8  n2v6 + n3v4  n4v2 + n5 = 0 other parameters are listed in Appendix A. 4.1.2. Point incremental force in x2 direction For a unit point incremental force in x2 direction at point ^pðn1 ; n2 ; n3 ; tÞ, the extended incremental displacements rate at point ^ qðx1 ; x2 ; x3 ; tÞ can be expressed as

1022

B.J. Zhu, T.Y. Qin / Applied Mathematical Modelling 33 (2009) 1014–1041

G_ 2J ¼

5 X

Di ðx2  n2 Þ

b 2 daJ R1 i R i ðna

 xa Þ 

j1 i1 vi

b 1 R1 i Ri

5 X

i¼1

n¼2

 ððx1 

b 2 n1 ÞR1 0 R 0 ððx1

! jin dðnþ1ÞJ 

ji5 vi R3 i d6J

! þ

b 1 d2J R i

!,

b 1 d2J Þ  n1 Þd2J  ðx2  n2 Þd1J Þ  R 0

4pc44

ð32Þ

in which Di ¼ V^ 1 ½ 0 0

0

0

T

1 ;

^vil ¼ ðjil vi dkl þ j55 di5 d4l Þj1 i1 þ vi d5l : 4.2. Time-domain hypersingular integral equations Consider a 3D crack in EMTE-CVP-MCs as shown in Fig. 1. A fixed global rectangular Cartesian system xi (i = 1, 2, 3) is chosen. Assume that a stochastic crack S(S+ [ S-) is subjected to remote the incremental mechanical loads p_ i , the incremental electrical loads p_ 4 ðq_ 0 Þ, the incremental magnetic loads p_ 5 ðb_ 0 Þ and the incremental thermal loads p_ 6 ð#_ 0 Þ, respectively. The local rectangular Cartesian system ni are chosen, the stochastic flaws are assumed to be in the n1 ni2 plane and normal to the n3-axis, the angle between the fixed global axis xi and local Cartesian ni is defined as wi(xi, ni, t). Using the method given in [53], Eq. (28), can be reduced to   5 5 P P 3 ^i2 ~u_ b;a ð^qÞ þ 3r4 r;a k33 v2i q ^i1 ~u_ 6;a ð^qÞdsð^qÞ r c44 v0 ð3r;a r;b  2dab Þ=4p þ ðdab  3r;a r;b Þ qi2 q Sþ i¼1 i¼1 ð33Þ 5 ^ R P p 3 2 c T M w r ð1  3r;a Þð_ejK þ e_ jK þ e_ jK þ e_ jK ÞdV ð^qÞ ¼ p_ a ð^pÞ þ X a  i¼1 5  5 P 5 5 P P P 2 3 4 3 1 _ _ _ ^ ^ ~ ~ ~ r r A q ð^ q Þ þ r q ð^ q Þ þ 3r k r v j j q ð^ q Þ dsð^qÞ q q u u u ;a i4 i2 i1 a;m 3a ;a im it n;m i 5i i1 im 6;m Sþ i¼1 n¼3 i¼1 i¼1 ð34Þ 5 R P ^ ^ m Þð_epjK þ e_ cjK þ e_ TjK þ e_ M _ r3 ðwð1  r2;a Þ þ w ÞdV ð^ q Þ ¼  p ð^ p Þ þ X m jK m  i¼1   5 5 P P þ 2 2 2 6_ _ ^ ub;6 ð^ r ðdab  3r;a r;b Þ~ qÞ Ai4 k3b qi2 3r k3a r;a Ai vi qi ~u6;6 ð^qÞ dsð^qÞ Sþ i¼1 i¼1 ð35Þ 5 R P p 3 ^ 2 c T M ^ m Þð_ejK þ e_ jK þ e_ jK þ e_ jK ÞdV ð^qÞ ¼ p_ 6 ð^pÞ r ðwð1  r;a Þ þ w þ X i¼1

m

. . . . σ 33∞ , D33∞ , B33∞ ,ϑ33∞ ⊗

⊗. σ31∞

⊗

⊗ x3

ξ3

⊗

ω3

σ32∞

ξ2

S±

ω2 O′

ω1

ξ1 x2

O

. σ32∞

•

x1

•

•

•

. . . . σ 33∞ , D33∞ , B33∞ ,ϑ33∞

. σ31∞ •

Fig. 1. A crack in electro-magneto-thermo-elastic coupled viscoplastic multiple composites.

B.J. Zhu, T.Y. Qin / Applied Mathematical Modelling 33 (2009) 1014–1041

1023

where p_ i , p_ 4 ðq_ 0 Þ, p_ 5 ðb_ 0 Þ, and p_ 6 ð#_ 0 Þ can be obtained from the solution for the incremental loads rate of the uncracked solids. The material constants are given in Appendix A. It shows that the time-domain hypersingular integral equations have the similar structures as that studied in [61–63]. 4.3. Extended singular incremental stresses gradient and extended incremental integral J_ gi In order to investigate the singularity of the crack front, consider a local coordinate system defined as x1x2x3, as shown in Fig. 2, the x1-axis is the tangent line of the crack front at point q0, x2-axis is the internal normal line in the crack plane, and x3-axis is the normal of the crack. Then the extended incremental displacement discontinuities gradient of the crack surface near a crack front point ^q0 can be assumed as e_ i;J ¼ gk ð^ U q0 ; tÞnk2k ;

0 < Reðkk Þ < 1;

ð36Þ

q0 ; tÞ are non-zero constants related to point ^q0 , kk are the singular indexes at the crack front. Conwhere gk ð^ sider a small semi-circle domain Se on the crack surface including point ^q0 , using the time-domain main-part analytical method [53], the following relationships can be derived: Se

u_ i;1 dn1 n2 ffi 2pk1 g1 xk21 1 cotðk1 pÞ; r3 ~

ð37Þ

Se

u_ i;1 dn1 n2 ffi 4pk1 g1 xk21 1 cotðk1 pÞ; 3r5 ðx2  n2 Þ2 ~

ð38Þ

Se

2_ 3r5 ðx1  n1 Þ ~ ui;1 dn1 n2 ffi 2pk1 g1 xk21 1 cotðk1 pÞ;

ð39Þ

Se

r3 ðx2  n2 Þ~u_ i;2 dn1 n2 ffi 2g2 xk22 cotðk2 pÞ:

ð40Þ

Using above relationships, from Eqs. (33)–(35), the singular index can be determined by cotðk1 pÞ ¼ 0;

cotðk2 pÞ ¼ 0;

cotðk3 pÞ ¼ 0;

cotðk4 pÞ ¼ 0;

cotðk5 pÞ ¼ 0;

cotðk6 pÞ ¼ 0:

ð41Þ

Based on the method [64], an incremental local path independent vector integrals J_ gi of relevance to electromagnetic thermoelastoplastic creep fracture mechanics valid for arbitrary extend incremental loading/unloading histories is derived, Z Z J_ gi ¼ ½ðrij e_ ij þ 0:5r_ ij e_ ij Þni  ðp_ j uj;i þ pj u_ j;i þ p_ j u_ j;i ÞdS  ½ðrij e_ ij þ 0:5r_ ij e_ ij Þ;i  ðrjk þ r_ jk Þ_ejk;i  r_ jk ejk;i dA S As Z ð42Þ  ½ðr_ j3 uj;i þ rj3 u_ j;i þ r_ j3 u_ j;i Þ;3 þ ðbj þ b_ j Þu_ j;i þ b_ j uj;i dA: As

Considering relationships Eqs. (37)–(40), for a point p near the crack front, the following relationships can be obtained by using the time-domain main-part analytical method:

x3

x′1

x3

x1

pˆ

x′2 As

rθ qˆ0

Sε

x2 Fig. 2. A small semi-circle domain Se on the crack surface.

1024

B.J. Zhu, T.Y. Qin / Applied Mathematical Modelling 33 (2009) 1014–1041

Z

0:5

2

2 R3 u_ i;1 dn1 n2 ffi ðrr0 Þ 0 ð1  3R0 ðx2  n2 Þ Þ~

pg1 cos h_ 0 ;

ð43Þ

se

Z Z

2 2 2 _ ui;1 dn1 n2 ffi ðrr0 Þ0:5 pg1 cos h_ 0 ; R3 0 ð1  3R0 v0 x3 Þ~

ð44Þ

0:5 3R5 u_ i;3 dn1 n2 ffi ðrri Þ pg3 sin h_ j ; j vj x3 ðx2  n2 Þ~

ð45Þ

se

se

Z Z

0:5

4pg6 cos h_ j ;

ð46Þ

2 u_ i;6 dn1 n2 ffi ðrri Þ1:5 pg6 sin1 hi ð2 cos1 h_ j  3 cos h_ j  cos 5h_ j Þ: 30R7 j x2 x3 ~

ð47Þ

3R5 u_ i;6 dn1 n2 ffi ðrri Þ j x2 x3 ~

se

se

Other parameters are listed in Appendix A. Using relationships Eqs. (43)–(47), the extended singular incremental stress gradient field around the crack front can be expressed as follows: r_ 13 ¼ 0:25c44 v0 g1 ðrr0 Þ0:5 cos h_ 0 ; 5 P r_ 23 ¼ pðrri Þ0:5 ðg6 k33 ð15k33 v2i cot hi ð3 cos 0:25hi  3 sin2 0:5hi cos 0:5hi  0:25 sin2 hi cos 1:5hi Þ

ð48Þ

i¼1 1 hi ð2 cos1 0:5hi  3 cos 0:5hi  cos 2:5hi ÞÞ 2k33 ðrri Þ2 sin1 0:5hi þ 0:5k32 r1 r1 i sin   5 P ^im gm þ w ^ 2 e_ njK cos 0:5h0 Þ; þqi2 p cos 0:5hi t2 g2 þ 4 q m¼3

_3n ¼ r





5 P 0:5 ^i2 g2 cos1 h_ i þ 6 cot h_ i ^gimm þ g6 v2i ji5 j1 qin vi p cot hi ððrri Þ0:5 q cos h_ i Þ q i1 ðrr i Þ i¼1 m¼3   5 P 0:5 _ _ ^im gm þ ðrri Þ0:5 v3i k33 ji5 j1 þ 2pAi4 i33 ððrri Þ cos hi g2 cot hi þ q i1 g6 cos hi Þ

5 P

ð49Þ

ð50Þ

m¼3

þ

5 P

^ n e_ njK cos h_ 0 ; pðrri Þ0:5 w

n ¼ 3–5;

i¼1

P P 2 _ ^im gm þ 4k32 q ^i2 g2 v3i ji5 j1 #_ 3 ¼ Ai4 pðrr0 Þ0:5 ðcos h_ i ðk32 þ 4k33 v2i Þ sin1 h_ i Þ q i1 ð2=3k32 cos hi þ k33 ðrr 0 Þ 5

5

i¼1

m¼3

_ þ15rri g6 k32 v3i ji5 j1 i1 cot hi ð3 cos 0:5hi 1 1 þ k32 ðrr0 Þ sin hi ÞÞ

 _3n ¼ r_ 33 where r

D_ 3

 6 cos h_ i sin2 h_ i  0:25 sin2 hi cos 3h_ i Þð2v2i k33

ð51Þ

 B_ 3 .

5. Numerical procedure A method proposed in [61,65] can be generalized to solve the hypersingular in integral equations (33)–(35), numerically. Using the behavior near the crack front, the extended incremental displacement discontinuities gradient unknowns functions can be written as e_ i;J ðn1 ; n2 ; tÞ ¼ F iJ ðn1 ; n2 ; tÞnkJ W ðn1 ; n2 ; tÞ: U 2

ð52Þ

Other parameters are listed in Appendix A, aiJmn(t) are unknown constants. Substitute Eq. (52) into Eqs. (33)– (35), a set of algebraic equations for unknown aiJmn(t) can be obtained: S X H X s¼0

alJsh ðtÞI ilJsh ðx1 ; x2 Þ ¼ p_ J ðx1 ; x2 ; tÞ;

h¼0

where I ilJsh ðx1 ; x2 Þ are given in Appendix B.

ð53Þ

B.J. Zhu, T.Y. Qin / Applied Mathematical Modelling 33 (2009) 1014–1041

1025

6. Numerical solutions and discussion In order to verify the above method and illustrate its application, numerical calculations are performed for a 3D rectangular crack growth problem in EMTE-CVP-MCs. As shown in Fig. 3, the width of the crack is a and the length of the crack is b, the material constants we used are listed in Table 1. In order to facility the computing and comparing, we use non-dimensional quantities as follows: cij eij ; eij ¼ pffiffiffiffiffiffiffiffiffiffiffiffi ; c11 e11 c11 rffiffiffiffiffiffi rffiffiffiffiffiffiffi e11 l11   /; u ¼ u; /¼ c11 c11

2ij ¼

cij ¼

eij ; e11

d ij dij ¼ pffiffiffiffiffiffiffiffiffiffiffiffi ; c11 l11

Di  i ¼ pffiffiffiffiffiffiffiffiffiffiffi ; D c11 e11

d ij gij ¼ pffiffiffiffiffiffiffiffiffiffiffiffi ; e11 l11

ij ¼ l

lij l11

ð54Þ

Bi  i ¼ pffiffiffiffiffiffiffiffiffiffiffiffi B c11 l11

ð55Þ

. . . . σ 33∞ , D33∞ , B33∞ ,ϑ33∞ ⊗

⊗

⊗

⊗ . σ 31∞

⊗

σ32∞

x3(ξ3) S±

b P

r

x 2(ξ2)

ξ0

O

a

x1 (ξ3 )

.∞

σ 32

.

•

•

∞ σ31 •

•

•

.

.

.

.

∞ ∞ ∞ ∞ σ 33 , D33 , B33 ,ϑ33

Fig. 3. 3D rectangular crack in electro-magneto-thermo-elastic coupled viscoplatic multiple composites.

Table 1 Material constants for Eshelby’s tensors calculations [66] Elastic (GPa)

Magnetoelectric (N s/VC)

Piezoelectric (C m2)

c11

c12

c13

c33

c44

g11

g33

e15

e31

e33

286

173

170

269.5

45.3

0.005E9

0.003E9

11.6

4.4

18.6

Piezomagnetic (N/Am)

Magnetic (N s2 C2)

Dielectric (C2/N m2)

Thermal stress

d15

d31

d33

l11

l33

e11

e33

i11

i22

i33

550

580.3

699.7

590E6

157E6

0.08E9

0.093E9

0

0

0

Table 2 Non-dimensional material constant for Eshelby’s tensors calculations [67] c11

c12

c13

c33

c44

g11

g33

e11

e33

1

0.6049

0.5944

0.9423

0.1584

0.2301E4i

0.1381E4i

1

1.1625

d15

d31

d33

l11

l33

e15

e31

e33

0.0423i

0.0447i

0.0539i

1

0.2661

2.4251

0.9199

3.8885

1026

B.J. Zhu, T.Y. Qin / Applied Mathematical Modelling 33 (2009) 1014–1041

Substituting material constants in Table 1 into (Eqs. 54, 55), we obtain the value of non-dimensional material constant as Table 2.

-0.0007 -0.0006 -0.0005 -0.0004 -0.0003 -0.0002 -0.0001

0

0.0001 0.0002 0.0003 0.0004 0.0005 0.0006 0.0007 0.0008

Z

X

Y

1.00E-03 7.50E-04 5.00E-04 2.50E-04 . 0.00E+00 σ c -2.50E-04 -5.00E-04 -7.50E-04 -1.00E-03 -1 -0.75

0 0.25 0.5

-0.5 -0.25

0.75 1

x1 / a

0 1.25

0.25 1.5

0.5 1.75

x2 / b

0.75 2

1

(a) M=N=9

-8E-06 -6E-06 -4E-06 -2E-06

0

2E-06 4E-06 6E-06 8E-06 1E-05 Z

X

Y

0.0001 5E-05 0 -5E-05 -0.0001 -1 -0.75

0 0.25 0.5

-0.5 -0.25

0.75 1

x1 / a

0 1.25

0.25 1.5

0.5 1.75

x2 / b

0.75 2

1

(b) M=N=13 Fig. 4. Compliance of boundary condition (r_ c ¼ ðr_ 33 =r_ 1 33 Þ þ 1Þ.

.

σc

B.J. Zhu, T.Y. Qin / Applied Mathematical Modelling 33 (2009) 1014–1041

-0.0005 -0.0004 -0.0003 -0.0002 -0.0001

0

0.0001 0.0002

0.0003 0.0004

1027

0.0005 0.0006 Z

X

Y

1.00E-03 7.50E-04 5.00E-04 . 2.50E-04 σ c 0.00E+00 -2.50E-04 -5.00E-04 -7.50E-04 -1.00E-03 -1 -0.75

0 0.25 0.5 0.75 1 x1 / a 1.25

-0.5 -0.25 0 0.25 1.5 1.75

0.5

x2 / b

0.75 2

1

(a) M=N=9

-1.8E-05-1.6E-05 -1.4E-05 -1.2E-05 -1E-0 5 - 8E-06 -6E-06 -4E-06 -2E-0 6

0

2E-06

4E-06

6E-06

8E-06

1E-05 1.2E-05

Z

X

Y

8.00E-05 4.00E-05 0.00E+00 -4.00E-05 -8.00E-05 -1 -0.75

0 0.25

-0.5 -0.25

0.5 0.75

0

1

x1 / a

0.25

1.25

0.5

1.5

x2 / b

0.75

1.75 2

1

(b) M=N=13 Fig. 5. Compliance of boundary condition (r_ 31 =r_ 1 31 þ 1 ¼ 0Þ.

.

σc

1028

B.J. Zhu, T.Y. Qin / Applied Mathematical Modelling 33 (2009) 1014–1041

.

σc -0.95 -0.9 -0.85 -0.8 -0.75 -0.7 -0.65 -0.6 -0.55 -0.5 -0.45 -0.4 -0.35 -0.3 -0.25 -0.2 -0.15 -0.1 -0.05 Z

Y

X

0.5

0

-0.5 -1 -0.75

0 0.25 0.5

-0.5 -0.25

0.75

x1 / a1 1.25

0 0.25

1.5

x2 / b

0.5

1.75

0.75

(a) M=N=7

.

-0.0007 -0.0006 -0.0005 -0.0004 -0.0003 -0.0002 -0.0001

0

0.0001 0.0002 0.0003 0.0004 0.0005 0.0006 0.0007 0.0008

σc

Z

X

Y

0.5

0

-0.5 -1 -0.75

0 0.25 0.5

-0.5 -0.25

0.75

x1 / a

1

0 1.25

0.25 1.5

0.5 1.75

x2 / b

0.75

(b) M=N=9 _ 23I =r_ 21 Fig. 6. Extended remaining stresses rate on the crack surface r_ 13I =r_ 11 3I þ 1 ¼ 0 and r 3I þ 1 ¼ 0 when b/a = 1, KK = LL = 20  20.

B.J. Zhu, T.Y. Qin / Applied Mathematical Modelling 33 (2009) 1014–1041

1029

.

-8E-06 -6E-06 -4E-06 -2E-06

0

2E-06

4E-06 6E-06

8E-06 1E-05

σc Z

Y

X

0.5

0

-0.5 -1 -0.75

0 0.25 0.5

-0.5 -0.25

0.75

x1 / a

1

0

1.25

0.25

1.5

x2 / b

0.5

1.75

0.75

(c) M=N=11 .

σc 0

Z

Y

X

0.5

0

-0.5 -1 -0.75

0 0.25 0.5

-0.5 -0.25

0.75

x1 / a

1

0 1.25

0.25 1.5

0.5 1.75

0.75

(d) M=N=13 Fig. 6 (continued)

x2 / b

1030

B.J. Zhu, T.Y. Qin / Applied Mathematical Modelling 33 (2009) 1014–1041

6.1. Compliance of boundary condition and convergence of numerical solutions Consider a rectangular crack in EMTE-CVP-MCs subjected to the normal incremental mechanical load _1 _1 r_ 1 33 , the incremental electric load D33 , the incremental magnetic load B33 and the incremental thermal load _#1 in infinity. Figs. 4 and 5 show the compliance of the boundary condition along the crack surface for 33 b/a = 1, where the collocation point number is 400. It is found that the remaining incremental stress rate

σr

480 440 400

σr33

360

σr23

320 280 240 200 160 120 80 40 0

t 0

60

120

180

240

300

360

420

480

540

600

660

720

780

Fig. 7. Creep stresses redistribution at crack front as a function of remote tensile stress (b/a = 1, M = N = 13, KK = LL = 20  20).

120 110 100

σc

σc33 σc23

90 80 70 60 50 40 30 r

20

1.9902 1.9904 1.9906 1.9908 1.991 1.9912 1.9914 1.9916 1.9918 1.992

Fig. 8. Creep stresses concentration in the region very closed to crack front (b/a = 1, M = N = 13, KK = LL = 20  20).

B.J. Zhu, T.Y. Qin / Applied Mathematical Modelling 33 (2009) 1014–1041

1031

_ c ¼ 5:03  105 when M = N = 9, less than r_ c ¼ ðr_ 33 =r_ 1 33 Þ þ 1 on the crack surface is less than r 6 r_ c ¼ 1:04  10 when M = N = 13, Fig. 5 shows the remaining stress r_ c ¼ ðr_ 31 =r_ 1 31 Þ þ 1 on the crack surface is less than r_ c ¼ 7:05  105 when M = N = 9, less than r_ c ¼ 1:46  106 when M = N = 13. From the results in Fig. 6, it can be shown that the extended remaining stresses rate on the collocation points (KK = LL = 20  20,) possess a stable value at the flaw surface (the central region), though they increase sharply at the corner points (the boundary line). The present numerical method for multiple 3D flaws is stable and convincing. 6.2. Stresses and extended incremental displacement In case of b/a = 1, the number of the total nodes is taken as 20  20.Creep stresses redistribution at crack front is shown in Fig. 7. From the observation, it can be concluded that the reduction of the creep stress concentration is higher as time is increasing and the result is good agreement with those obtained in[54,10]. The creep stresses concentration in a small region closed to the crack front under the effect of a remote tensile stress loading reaching at a value of 1 is shown in Fig. 8. It shows that the creep stress component increase with r within a short distance from the crack front and a substantial reduction hereafter. Fig. 9 shown the extend incremental displacement ~ u_ as a function of the plastic strain rate e_ n . The change rule is good agree with those obtained in [54,55]. 6.3. Solutions for general cases For general cases, the polynomial exponents are taken as M = N = 13, the collocation point number is 20  20, and time is 800 s for the following results. Figs. 10 and 11 give the creep stresses redistribution rr23 and rr33 at crack front on the interface varied with the ration of b/a for different of t. The results show that both rr23 and rr33 increase b/a and converge to a constant value when b/a = 5. This means that the numerical solutions vary only slightly when b/a P 5, so the creep stresses redistribution at the crack front for the case of b/a P 5 is close to those for the 2D case. Figs. 12 and 13 give the creep stresses concentration of rc23 and rc33 in the region very closed to crack front on the interface as a function of b, a, and r. From the . u~

0.28 0.26 0.24

. u~1

0.22

. u~2

0.2 0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0

ε 0.006

0.008

0.01

0.012

0.014

0.016

0.018

n

0.02

Fig. 9. The extended incremental displacement as a function of plastic strain rate (b/a = 1, M = N = 13, KK = LL = 20  20).

1032

B.J. Zhu, T.Y. Qin / Applied Mathematical Modelling 33 (2009) 1014–1041 σr23 350 325 300 b/a=1 b/a=2 b/a=3 b/a=5

275 250 225 200 175 150 125 100 75 50

t

25 0

100

200

300

400

500

600

700

800

σr23

b/ a

40 60 80 100 120 140 160 180 200 220 240 260 280 300 320

9 8 7 6 5 4 3 2 2

t (×10 )

1

0

1

2

3

4

5

6

7

8

Fig. 10. Creep stresses redistribution rr23 at crack front as a function of b, a and t (KK = LL = 20  20, M = N = 13).

results in Figs. 12 and 13, we also found that the creep concentration rc23 and rc33 not only depend on the r but also on the ration of b/a, and they vary more gently when b/a P 5. 7. Concluding remarks The modeling of 3D crack growth in EMTE-CVP-MCs under electro-magneto-thermo-elastic coupled incremental loads rate is investigated by a time-domain hypersingular integral equation method. The conclusions to be drawn are as follows:

B.J. Zhu, T.Y. Qin / Applied Mathematical Modelling 33 (2009) 1014–1041

1033

σ r33 550 500 b/a=1 b/a=2 b/a=3 b/a=5

450 400 350 300 250 200 150 100

t

50 0

b/ a

100

200

300

400

500

600

700

800

σ r33 60 90 120 150 180 210 240 270 300 330 360 390 420 450 480 510 540

9 8 7 6 5 4 3 2 2

t (×10 )

1

0

1

2

3

4

5

6

7

8

Fig. 11. Creep stresses redistribution rr33 at crack front as a function of b,a and t (KK = LL = 20  20, M = N = 13).

1. Using the method involving of time-domain finite-part integrals and the extended fundamental incremental displacement rate solutions, the problem is analyzed through a set of TD-HIEs coupled with time-domain boundary integral equations. The behaviors of extended stress singularities near the front of crack are obtained by the time-domain main-part analysis of two-dimensional TD-HIE method, and the extended singular orders are analyzed. The results show that the singularity near the crack front in EMTE-CVPMCs is similar to those in general homogenous materials. The singular incremental stresses and extended J_ gi integral near the crack front in EMTE-CVP-MCs are given.

1034

B.J. Zhu, T.Y. Qin / Applied Mathematical Modelling 33 (2009) 1014–1041

σc23 85 80 75

b/a=1 b/a=2 b/a=3 b/a=5

70 65 60 55 50 45 40 35 30

r

25 1.9902 1.9904 1.9906 1.9908

1.991

1.9912 1.9914 1.9916 1.9918

1.992

σc23 25

9

29.6154 34.2308 38.8462 43.4615 48.0769 52.6923 57.3077 61.9231 66.5385 71.1538 75.7692 80.3846

85

b/ a

8 7 6 5 4 3 2

r 1

1.9902 1.9904 1.9906 1.9908

1.991

1.9912 1.9914 1.9916 1.9918

1.992

Fig. 12. Creep stresses concentration rc23 in the region very closed to crack front as a function of b, a and r (KK = LL = 20  20, M = N = 13).

2. The TD-HIEs are translated into a set of algebraic equations. In the process, the singular integrals of various types are specially treated and the specific calculating formulae are carried out. The numerical calculation expressions of the extended crack incremental opening displacements of the crack front are obtained. Some numerical solution are calculated and compared with those obtained in [54,55,10]. It shows that the present results are satisfied. 3. The creep stresses redistribution increase b/a and converge to a constant value when b/a = 5. This means that the numerical solutions vary only slightly when b/a P 5, so the creep stresses redistribution at the crack front for the case of b/a P 5 is close to those for the 2D case. The creep stresses concentration in

B.J. Zhu, T.Y. Qin / Applied Mathematical Modelling 33 (2009) 1014–1041

1035

σc33 130 120 b/a=1 110

b/a=2 b/a=3

100

b/a=5

90 80 70 60 50

r

40 1.9902 1.9904 1.9906 1.9908

1.991

1.9912 1.9914 1.9916 1.9918

1.992

σc33 40

9

45

50

55

60

65

70

75

80

85

90

95 100 105 110 115 120 125 130

b/ a

8 7 6 5 4 3 2

r 1

1.9902 1.9904 1.9906 1.9908 1.991 1.9912 1.9914 1.9916 1.9918 1.992

Fig. 13. Creep stresses concentration rc33 in the region very closed to crack front as a function of b, a and r (KK = LL = 20  20, M = N = 13).

the region very closed to crack front on the interface not only depend on the r but also on the ration of b/a, and they vary more gently when b/a P 5. 4. In conclusion, whenever there is an electro-magneto-thermo-elastic coupled structure containing 3D crack growth problem, an analysis of the type described in this paper can be utilized in order to find the critical configurations under which the structure may be most vulnerable. In such cases, the strength predictions would be safer if the crack growth to be taken into account.

1036

B.J. Zhu, T.Y. Qin / Applied Mathematical Modelling 33 (2009) 1014–1041

Acknowledgement The authors wish to thank the reviewers for the useful comments and suggestions.

Appendix A ni ¼ c44 l24 di1 þ ðc11 l24 þ ðc244  l221 Þl24 þ ðc44 e11 þ l232 Þl34 þ ðc44 l11 þ l232 Þl44 þ 2ðc44 e15  l21 l31 Þl34 þ 2ðc44 d 15  l21 l32 Þl41 þ 2ðc44 g11 þ l31 l32 Þl43 Þdi2 þ c11 l24 di5 þ ðc11 ðc44 l24 þ e15 l34 þ g15 l41 Þ þ c44 ðc33 l25 þ e33 l35 þ g33 l45 Þ þ ðc11 c33 þ c244  l221 Þl11 þ l231 l22 þ l232 l33 þ ðc11 e33 þ c44 e15  2l21 l32 Þl12 þ ðc11 d 33 þ c44 d 15  2l21 l32 Þl13 þ 2l31 l32 l23 Þdi3 þ ðc44 l14 þ ðc11 c33  l221 Þl24 þ ðc1133 þ l231 Þl42 þ ðc11 l33 þ l232 Þl44 þ 2ððc11 e33  l21 l31 Þl12 þ ðc11 d 33  l21 l32 Þl34 þ ðc11 d 33  l21 l31 Þl41 þ ðc11 g33  l31 l32 Þl43 ÞÞdi4 ;

ðA:1Þ

n4l ¼ ði1 l14 þ i3 l21 l24 þ ði3 l31  13 l21 Þl24 þ ði3 l32  k3 l21 Þl41 þ 13 l31 l42 þ ð13 l32 þ k3 l31 Þl52 þ k3 l32 l53 Þdl1 þ ði1 ðc44 l24 þ e15 l34 þ d 15 l31 Þ  ði1 c33  i3 l34 Þl11  ði1 e33  i3 l31 þ d 3 l21 Þl12  ði1 d 33  i3 l32 þ k3 l21 Þl13 þ ð13 l32 þ k3 l31 Þl23 þ 13 l31 l22 þ k3 l32 l33 Þdl2 þ ði1 ðc44 l11 þ e15 l12 þ d 15 l13 Þ  ði1 c33  i3 l21 Þl25  ði1 e33  ib3 l31 þ 13 l21 Þl35  ði1 d 33  i3 l32 þ k3 l21 Þl45 þ ð13 l32 þ k3 l31 Þl52 þ 13 l31 l51 þ k3 l32 l33 Þdl3  i1 l15 dl4 ;

ðA:2Þ

n5l ¼ c44 ði3 l24  13 l34  k3 l41 Þd1l þ ðð13 c11  i1 l21 Þl41 þ c44 ði3 l11  13 l12  k3 l13 Þ þ ði3 l31 þ 13 l21 Þðl33 l31  g33 l32 Þ  ði3 l32 þ k3 l21 Þðg33 l32  e33 l32 Þ þ ð13 l32  k3 l31 Þðd 33 l32  e33 l12 ÞÞd2l þ ðði3 c11  i1 l21 Þl11  c44 ð13 c11 þ i1 l32 Þl12  ðk3 c11 þ i1 l32 Þl13  c44 ð13 l35  i3 l25  k3 l45 Þ þ ði3 l32 þ 13 l21 Þð11 l32  g11 l32 Þ  ði3 l32 þ k3 l21 Þðg11 l32  e11 l32 Þ þ ð13 l32  k3 l31 Þðd 15 l32  e15 l12 ÞÞd3l þ ðði3 c11  i1 l21 Þl25  ð13 c11 þ i1 l32 Þl35  ðk3 c11 þ i1 l32 Þl45 Þd4l ;

ðA:3Þ

n6l ¼ c44 ði3 l34 þ 13 l42 þ k3 l43 Þd1l þ ðði3 c11  i1 l21 Þl34 þ l42 ði3 l11  13 l12  k3 l13 Þ þ l43 ðk3 c11  i1 l21 Þ  c44 ði3 l12 þ 13 l22 þ k3 l23 Þ  ði3 l31 þ 13 l21 Þðl33 l21  g33 l32 Þ  ði3 l32 þ k3 l21 Þðg33 l21 þ e33 l32 Þ þ ð13 l32  k3 l31 Þðd 33 l32  e33 l12 ÞÞd2l þ ðði3 c11  i1 l21 Þl12  ð13 c11 þ i1 l42 Þl22  ðk3 c11 þ i1 l32 Þl23 þ c44 ð13 l42 þ i3 l35 þ k3 l43 Þ  ði3 l32 þ 13 l21 Þðl11 l21 þ d 15 l32 Þ  ði3 l32 þ k3 l21 Þðg11 l21 þ e15 l32 Þ þ ð13 l32  k3 l31 Þðd 15 l21  c15 l32 ÞÞd3l þ ðði3 c11  i1 l21 Þl35 þ ð13 c11 þ i1 l32 Þl42 þ ðk3 c11 þ i1 l32 Þl52 Þd4l ;

ðA:4Þ

n7l ¼ c44 ði3 l41 þ 13 l43 þ k3 l44 Þd1l þ ðði3 c11  i1 l21 Þl41 þ l43 ði1 l31 þ 13 c11 Þ þ l44 ðk3 c11 þ i1 l32 Þ þ c44 ði3 l13 þ 13 l23 þ k3 l33 Þ þ ði3 l31 þ 13 l21 Þðd 33 l21 þ g33 l32 Þ  ði3 l32 þ k3 l21 Þðe33 l21 þ e33 l31 Þ þ ð13 l32  k3 l31 Þðe33 l21  e33 l31 ÞÞd2l þ ðði3 c11  i1 l21 Þl13  ð13 c11 þ i1 l32 Þl23 þ ðk3 c11 þ i1 l32 Þl23 þ c44 ð13 l52 þ i3 l45 þ k3 l53 Þ  ði3 l31 þ 13 l21 Þðd 11 l21 þ g15 l31 Þ  ði3 l32 þ k3 l21 Þðe11 l21 þ e15 l31 Þ þ ð13 l32  k3 l31 Þðe15 l21  c44 l31 ÞÞd3l þ ðði3 c11  i1 l21 Þl45 þ ð13 c11 þ i1 l31 Þl52 þ ðk3 c11 þ i1 l32 Þl53 Þd4l ;

ðA:5Þ

B.J. Zhu, T.Y. Qin / Applied Mathematical Modelling 33 (2009) 1014–1041

1037

lIJ ¼ d11 ðe11 l33 þ e33 l11  2g11 g33 Þ þ d12 ðe33 l11 þ e15 l33  g11 d 33  g33 d 15 Þ þ d13 ðd 33 e11 þ d 15 e33  g11 e33  g33 e15 Þ þ d21 ðc13 þ c44 Þ þ d22 ðc44 l33 þ c33 l11 þ 2d 15 d 33 Þ  d23 ðc44 g33 þ c33 g11 þ e15 d 33 þ e33 g15 Þ þ d31 ðe15 þ e31 Þ þ d33 ðc44 e33 þ c33 e11 þ 2e15 e33 Þ þ d32 ðd 15 þ d 31 Þ þ d14 ðg233 e33 þ ðe233 þ c33 e33 Þl33  2g33 e33 d 33  c33 d 233 Þ þ d42 ðc33 l33 þ d 233 Þ  d43 ðc33 g33 þ e33 d 33 Þ þ d24 ðl33 e33  g233 Þ þ d34 ðd 33 g33  e33 l33 Þ þ d41 ðe33 d 33  g33 e33 Þ þ d25 ðl11 e11  g211 Þ þ d35 ðd 15 g11  e15 l11 Þ þ d15 ðg233 e33 þ ðe233 þ c33 e33 Þl33  2d 15 e15 g11  c44 g211 Þ þ d51 ðc44 l11 þ g215 Þ  d52 ðc44 g11 þ d 15 e15 Þ þ d44 ðe233  c33 e33 Þ þ d45 ðe15 g11  d 15 e11 Þ;

ðA:6Þ

1 wiJ ¼ di1 þ di2 ðc13 þ ðc33 ji2 þ e33 ji3 þ d 33 ji4 Þv2i j1 i1  i3 j55 di5 j51 Þ þ di3 ðe31 þ ðe33 ji2  e33 ji3 1 2 1 1  g33 ji4 Þv2i j1 i1  13 j55 di5 j51 Þ þ di4 ðd 31 þ ðd 33 ji2  g 33 ji3  l33 ji4 Þvi ji1  g3 j55 di5 j51 Þ 1 þ 2di5 ðc11  c66 þ ðc13 ji2 þ e31 ji3 þ g31 ji4 Þv2i j1 i1  i3 j55 di5 j51 Þ;

jIJ ¼ n4l v6i  n5l v4i þ n6l v2i  n7l ;

K IJ ¼

5 X

ðA:7Þ

^i2 þ dI1 dJ 6 k33 v2i q ^i1 þ dmI ðdJt qim q ^i1 þ daJ Ai4 q ^i2 qi1 þ d6J Ai4 v3i ji5 j1 ðdI1 dJ 2 qi2 q i1 k33 qi6 ÞÞ

i¼1

b IJ ¼ K

 dI1 dJ 1 c44 v0 =4p;

ðA:8Þ

5 X ^i2 þ dmI d6J v3i ji5 j1 ðdIa dJ b Ai4 k3b q i1 qim Þ;

ðA:9Þ

i¼1

^iJ ¼ di1 ðc244 Di vi þ c44 Ai1 þ e15 Ai2 þ d 15 Ai3 Þ þ di2 vi ðc44 þ j1 q i1 ðc44 ji2 þ e15 ji3 þ d 15 ji4 ÞÞ  di3 c13  di4 e31  di5 d 31 þ v2i j1 i1 ðdi3 ðc33 ji2 þ e33 ji3 þ d 33 ji4 þ i33 ji5 ÞÞ þ di4 ðe33 ji2i  e33 ji3  g33 ji4 þ 133 ji5 Þ þ di5 ððd 33 ji2  g33 ji3  l33 ji4 þ g33 ji5 ÞÞ;

ðA:10Þ

qiJ ¼ di1 ðd3m i33 þ d4m 133 þ d5m g33 Þ þ di2 ðc44 Di vi þ c44 Ai1 þ e15 Ai2 þ d 15 Ai3 Þ þ di3 ðvi ðc33 Ai1 þ e33 Ai2 þ d 33 Ai3 Þ  Di c13 Þ þ di4 ðvi ðe33 Ai1  e33 Ai2  g33 Ai3 Þ  Di e31 Þ þ di5 ðvi ðd 33 Ai1  g33 Ai2  l33 Ai3 Þ  Di d 31 Þ þ di6 ðd3m i33 þ d4m 133 þ d5m g33 þ d6m Þ;

ðA:11Þ

^

1 1 w ¼ 0:5Ai1 vi ðdaJ ð2c44 þ ðc44 ð1 þ ji2 j1 i1 Þ þ e15 ð1 þ ji3 ji1 Þ þ d 15 ð1 þ ji5 ji1 ÞÞ J

1 1 þ d3J c13 ð1 þ ji2 j1 i1 Þ þ d4J e31 ð1 þ ji2 ji1 Þ þ d5J d 31 ð1 þ ji2 ji1 ÞÞÞ;

ðA:12Þ

2 2 2 2 2 ^ J ¼ 1:5Ai1 j1 w i1 vi ðd3J ðe33 ðji2 si þ ji3 Þ þ d 33 ðji2 si þ ji4 Þ þ 2c33 ji2 Þ þ d4J vi ðe33 vi ðji2 vi þ ji3 Þ 2 2 2  ðg33 ðji2 v2i þ ji4 Þ  2e33 ji5 vi ji2 j1 i1  f33 ji4 ÞÞ þ d5J vi ðd 33 ðji2 vi þ ji3 Þ  l33 ðji2 vi þ ji4 Þ

 g33 ji5 ji2 j1 i1 vi  g33 ji5 ÞÞ; h_ i ¼ 0:5hi ;

ðA:13Þ

r0 ¼ ðcos2 h þ v20 sin2 hÞ0:5 ; 0:5

W ðn1 ; n2 Þ ¼ ðða2  n21 Þð2b  n2 ÞÞ ;

h0 ¼ tan1 ðv0 tan hÞ;

F iJ ðn1 ; n2 ; tÞ ¼

S X H X s¼0

h¼0

rj eihj ¼ cos h þ ivj sin h;

aiJsh ðtÞns1 nh2 :

ðA:14Þ

ðA:15Þ

1038

B.J. Zhu, T.Y. Qin / Applied Mathematical Modelling 33 (2009) 1014–1041

Appendix B Z I i12sh ¼ 3 r3 ðK 11 þ K 12 Þr;1 r;2 ns1 nk22 þh W dn1 dn2 S Z i 3 r ðK 11 ð1  3r2;2 Þ þ K 12 ð1  3r2;1 ÞÞns1 nk21 þh W dn1 dn2 I 11sh ¼ S Z r3 K mn ns1 nk2n þh W dn1 dn2 I imnsh ¼ S Z i I 21sh ¼  3r3 ðK 11 þ K 12 Þr;1 r;2 ns1 nk21 þh W dn1 dn2 S Z r3 ðK 11 ð1  3r2;1 Þ þ K 12 ð1  3r2;2 ÞÞns1 nk22 þh W dn1 dn2 I i22sh ¼ I i61sh ¼ I i66sh ¼

Z

S

Z

S

ðB:1Þ ðB:2Þ ðB:3Þ ðB:4Þ ðB:5Þ

b 11 ð1  3r2 Þ  K b 12 r;1 r;2 Þns nk21 þh W dn1 dn2 r3 ð K ;1 1

ðB:6Þ

3r4 K 66 k3a r;a ns1 nk26 þh W dn1 dn2

ðB:7Þ

b m6 k3a r;a ns nk26 þh W dn1 dn2 3r4 K 1

ðB:8Þ

b 12 ð1  3r2 Þ  3 K b 11 r;1 r;2 Þns nk22 þh W dn1 dn2 r3 ð K ;2 1

ðB:9Þ

S

I im6sh ¼ I i62sh ¼ I ia6sh ¼

Z Z

S

Z

S

3r4 K a6 r;a ns1 nk26 þh W dn1 dn2

ðB:10Þ

r2 K m1 r;a ns1 nk2a þh W dn1 dn2

ðB:11Þ

S

I imash

¼

Z

S

The integral equations ((B.1)–(B.11)), is hypersingular and must be treated before being numerically evaluated. Using the Taylor’s expansion and the local polar coordinates n1  x1 = r1cos h1, n2  x2 = r1sin h1 as shown in Fig. B.1, following relationships can be obtained: ( 0:5 0:5 ða2  x21 Þ  x1 r1 cos h1 =ða2  x21 Þ  Q1 r21 ; s ¼ 0; s 2 2 0:5 ðB:12Þ n1 ða  x1 Þ ¼ 2 2 2 2 0:5 xs1 ða2  x21 Þ0:5 þ xs1 þ Q2 r21 ; s > 0; 1 ðsa  ðs þ 1Þx1 Þr1 cos h1 =ða  x1 Þ nk2J þh ð2b  n2 Þ

0:5

¼ x2kJ þh ð2b  x2 Þ

0:5

þ ðx2kJ þh1 ð4bðkJ þ hÞ

ð2kJ þ 2h þ 1Þx2 Þr1 sin h1 þ QJ3 r21 Þ=2ð2b  x2 Þ

0:5

R (θ1 ) a

r1 ( x1 , x 2 )

(ξ 1 , ξ 2 ) θ1

x2 (ξ 2 ) a

2b

x1 (ξ1 )

Fig. B.1. Integral parameters in local coordinate.

ðB:13Þ

B.J. Zhu, T.Y. Qin / Applied Mathematical Modelling 33 (2009) 1014–1041

1039

in which Q1 ¼ a2 ðx1 þ n1 Þðða2  x21 Þ0:5 ðða2  x21 Þ0:5 þ ða2  n21 Þ0:5 Þðn1 ða2  x21 Þ0:5 þ x1 ða2  n21 Þ0:5 ÞÞ1 cos2 h1 ðB:14Þ ( 3=2 sþ2  ðs2 þ 1Þa2 xs1 þ sðs þ 1Þx1 Þ cos2 h1 =2ða2  x21 Þ ; n1 ¼ x1 ðsðs  1Þa4 xs2 1 Q2 ¼ ðB:15Þ 2 0:5 2 s 2 2 2 2 2 0:5 r1 ðn1 ða  n1 Þ  xs1 ða2  x21 Þ0:5  xs1 n1 –x1 1 ðsa  ðs þ 1Þx1 Þðn1  x1 Þ=ða  x1 Þ Þ; 8 2ða2  x21 Þ1:5 x2kJ þh2 sin2 h1 ð16ðkJ þ hÞðkJ þ h  1Þb2  8bðk þ nÞð2kJ þ 2h  1Þx2 þ ð4ðkJ þ hÞ2  1Þx22 Þ; > > > < n 2 ¼ x2 QJ3 ¼ ðB:16Þ kJ þh > r2 ð2b  n2 Þ0:5  x2kJ þh ð2b  x2 Þ0:5  ð2b  x2 Þ0:5 x2kJ þh1 ð2bðkJ þ hÞ  ðkJ þ h þ 0:5Þx2 Þðn2  x2 ÞÞ; > 1 ðn2 > : n2 –x2

Using relationships (B.12, B.13), the kernel of integral equations ((B.1)–(B.11)), can be written as follows: nm1 nk2J þh ðða2  n21 Þð2b  n2 ÞÞ

0:5

¼ DJ0 þ DJ1 r1 þ DJ2 r21 ;

ðB:17Þ

where DJ0 , DJ1 , and DJ2 are known functions, and can be obtained by Taylor’s expansion. Using the finite part integral method and relationships (B.17), the hypersingular integral ((B.1)–(B.11)) can be reduced as   Z 2p Z R 2 i 2 1 1 1 1 ðB:18Þ ðK 11 ð1  3 sin hÞ þ K 12 ð1  3 cos hÞÞ D1 ln R  R D0 þ D2 dr1 dh I 11sh ¼ 0 0   Z 2p Z R ðB:19Þ I i12sh ¼ 1:5ðK 11 þ K 12 Þ sin 2h D21 ln R  R1 D20 þ D22 dr1 dh 0 0   Z 2p Z R I i21sh ¼ ðB:20Þ 1:5ðK 11 þ K 12 Þ sin 2h D11 ln R  R1 D10 þ D12 dr1 dh 0 0   Z 2p Z R 2 i 2 2 1 2 2 I 22sh ¼ ðB:21Þ ðK 11 ð1  3 cos hÞ þ K 12 ð1  3 sin hÞÞ D1 ln R  R D0 þ D2 dr1 dh 0 0   Z 2p Z R b 11 ð1  3 cos2 hÞ  0:5 K b 12 sin 2hÞ D1 ln R  R1 D1 þ ðB:22Þ I i61sh ¼ ðK D12 dr1 dh 1 0 0 0   Z 2p Z R I ia6sh ¼ ðB:23Þ 3ðK 16 cos hd1a þ K 26 sin hd2a Þ D61 ln R  R1 D60 þ D62 dr1 dh 0 0   Z 2p Z R 2 b 12 ð1  3 sin2 hÞ  0:5 K b 11 sin 2hÞ D2 ln R  R1 D2 þ I i62sh ¼ ðB:24Þ ðK D dr 1 dh 1 0 2 0 0   Z 2p Z R i 6 1 6 6 ðB:25Þ I 66sh ¼ 3K 66 ðk31 cos h þ k32 sin hÞ D1 ln R  R D0 þ D2 dr1 dh 0 0   Z 2p Z R b m6 ðk31 cos h þ k32 sin hÞ D6 ln R  R1 D6 þ I im6sh ¼ ðB:26Þ 3K D62 dr1 dh 1 0 0 0   Z 2p Z R I imash ¼ ðB:27Þ K m1 cos h Da1 ln R  R1 Da0 þ Da2 dr1 dh 0 0   Z 2p Z R ðB:28Þ I imnsh ¼ K mn Dn1 ln R  R1 Dn0 þ Dn2 dr1 dh 0

0

References [1] P.M. Davis, Piezomagnetic computation of magnetic-anomalies due to ground loading by a man-made lake, Pure Appl. Geophys. 112 (1974) 811–819. [2] N.F. Jordan, A.C. Eringen, On the static nonlinear theory of electromagnetic thermoelastic solids—I, Int. J. Eng. Sci. 2 (1964) 59–95. [3] N.F. Jordan, A.C. Eringen, On the static nonlinear theory of electromagnetic thermoelastic solids—II, Int. J. Eng. Sci. 2 (1964) 97– 114. [4] M. Tinkham, Electromagnetic properties of superconductors, Rev. Mod. Phys. 46 (1974) 587–596.

1040

B.J. Zhu, T.Y. Qin / Applied Mathematical Modelling 33 (2009) 1014–1041

[5] J. Vandenboomgaard, A.M.J.G. Vanrun, J. Vansuchtelen, Magnetoelectricity in piezoelectric-magnetostrictive composites, Ferroelectrics 10 (1976) 295–298. [6] X.J. Zheng, L. Sun, A nonlinear constitutive model of magneto-thermo-mechanical coupling for giant magnetostrictive materials, J. Appl. Phys., 100 (2005) 063906. [7] J. Aboudi, J.D. Achenbach, Rapid mode-III crack propagation in a strip of viscoplastic work-hardening material, Int. J. Solids Struct. 17 (1981) 879–890. [8] J. Aboudi, J.D. Achenbach, Arrest of fast mode-i fracture in an elastic–viscoplastic transition zone, Eng. Fract. Mech. 18 (1983) 109– 119. [9] H. Riedel, Creep deformation at crack tips in elastic–viscoplastic solids, J. Mech. Phys. Solids 29 (1981) 35–49. [10] M. Morjaria, S. Mukherjee, Numerical analysis of planar, time-dependent inelastic deformation of plates with cracks by the boundary element method, Int. J. Solids Struct. 17 (1981) 127–143. [11] J.D. Achenbach, N. Nishimura, J.C. Sung, Crack-tip fields in a viscoplastic material, Int. J. Solids Struct. 23 (1987) 1035–1052. [12] H.D. Bui, K. Dang Van, E. de Langre, A simplified analysis of creep crack growth using local approach, Nucl. Eng. Des. 105 (1987) 147–153. [13] R. Krishna Kumar, R. Narasimhan, O. Prabhakar, Finite element analysis of dynamic crack growth under plane stress in a viscoplastic material, Eng. Fract. Mech. 37 (1990) 1071–1084. [14] K.C. Wu, J.H. Lai, Anti-plane steady state crack propagation in elastic–viscoplastic material, Theor. Appl. Fract. Mech. 13 (1990) 239–249. [15] Y.C. Gao, Further study on strain singularity behavior of moving cracks in elastic–viscoplastic materials, Theor. Appl. Fract. Mech. 14 (1990) 233–242. [16] Z.X. Feng, An unified treatment of the singular integrals in BEM model for 3-D elasto-plastic crack analysis, in: Proceeding of the Fourth Japan–China Symposium on the Boundary Element Methods, 1991. [17] S. Baveendra, P.K.T. Banerjee, Boundary element analysis of cracks in thermally stress planar structures, Int. J. Solids Struct. 29 (1992) 2301–2317. [18] J.C. Sung, J.Y. Liou, Transient crack growth in a viscoplastic solid, Eng. Fract. Mech. 46 (1993) 399–411. [19] Y.Y. Zhu, S. Cescotto, A fully coupled elasto-visco-plastic damage theory for anisotropic materials, Int. J. Solids Struct. 32 (1995) 1607–1641. [20] C. Leppin, P. Wriggers, Numerical simulation of rapid crack propagation in viscoplastic materials, Compu. Struct. 61 (1996) 1169– 1175. [21] T. Siegmund, A. Needleman, A numerical study of dynamic crack growth in elastic–viscoplastic solids, Int. J. Solids Struct. 34 (1997) 769–787. [22] X. Yan, Analysis for a crack emanating from a corner of a square hole in an infinite plate using the hybrid displacement discontinuity method, Appl. Math. Model. 28 (2004) 835–847. [23] X. Yan, Multiple crack fatigue growth modeling by displacement discontinuity method with crack-tip elements, Appl. Math. Model. 30 (2006) 489–508. [24] S.L. Wang, M. Zhou, Y.Y. Zhang, BEM for elasto-plastic analysis of 3-D crack problems, in: Proceeding of the Third Japan–China Symposium on the Boundary Element Methods, 1990, pp. 13–21. [25] A.P. Cisilino, M.H. Aliabadi, J.L. Otegui, A three-dimensional boundary element formulation for the elastoplastic analysis of cracked bodies, Int. J. Numer. Meth. Eng. 42 (1998) 237–256. [26] J.Z. Zhang, M.D. Halliday, P. Bowen, P. Poole, Three dimensional elastic–plastic finite element modelling of small fatigue crack growth under a single tensile overload, Eng. Fract. Mech. 63 (1999) 229–251. [27] L.H. Liang, Y. Liu, G.S. Jia, J. Yang, Elastoplastic and creep line spring models for surface flaws by BEM Int. J. Pressure Vessels Piping 76 (1999) 781–787. [28] D. Bigoni, D. Capuani, Green’s function for incremental nonlinear elasticity: shear bands and boundary integral formulation, J. Mech. Phys. Solids 50 (2002) 471–500. [29] J. Sladek, V. Sladek, Z.P. Bazant, Non-local boundary integral formulation for softening damage, Int. J. Numer. Meth. Eng. 57 (2003) 103–116. [30] G.E. Oleaga, On the dynamics of cracks in three dimensions, J. Mech. Phys. Solids 51 (2003) 169–185. [31] J.L. Swedlow, T.A. Cruse, Formulation of boundary integral equations for three-dimensional elasto-plastic flow, Int. J. Solids Struct. 7 (1971) 1673–1683. [32] P. Riccardella, An Implementation of the Boundary Integral Technique for Planar Problems of Elasticity and Elastoplasticity, Ph.D. Thesis, Carnegie-Mellon University, Pittsurgh, PA, 1973. [33] V. Kumar, S. Mukherjee, A boundary-integral equation formulation for time-dependent inelastic deformation. in metals, Int. J. Mech. Sci. 19 (1975) 713–724. [34] A. Mendelson, Boundary Integral Methods in Elasticity and Plasticity, NASA TN-D-7418, 1973. [35] S. Mukherjee, Corrected boundary-integral equations in planar thermoelastoplasticity, Int. J. Solids Struct. 13 (1977) 331–335. [36] H.D. Bui, Some remarks about the formulation of three-dimensional thermoelastoplastic problems by integral equations, Int. J. Solids Struct. 14 (1978) 935–939. [37] S. Mukherjee, V. Kumar, Numerical analysis of time-dependent inelastic deformation in metallic media using. The boundary-integral equation method, J. Appl. Mech. Trans. ASME 45 (1978) 785–790. [38] J. Telles, C. Brebbia, On the application of the boundary element method to plasticity, Appl. Math. Model. 3 (1979) 466–470. [39] C. Brebbia, J. Telles, L. Wrobel, Boundary Element Techniques, Springer-Verlag, 1984.

B.J. Zhu, T.Y. Qin / Applied Mathematical Modelling 33 (2009) 1014–1041

1041

[40] J.C.F. Telles, C.A. Brebbia, The boundary element method in plasticity, Appl. Math. Model. 5 (1981) 275–281. [41] P. Banerjee, T. Davies, Advanced implementation of boundary element methods for three-dimensional problems of elastoplasticity and viscoplasticity, Dev. Bound. Element Meth. 3 (1984) 1–26. [42] G. Dargush, P.K.F. Banerjee, Boundary element methods for three-dimensional thermoplasticity, Int. J. Solids Struct. 29 (1991) 549– 565. [43] O. Huber, R. Dallner, P. Partheym, G. Kuhn, Evation of the stress tensor in 3-D elastoplasticity by direct solving of hypersingular integrals, Int. J. Numer. Meth. Eng. 39 (1998) 2555–2573. [44] V. Sladek, J. Sladek, Displacement gradients in BEM formulation for small strain plasticity, Eng. Anal. Bound. Elem. 23 (1999) 471– 477. [45] A. Frangi, G. Maier, Dynamic elastic–plastic analysis by a symmetric Galerkin boundary element method with time-independent kernels, Comp. Meth. Appl. Mech. Eng. 171 (1999) 281–308. [46] Y. Liu, L.H. Hong, Q.C. Hong, H. Antes, Non-linear surface crack analysis by three dimensional boundary element with mixed boundary conditions, Eng. Fract. Mech. 63 (1999) 413–424. [47] X.-W. Gao, T.G. Davies, An effective boundary element algorithm for 2D and 3D elastoplastic problems, Int. J. Solids Struct. 37 (2000) 4987–5008. [48] G.D. Hatzigeorgiou, D.E. Beskos, Dynamic elastoplastic analysis of 3-D structures by the domain/boundary element method, Comp. Struct. 80 (2002) 339–347. [49] M. Brun, D. Bigoni, D. Capuani, A boundary element technique for incremental, non-linear elasticity Part II: bifurcation and shear bands, Comp. Meth. Appl. Mech. Eng. 192 (2003) 2481–2499. [50] M. Brun, D. Capuani, D. Bigoni, A boundary element technique for incremental, non-linear elasticity Part I: formulation, Comp. Meth. Appl. Mech. Eng. 192 (2003) 2461–2479. [51] D. Soares, J.A.M. Carrer, W.J. Mansur, Non-linear elastodynamic analysis by the BEM: an approach based on the iterative coupling of the D-BEM and TD-BEM formulations, Eng. Anal. Bound. Elem. 29 (2005) 761–774. [52] K. Bertoldi, M. Brun, D. Bigoni, A new boundary element technique without domain integrals for elastoplastic solids, Int. J. Numer. Meth. Eng. 64 (2005) 877–906. [53] T.Y. Qin, R.J. Tang, Finite-part integral and boundary-element method to solve embedded planar crack problems, Int. J. Fract. 60 (1993) 373–381. [54] C.P. Providakis, Viscoplastic BEM fracture analysis of creeping metallic cracked structures in plane stress using complex variable techniques, Eng. Fract. Mech. 70 (2003) 707–720. [55] Z.S.K.H. Li, Engineering treatment model for creep crack driving force estimation: CTOD in terms of d5, Eng. Fract. Mech. 68 (2001) 221–233. [56] D. Bigoni, D. Capuani, Time-harmonic Green’s function and boundary integral formulation for incremental nonlinear elasticity: dynamics of wave patterns and shear bands, J. Mech. Phys. Solids 53 (2005) 1163–1187. [57] W.Q. Chen, K.Y. Lee, H. Ding, General solution for transversely isotropic magneto-electro-thermo-elasticity and the potential theory method, Int. J. Eng. Sci. 42 (2004) 1361–1379. [58] H.J. Ding, A.M. Jiang, P.F. Hou, W.Q. Chen, Green’s functions for two-phase transversely isotropic magneto-electro-elastic media, Eng. Anal. Bound. Elem. 29 (2005) 551–561. [59] E. Pan, Three-dimensional Green’s functions in anisotropic magneto-electro-elastic bimaterials, Z. Angew. Math. Phys. 53 (2002) 815–838. [60] H.H. Zhu, Q.J. Chen, L.D. Yang, Boundary Element Method and its Application in Rock Soil Engineering, Tong ji University Publication, 1996, pp. 162–169. [61] T.Y. Qin, N.A. Noda, Analysis of a three-dimensional crack terminating at an interface using a hypersingular integral equation method, J. Appl. Mech. – Trans. ASME 69 (2002) 626–631. [62] T.Y. Qin, N.A. Noda, Application of hypersingular integral equation method to a three-dimensional crack in piezoelectric materials, JSME Int. J. Ser. A – Solid Mech. Mater. Eng. 47 (2004) 173–180. [63] T.Y. Qin, Y.S. Yu, N.A. Noda, Finite-part integral and boundary element method to solve three-dimensional crack problems in piezoelectric materials, Int. J. Solids Struct. 44 (2007) 4770–4783. [64] R. Zong, Boundary Element Method for Analysis of Thermoelastoplastic-creep Fracture [Dr. Dissertation], Shanhai Jiaotong University, 1996. 36–40. [65] T.Y. Qin, N.A. Noda, Stress intensity factors of a rectangular crack meeting a bimaterial interface, Int. J. Solids Struct. 40 (2003) 2473–2486. [66] J.Y. Li, Magnetoelectroelastic multi-inclusion and inhomogeneity problems and their applications in composite materials, Int. J. Eng. Sci. 38 (2000) 1993–2011. [67] B.J. Zhu, T.Y. Qin, Application of hypersingular integral equation method to three-dimensional crack in electromagnetothermoelastic multiphase composites, Int. J. Solids. Struct. 44 (2007) 5994–6012.