A hybrid finite element method for heterogeneous materials with randomly dispersed elastic inclusions

FINITE ELEMENTS "IN ANALYSIS A N D DESIGN

ELSEVIER

Finite Elements in Analysis and Design 19 (1995) 45-55

A hybrid finite element method for heterogeneous materials with randomly dispersed elastic inclusions J. Zhang, N. K a t s u b e Department of Engineering Mechanics, The Ohio State University, Columbus, OH 43210, USA

Abstract A new hybrid finite element method is presented for the mechanical analysis of heterogeneous materials with randomly dispersed inclusions. A special n-sided polygonal element with an elastic inclusion is developed on the basis of the Hellinger-Reissner principle. The element formulations are derived by decomposing the original problem into inclusion and matrix problems and relating each other through consistency conditions at their interface. For circular inclusions, the proposed method is verified against simple analytical solutions and shown to be suitable even for extreme cases such as porous materials and materials with rigid inclusions. It is also demonstrated that the effect of random packing of inclusions on stress concentration factors can be easily evaluated using the proposed method.

1. Introduction In the conventional finite element analysis of a heterogeneous material, a relatively complicated mesh [1] needs to be generated. F o r many heterogeneous materials, the inclusion phases often exhibit irregularities in geometry and randomness in distribution. This gives rise to great difficulty in generating the finite element mesh, especially when the number of the inclusion phases is large. G h o s h and M u k h o p a d h y a y [2] recently proposed a finite element model by generating an n-sided polygonal network for a two-dimensional heterogeneous medium. In their model, each polygon contains an inclusion and serves as an element. The complexity of generating a conventional mesh is thus avoided. The notion of transformation strain is introduced into their formulation as was proposed by Accorsi [3, 4]. However, it has been pointed out by Zhang and K a t s u b e [5] that the transfe, rmation strain does not play any role in this application. While generation of an n-sided network i,; indeed a notable contribution, some difficulties still exist in prediction of irregular stress distributions. Zhang and K a t s u b e [6-1 recently focused their attention on the problems where randomly dispersed inclusions are rigid and circular. The n-sided polygonal element by G h o s h and 0168-874X/95/$09.50 ~ 1995 Elsevier Science B.V. All rights reserved SSDI 0 1 6 8 - 8 7 4 X ( 9 4 ) 0 0 0 5 6 - 5

46

J. Zhang, N. Katsube / Finite Elements in Analysis and Design 19 (1995) 45-55

Mukhopadhyay is employed in order to deal with randomness. In numerically approximating irregular stress distribution, however, their approach is similar to the works by Tong et al. [7] and also by Piltner [8]. In its element formulation, classical elasticity solutions are directly used. Because of the simplification that the stress distribution only in the matrix needs to be analyzed, the resulting approach is very effective. Only one polygonal element with a rigid circular inclusion replaces many conventional elements in obtaining the same numerical accuracy. In this work, Zhang and Katsube's approach is extended to the cases where inclusions are also linearly elastic. The difficulties associated with approximating irregular stress distribution both in the matrix and the inclusion within an element are overcome based on the notion of superposition. The element formulations for the matrix and those for the inclusion region are performed separately and are later combined on the basis of consistency conditions at the interface. As observed in the limited case with rigid inclusions, only one polygonal element with an elastic circular inclusion (a hole or a rigid inclusion as extreme cases) is needed to replace many conventional elements in obtaining the same numerical accuracy.

2. Element decomposition and variational formulations Let us consider a special polygonal element embedded with an elastic inclusion of arbitrary geometry as shown in Fig. l(a). Our goal is to establish the relationship between the element's nodal forces and displacements, or simply, to formulate the element stiffness matrix. The key idea in formulating this special element is to decompose the original problem into two different boundary value problems as shown in Figs. l(b) and (c): (i) a specified boundary displacement problem in the matrix region t21 with boundaries C1 and C2; (ii) a specified boundary traction problem in the inclusion region ~'~2 with boundary C2. Cx and C2 are the element's outer boundary with neighboring elements and the inner interface between the matrix and the inclusion, respectively. Following the Hellinger-Reissner principle, we define the following two separate hybrid functionals for these two problems:

For problem i, ,, ~me = fa r±ar k2 m~mO'm --

a.Tm(DUm)-I dQ1 --

1

~ctTm(~(1)Um) --

-- (~ tTm(/I(2) -- l/m) dS 2, dc

dS1

1

(1)

2

For problem ii e = ; t ~ r-l-aT" 7~p 1_2 p dbpO'p - a p T ( D u p ) ] d Q 2

+~c

upTt(2) dS2,

(2)

2

where the subscripts m and p indicate the matrix and inclusion regions, respectively; the superscript T denotes the transpose of a vector or a matrix, a, u, t, and S are, respectively, the stress,

J. Zhang, N. Katsube / Finite Elements in Analysis and Design 19 (1995) 45-55

47

^(1)

m D

° S_j (a)

(b)

(c)

Fig. 1. Element decomposition: (a) original element; (b) matrix region with specified boundary displacements; (c) inclusion region with specified boundary traction.

displacement, and boundary traction vectors, and the elastic compliance matrix; D is the matrix differential operator relating strains to displacements; and the symbol -,~ represents a specified quantity. Based on the hybrid finite element method [9, 10], the boundary displacements ti m are assumed separately from Um and arc expressed in terms of the nodal displacements of the clement. The boundary displacements /~(2) in problem i and the boundary traction t~2) in problem ii are assumed to be specified and are actually unknown. The stationary values of the two functionals defined by Eqs. (1) and (2) yield the following two sets of equations:

Problem i SmO"m = Du m

DTO'm = 0,

t m = ~(1)O'm,

11m = li ")

on

in ~C'~l,

(3, 4)

Ci, i = 1, 2.

(5, 6)

Problem ii DTap = 0,

-

li<2)ap =

Sptrp =

Dup

in f2z,

(7, 8)

t~'.> on C2,

(9)

where ~i¢i) is a 2 x 3 matrix of the unit normal to boundary Ci (i = 1, 2). It is seen from the above equations that the true solutions of the two problems minimize the two hybrid functionals. Following the same argument as in the case for rigid inclusion [6], we may simplify the two functionals by constructing the stress ~ and displacement u fields in such a manner that Eqs. (3)-(5), (7), and (8) are automatically satisfied. By doing this and using the divergence theorem o v e r ~r~1 and f22, the two functionals are reduced to T tmUm dS1 +

= 1

tTum dS2

tTu tl) dS1 --

1

2

/Tm/~(2) dS2,

(10)

2

and

1~c tTUp dS2 + ~C

UpT~2)dS2,

(11)

J. Zhan9, N. Katsube/ Finite Elements in Analysis and Design 19 (1995) 45 55

48

where tp =

(12)

- - /~(2)~.p.

In order to recover the solution for the original problem from those of the two decomposed problems, the following consistency conditions need to be imposed: /~(2) = Up ~2) =

C2,

on

__ tin

on

(13)

C2.

(14)

Thus, Eqs. (10)-(14) consist of the fundamental variational formulations for the special elastic inclusion-embedded polygonal element.

3. Construction of special approximating functions The performance of an element largely depends on the selection of approximating functions for the element. In the conventional displacement-based or hybrid finite element methods, approximating functions for displacements and stresses are constructed in the form of a complete polynomial of a certain degree. As long as the displacement and stress fields are regular (smooth) in the element and the size of the element is sufficiently small, a good approximation can be achieved by assuming such a polynomial. In the present case, however, the presence of the inclusion severely perturbs an otherwise regular distribution of the stresses and displacements in the matrix. Moreover, the element mesh is not allowed to be refined. Thus, one has to search for special approximating functions such that the actual stress and displacement fields in the element can be properly characterized. In addition, in order to satisfy the requirements for using the two hybrid functionals of Eqs. (10) and (11) mentioned in the previous section, special approximating functions are needed such that the assumed stresses and displacements are able to satisfy a priori both the equilibrium equations (3) and (7) and the constitutive equations (4) and (8). In Zhang and Katsube [6], two special approximating functions were constructed for stresses and displacements based on classical elasticity solutions [11] for a multiply connected circular region. The two holomorphic functions in the general expressions of the elasticity solutions are expanded into two complex Laurent series respectively. The resulting approximating functions, which are suitable for the matrix region f21 in this case, are summarized as follows:

{ O-~nl)

mu

o.m2J k=mb

U~} =

u~

1

~

2A1

-- A 3

-

2A2 + A4

2A1

+ A 3

-

2A2 -

A4

[KB1 -- n 3

2Gm k~m~[_tcB2 + B4

A4

A3

-- l¢B2 + B4

KB1 + B3

-

Aa

A2

A1

-- A2

A2

A1

B1 B2

--

or

o" m

= em#m

(15)

br~)

or

Um= Um~m

(16)

49

J. Zhang, N. Katsube/ Finite Elements in Analysis and Design 19 (1995) 45-55

where Al(k,r,O)

= krk-Xcos(k

A3(k,r,O)

= k(k -

Bl(k,r,O)

= rkcoskO,

B3(k,r,O)

= krkcos(k

--

1)8,

1)rk-lcos(k

A2(k,r,O)

3)8,

--

B2(k,r,O) -

2)0,

= kr k-1

sin(k - 1)0,

A4(k,r,O) = k(k -

= krksin(k

--

3)0,

(19,20) (21,22)

= rksinkO,

B4(k,r,O)

1)rk-Xsin(k

(17, 18)

-

2)0.

(23, 24)

In the above expressions, x = (3 - Vm)/(1 q- Vm) for plane stress and x = 3 - 4vm for plane strain; Gm and Vmare the shear modulus and Poisson's ratio, respectively; r and 0 are polar coordinates with the origin located inside the inner boundary; m, and mb are the upper and lower limits of the Laurent series and mb is generally set to be - m, for symmetry; and 8k, ak, b'k, and /~k are real unknown coefficients to be determined. The superscripts and subscripts m in the expressions of the stress, displacement, and coefficient components indicate the matrix. As discussed in [6], the four coefficients 80, 40, b-o, and /~o are zero due to the necessity of eliminating the rigid body movement in the above expressions. In addition, the coefficient/~_ 1 has to be zero to guarantee zero resultant moment exerted on the inclusion. The expressions of Eqs. (15) and (16) can also be used as approximating functions for the elastic inclusion region ~22 if the terms involving negative k are excluded. Thus, by setting the value of mb to be one, the following approximating functions similar to those for the matrix region are obtained for the inclusion region: ~p = e p ~ p ,

(25)

Up = Up~p.

(26)

4. Element stiffness matrix

In this section, the stiffness matrix of the special polygonal element is derived by incorporating the results in the previous two sections. With Eqs. (5), (12), (15), and (25), the boundary traction on Ca and C2 can be expressed as follows: t m = (~(i)pm)flm

= R(im)flm

tp ~- ( -- ~(2)jl~p)~p =

on Ci,

i = 1, 2

R(2)~p o n C 2.

(27) (28)

By substituting Eqs. (16), (26), (27), and (28) into Eqs. (10) and (11), the two hybrid functionals become 1 ~T i'/..~(1) q'- H ( 2 ) ) ~ m - /l~e = 2 1 =

-

T

i #pHp#.

T + #,,Dp,

~TGmq

T

(29)

(30)

J. Zhan#, N. Katsube/ Finite Elements in Analysis and Design 19 (1995) 45-55

50

where 1 ~c [(R(ml))TUm "1"-~/'fT H(m1) = "~ m -R(1)-I ' m J dS1,

H(2) = ~1 ~c [(R(m2))TUm +

1

Up = ~1 ~c

~/ fmT=~- (m2 ) - ]j

dS2,

(31, 32)

2

F(R(2)) T

Up -'~ U T ; e p(2)] dS2,

Gm = ~c, (R~))TL dS1,

(33, 34)

2

(35, 36) L is the matrix relating the boundary displacements/~(1) to the nodal displacements q. According to the hybrid finite element method, L is assumed separately from the approximating functions inside the element in order to satisfy the inter-element displacement compatibility. For the specified boundary displacement problem i, the stationary value of n~mof Eq. (29) with respect to//m gives (H(m1) -1- H(rn2))flm --

Gmq -Dm = 0.

(37)

Similarly, for the specified boundary traction problem ii, we have -

np~p --~Dp =

0.

(38)

Using the consistency conditions of Eqs. (13) and (14), we may solve Eqs. (37) and (38) simultaneously for the unknown coefficients ~m and #p in terms of the nodal displacements q as follows:

tim

=

(H~) + H~ ) + H m p ) - l G m q ,

/]p = _ H p l

(39)

T Gmpflm,

(40)

Hm p• = GrapHp 1Gmp, T

(41)

Gmp = ~_ (e(m2))TUpdS2. .it:

(42)

where

2

As we know, within the linearly elastic range, an element stiffness matrix can be obtained by expressing the strain energy of the element in terms of the element's nodal displacements. For the special polygonal element, the strain energy Vst is expressed as Vst -- ½ # Tm ( n (11) -}- H(m2))flm+ ½~TUp#p"

(43)

The substitution of Eqs. (39) and (40) into Eq. (43) leads to Vs t = 12Pm[a"tm aT/~t)

+ H~) + Hmp)ffm = laT/~T [/.I(1) "Jr-U(m2) + Ump)- 1Gmq. z"z ~m~.==m

(44)

J. Zhang, N. Katsube/ Finite Elements in Analysis and Design 19 (1995) 45-55

51

0.5"

Co

"''[1''"

Fig. 2. A large plate with a small circular inclusion subjected to uniform tension.

Then, the stiffness matrix of the special element is T 1} Ke -- Gm(H(m 1 -4- H(2) -4- Hmp)- 1Gin.

(45)

It is seen from the above equations that the element stiffness matrix Ke can be numerically evaluated by simple matrix operations and numerical integration only along two boundaries. Because of the elimination of area integration, the proposed method has great advantage in dealing with high-order approximating functions. The stresses and displacements in the element can be calculated from Eqs. (15), (16), (25), and (26) once the unknown coefficients tim and pp are determined.

5. Numerical examples

5.1. Verification against simple analytical solutions As shown in Fiig. 2, an elastic circular inclusion of radius 0.25" is embedded at the center of a 4"× 4" square plate. The plate is subjected to a remote uniform tension ao in the vertical direction. The material properties of the plate and the inclusion, respectively, are: Young's moduli E m = 2.5 x 105 psi and E p = k x Era; and shear moduli G m : 0.96 x 105 psi and Gp : k x Gin. k is the parameter representing the ratio of the elastic moduli of the plate and the inclusion materials. In the finite element model, the plate is partitioned into 60 conventional 4-node and one special

J. Zhang, N. Katsube/ Finite Elements in Analysis and Design 19 (1995) 45 55

52

3 H Y B R I D (61 elements) :

k=

2.5

= Gp / G m ¢1 r~ t~

"0 O)

: k = 0.01

/

Ep/E m

2 1.5

•

: k=0.5

i

s

: k=l.0

I

~

: k=2.0

i

o

: k = 100.0

~"= ~

Analytical: - -

.N

11 0

"~l

z

0.5

~

3~xx / ~o

i

0 -0.5~

I

-1.5

I

-1

l

-0.5 0 0.5 1 X coordinateat cross section A-A

1.5

2

Fig. 3. N o r m a l i z e d stresses c%/ao a n d 3axx/ao at cross section A-A for various cases.

8-node elements. Doubly linear interpretation is used for the boundary displacement approximation on each side of the special element to satisfy the inter-element displacement compatibility. The value of the upper limit of the Laurent series m, is set to be five in accordance with the criterion discussed in [6]. The 4-node elements can either be conventional displacement-based quadrilateral elements or homogeneous hybrid ones. Fig. 3 shows the normalized stresses ayr/aO and 3axx/ao at cross section A - A for five cases with different values of k. The analytical results are shown in solid line for each case. A small value of k represents the case where the inclusion approaches a hole while a large value of k represents the case where the inclusion approaches a rigid one. In all five cases, excellent agreement is observed for both 6ry and axx between the numerical results obtained by the proposed method and the analytical predictions. This implies that the specially formulated element can be used for any inhomogeneous material including porous materials and materials with rigid inclusions. The efficiency of the proposed special element is clearly demonstrated since many conventional elements must otherwise be used to model the presence of one inclusion.

5.2. Effect of random packing on stress concentration factors A unit cell model is often used in the micromechanical analysis of fibrous composites. This is based on the assumption that the fibers are perfectly spaced and periodically arranged in the matrix such as square in Fig. 4(a) and hexagon in Fig. 4(b). The unit cell method works reasonably well for predicting some microstructure-insensitive properties such as elastic moduli in the axial direction. However, many other physical quantities including elastic moduli in transverse directions are

53

J. Zhang, N. Katsube / Finite Elements in Analysis and Design 19 (1995) 45-55

0000 0000 0000 0000 (a)

---.

??.. . : . ) •

•

uO|

",'.','))..-.I"'"' (b)

(c)

Fig. 4. Three fiber packing patterns: (a) square; (b) hexagonal; (c) random.

1~o

1~o

Fig. 5. Cross section of a Ti/SCS-9 composite specimen with its n-sided polygonal mesh.

microstructure sensitive. Failure mechanisms such as microcrack initiation and propagation also depend on the microstructure. Since the distribution of fibers in a real composite is often random as shown in Fig. 4(c), it is important to examine the effect of random packing on the local stress concentration factors. In this example, a unidirectional Ti/SCS-9 composite specimen 1-12-] is analyzed using the proposed finite element method. A portion (0.15" x 0.04") of the cross section of the specimen is shown in Fig. 5 with the generated polygonal mesh. The specimen is subjected to a uniformly transverse tension ao and a plane strain deformation model is used in the analysis. The material properties of the composite are: E m = 12.3 × 106 psi; Vm = 0.32; Ep = 62.0 × 106 psi; and Vp = 0.30. The generated polygonal mesh consists of 24 four-sided and 41 five-sided elements. There are three nodal points on each side of the elements and quadratic or doubly linear interpolation can be used for the displacement approximation on each side. The value of m, is set to be five for the four-sided elements and six for the five-sided ones. Fig. 6 shows the normalized radial, circumferential, maximum principal, and maximum shear stresses around includion # A. These numerical results are useful in determining the initiation of matrix cracking and matrix-fiber debonding. In Fig. 7, a comparison is made among the maximum principal stresses for three different fiber packing patterns: (a) random distribution (inclusion # A in Fig. 5), (b) square distribution, and (c) hexagonal distribution. The peak value of the principal stresses obtained from the random (real) distribution is about 42% higher than those based on periodic distribution assumptions.

J. Zhang, N. Katsube/ Finite Elements in Analysis and Design 19 (1995) 45-55

54

3

,

,

,

,

,

,

,

.N ¢0

E. 1

o Z

-2

: radial stress • : hoop stress : maximum principal stress o : maximum s h e a r stress •

_.~I

I

I

-0

50

100

I

principal direction

I

150 200 Angle (degree)

!

I

I

250

300

350

Fig. 6. Normalized stresses around inclusion # A.

i

i

i

i

i

i

i

• : random(Fig. 5)

2.51

= :

square

2

~

1.5

F-,

1

~0.5

o -0.5 10

10

I

100

I

I

150 200 Angle (degree)

I

I

I

250

300

350

Fig. 7. Comparison of maximum principal stresses for three different fiber distributions.

6. C o n c l u s i o n s

In this work, a new hybrid finite element is presented for the mechanical analysis of heterogeneous materials with elastic matrix and randomly dispersed elastic inclusions. An n-sided polygonal mesh is used to model the heterogeneity and a special polygonal element containing an elastic inclusion is developed. The element formulations are derived by separating the inclusion from the

J. Zhang, N. Katsube / Finite Elements in Analysis and Design 19 (1995) 45-55

55

matrix and using two different hybrid functionals for the matrix and the inclusion, respectively. Consistency conditions are imposed at the interface of the matrix and inclusion to recover the solution for the original element. For circular inclusions, special approximating functions are constructed based on the classical elasticity solutions for a multiply connected circular region. The resulting element stiffness matrix only involves integration over the element's outer and inner boundaries. This makes numerical evaluation of the element stiffness matrix relatively easy. The proposed method is verified against analytical solutions of problems with one circular inclusion in the infinite (large) plate. By using one special element with an inclusion, the analytical stress distribution surrounding the inclusion is reproduced based on the proposed method. Excellent agreement is obtained even in extreme cases where the inclusion becomes rigid or a hole. Versatility of the proposed method is also demonstrated by examining the effect of random packing on the local stress concentration factors.

Acknowledgments This work was supported by the National Science Foundation under the grant No. MSS9009347.

References [1] T. Nakamura and S. Suresh, "Effects of thermal residual stresses and fiber packing on deformation of metal-matrix composites", Acta Metall. Mater. 41 (6), pp. 1665-1681, 1993. I'2] S. Ghosh and S.N. Mukhopadhyay, "A material based finite element analysis of heterogeneous media involving Dirichlet tessellations', Comput. Methods Appl. Mech. Eng. 104, pp. 211-247, 1993. 1'3] M.L. Accorsi, "A method for modeling microstructural material discontinuities in the finite element analysis", Int. J. Numer. Methods Eng. 26, pp. 2187-2197, 1988. [41 M.L. Accorsi and R. Chamarajanagar, "Numerical validation of a hybrid finite element method using eigenstrain", Comput. Struct. 41, pp. 1065-1071, 1990. 1'5-1 J. Zhang and N. Katsube, "Problems related to application of transformation strains in a finite element analysis", Int. J. Numer. Methods Eng. 37, pp. 3185-3193, 1994. 1'6] J. Zhang and N. Katsube, "A hybrid finite element method for heterogeneous materials with randomly dispersed rigid inclusions", Int. J. Numer. Methods Eng. (to appear). [7] P. Tong, T.H.H. Pian and S.J. Lasry, "A hybrid-element approach to crack problems in plane elasticity", Int. J. Numer. Methods Eng. 7, pp. 297-308, 1973. [8] R. Piltner, "Special finite elements with holes and internal cracks", Int. J. Numer. Methods Eng. 21, pp. 1471-1485, 1985. [9] T.H.H. Pian, "Derivation of element stiffness matrices by assumed stress distributions", AIAA J. 2, pp. 1333-1336, 1964. [10] P. Tong and T.H.H. Pian, "A variational principle and the convergence of a finite-element method based on assumed stress distribution", Int. J. Solids Struct. 5, pp. 463-472, 1969. [11] N.I. Muskhelishvili, Some Basic Problems of the Mathematical Theory of Elasticity, P. Noordhoff, Groningen, Holland, pp. 120-129, 1953. [12] B.M. Rabeeh and W.O. Soboyejo, unpublished research work.