Strain singular field near the edge of interface in material with cuboidal inclusion based on experiment and numerical analysis

(1)

Strain singular field near the edge of interface in material with cuboidal inclusion based on experiment and numerical analysis

Makoto KASAI

¹

, Hideo KOGUCHI

²

, Takahiko KURAHASHI

²

1

Graduate school of Nagaoka University of Technology, 1603-1 Kamitomiokamachi, Nagaoka, Niigata

2

Department of Mechanical Engineering , Nagaoka University of Technology, 1603-1 Kamitomiokamachi, Nagaoka, Niigata [email protected]

Abstract

We present strain singular field near the interface edge in three-dimensional joint and experimental and numerical study are carried out. Dissimilar material joints are used in various fields. Singularity fields for stress and strain occur at vertex on interface in dissimilar material joints, when external load is applied to the dissimilar material joints, and it might be induced delamination or crack. Therefore, it is necessary to investigate singularity field near the interface edge. In this study, a coupon of specimen which silicon chip was embedded in resin is employed, and tensile load is applied to the specimen. Surface shape of specimen is measured by digital microscope before and after loadings.

The Digital Image Correlation Method (DICM) using surface shape on specimen and cross correlation coefficients for surface pattern was employed for evaluating displacement on the surface. Strain on surface of specimen is computed by using the Moving Least Square Method (MLSM). On the other hand, the Element Free Galerkin Method (EFGM) is applied to compute the displacement and strain distribution in the three-dimensional model of specimen. In this study, the specimen model employed in experiment is employed as computational model, and strain distribution near the edge of interface is computed based on the EFGM. Finally, discussions of significant digit for computation of strain are carried out based on results by the EFGM. In addition, results of displacement and strain distributions by the DICM are shown.

Discussion 1 Introduction

When external load act on dissimilar material joints, it is well known that stress and strain singularity field occurs near the edge of interface in dissimilar material joint. [1] In this study, we discuss strain field near the edge of interface in three-dimensional joints by experiment and numerical analysis.

In this study, tensile load is applied to specimen, and surface image of specimen is measured by digital microscope.

Light values of each pixel are obtained from surface image.

Displacement on surface of specimen is computed by using the DICM for light value of each pixel. [2], [3], [4] In the DICM, displacement is computed by the comparison of cross-correlation coefficients at each subset. The DICM is consists of rough search and fine search. In rough search, value of rigid body motion at each subset is computed. In fine search, value of subset deformation is computed. Strain is computed by using the MLSM for displacement. The MLSM can compute strain from distributed displacement.

Advantage of this method is not necessary to prepare elements for computation of strain.

On the other hand, the EFGM is applied to compute the displacement and strain distribution in the three-dimensional model of specimen. [5][6] The BEM and the FEM are frequently employed for stress analysis for three-dimensional joints. However, it is difficult to obtain surface strain distribution by the BEM. Furthermore, when the FEM is employed to obtain the strain and stress distribution, it is necessary to take care for mesh generation around interface considering mesh connectivity condition. In the EFGM, it isn’t necessary to generate nodes considering mesh connectivity condition. In case of the EFGM, computational model can be easily constructed comparing with the FEM. In addition, surface strain distribution can be obtained if the EFGM is employed. Based on the above reasons, the EFGM is applied to compute strain distribution near edge of interface in this study.

Discussion 2 Experiment and image analysis Tensile test and specimen

In this study, tensile test for a specimen, i.e., silicon embedded in resin (Fig.1) is attached to tensile test machine (Fig.2), and surface image of specimen is taken by digital microscope. In this experiment, the tensile load 25N act on specimen and surface image before deformation is measured.

Further, the tensile load 115N act on specimen and surface image after deformation is measured. Load for specimen is measured by load cell, and this value is employed in numerical analysis. Material properties are shown in Table.1.

Fig.1 Schematic diagram of specimen Table.1 Material properties

Young’s modulus [GPa] Poisson ratio

Si 166 0.26

Re 5.49 0.32

(2)

Fig.2 Tensile test machine

Measurement of displacement field based on the DICM In the DICM, the process of computation of displacement vector is divided into rough search and fine search. Rough search is carried out at first stage, and the value of rigid body motion is estimated. After that, deformation of subset is estimated by fine search, and the displacement vector is obtained for each subset. The detail of each search is such as follows.

Rough search in the DICM

In rough search of the DICM, value of rigid body motion at each subset is estimated by cross-correlation coefficient.

Image of rough search is shown in Fig.3. The cross- correlation is calculated by Eq. (1).

!x,y"

!x

^*

,y

^*

"

!u,v"

x y

Before Deformation After Deformation Fig3. Image diagram of rough search

, (1)

where is cross-correlation coefficient, and are number of points for line and column directions at each subset, is light value on specimen surface at points for subset at before deformation, and is average of specimen surface in subset . In addition, is light value on specimen surface at point for subset at after deformation, and is average of specimen surface in subset .

Fine search in the DICM

In the fine search of the DICM, deformation of each subset is estimated by iterative computation. Image of fine search is shown in Fig.4. The problem is to find appropriate subset deformation value i.e., displacement and strain components, such that the highest correlation coefficient is obtained. To solve this problem, the performance function is

defined as Eq. (2). Where, it is necessary to obtain the appropriate subset deformation values i.e., displacement and

strain components and , so as to

minimize the performance function .

, (2)

!x,y"

!x

^*

,y

^*

"

!u,v,w"

!u

!x ! !u

!y

!v

!x ! !v

!y

!w

!x ! !w

!y

x y

Before Deformation After Deformation Fig.4 Image diagram of fine search

where and are light value of specimen surface at before and after deformation. Where, value of

and are defined as Eqs.(3) and (4).

, (3)

, (4) where, and are light value on surface before and after deformation, is fluctuation for light value at before and after deformation, and and are average of light value of specimen surface in subset before and after deformation. In addition and are the point of transferred from arbitrary point of , and those are written as Eqs. (5), (6) and (7).

, (5)

, (6)

. (7)

In Eqs. (5), (6) and (7), and are distance between center point and target point in subset. In addition,

is obtained by interpolation using light value of specimen surface at surrounding pixel (Fig.5). The equation of bi- linear interpolation is written as Eq. (8). [7]

. (8)

Where, and are indicated by written as

Eq.(9), (10), (11) and (12).

(3)

y’

x'

dy

g(i,j) (X,Y)

dx g(x

^*

,y

^*

)

g(i,j+1) g(i+1,j+1)

g(i+1,j)

Fig.5 Image diagram of bi-linear interpolation

, (9)

, (10)

, (11)

. (12)

In Eq. (8), is distance between point and point , and it is necessary to be smaller than 1.0.

Therefore written as Eqs. (13), (14).

(13)

(14)

In this study, steepest descent method is applied to solve this problem. The update equation with respect to unknown variables is shown in Eq. (15).

, (15)

where and are expressed as shown in Eqs. (16) and (17).

, (16)

, (17)

where, is step length, and is number of iterations.

Measurement of strain field based on the MLSM

In the MLSM, the interpolation function for unknown variable i.e. each component of displacement vector, is written as Eq. (18).

, (18) where, and are basis function and unknown coefficient vector and are and

respectively. In addition, is shape function, and , and are computed by Eqs.

(19) - (21).

, (19) , (20)

, (21)

where is number of evaluation nodes in domain of influence (See Fig.6), and , and are weighting function, coordinates at evaluation and reference nodes. In this study, two dimensional first order basis functions are applied to the shape function shown in Eq. (18). In addition, fourth order spline function shown in Eq. (22) is employed as weighting function.

, (22)

where and are distance from to , and radius of domain influence. Consequently, strain components

, are calculated as Eqs.(23) and (24).

, (23)

. (24)

!

X r

_i

r

_mi

:Evaluation node :Referred nodes :Other nodes

!

_a

Fig.6 Domain influence

Discussion 3 Numerical analysis

Computation of strain field based on the EFGM

In this study, numerical analysis for specimen model

(Fig.7) is carried out. Specimen model used by numerical

analysis is the same as specimen used in experiment. In this

analysis, strain singular field near the edge of interface in

specimen model is computed based on the EFGM. Where,

tensile stress that acts on the bottom surface of specimen

model is the same as experiment. Finally, strain is changed

into polar coordinate and influence of cuboidal inclusion is

compared with singular order obtained by eigen analysis.

(4)

x z y

Tensile stress Number of nodes : 17300!

Number of cells : 12899 (a)Specimen (b)1/2 specimen model (c)Model for analysis

Unit[mm]

Re Si

12.6 0.4

6 12 1

x z y Tensile stress x

z y

Si 12

0.63 0.2

1 6 Re

Fig.7 Specimen model Discretization of elastic equation by the EFGM

The equilibrium equation, the strain-displacement relation and the stress-strain relation are written as Eqs. (25), (26), and (27).

, (25)

, (26)

, (27)

where , , and are stress, strain, displacement and elastic coefficient respectively. Where Eqs. (25), (26) and (27) are represented as Eq. (28), (29) and (30). , (28)

, (29)

, (30)

where , , , and are shown in Eq. (31), (32), (33), (34) and (35) respectively. , (31)

, (32)

, (33)

, (34)

, (35)

where and are Lame’s constants, and the constants are written as Eqs. (36), (37). , (36)

, (37)

where, and indicate Young’s modulus and Poisson ratio respectively. Multiplying weighting function for both sides of equilibrium equation and integrating a domain influence (See Fig.6), Eq. (38) is obtained. . (38)

Applying the Green theorem to Eq. (38), Eq. (39) is obtained. , (39)

where is traction force vector, written as Eq.(40). . (40)

Substituting the stress-strain relation to Eq. (39), the Eq. (39) is rewritten as Eq. (41). . (41)

If the weighting function and displacement at an arbitrary point are interpolated by each values at referred nodes in domain of influence based on Galerkin procedure, the interpolation function for each values are written as Eqs. (42) and (43). , (42) , (43)

where is shape function, and is number of referred

nodes in the domain influence. The shape function is

determined by the MLSM, and linear basis and forth order

spline functions are employed as the basis function and

weighting function. Applying the interpolation functions to

the Eq.(41), Eq.(44) is finally obtained.

, (44)

(5)

where the left hand side coefficient matrix and the right hand side vector are the stiffness matrix and external force respectively. In addition is the boundary on domain . The Legendre-Gauss formula is employed as the numerical integration for Eq.(44). In addition, the penalty function method is employed as treatment of essential boundary condition.

Discussion 4 Numerical Examples

A specimen that silicon chip was embedded in resin is employed as the computational model (See Fig.7). In this study, strain distribution near interface edge is computed in enlarged region shown in Fig.8.

Si Re

y x

z z

x Enlargement

Enlargement Region1

Region2

112.5

87.5 Nodal distribution of enlarged region Number of nodes:1073

Minimum cell length:3.125(µm) Interface

Unit[µm]

Fig.8 Computational model

In this study, tensile stress 40MPa is applied to the bottom surface of computational model as a primary load, and this state is treated as before deformation. Furthermore, the tensile stress of 190MPa is applied to the bottom surface of computational model, and this state is treated as after deformation. Displacement distribution near vertex is shown in Fig.9. Distribution of relative displacement with respect to vertex on interface edge is shown in Fig.10. The relative displacement is obtained by , and

indicates displacement at vertex on silicon interface.

y [ ! m ]

x[!m]

Displacement[!m]

Re Si x

y

Interface edge

(a)Region1

y [ ! m ]

x[!m]

Displacement[!m]

Re Si

x

y

Interface edge

(b)Region2

Fig.9 Displacement distribution

y [ ! m ]

x[!m]

Displacement[!m]

Re Si x

y

(a)Region1

y [ ! m ]

x[!m]

Displacement[!m]

Re Si

x

y

(b)Region2

Fig.10 Distribution of relative displacement

In addition, computation of strain field is carried out based on the MLSM using displacement of second, third, fourth, fifth and eighth digit in significant digit. Radius of is 4.7μm. Figs. 11 and 12 shows strain and that computed by the displacement of eighth digit in significant digit.

In Figs.11 and 12, the broken line indicates interface, and

the region surrounded by solid line indicate the region

changed into the polar coordinate system.

(6)

! _x

x[!m]

Re

Si

x z

z [ ! m ]

(a) Region1

! _x

Re Si

x z

z [ ! m ]

x[!m]

(b)Region2

Fig.11 Distribution of strain (In case of eighth digit)

! _z

Re

Si

x z

z [ ! m ]

x[!m]

(a) Region1

! _z

Si Re

x z

z [ ! m ]

x[!m]

(b)Region2

Fig.12 Distribution of strain (In case of eighth digit) Strain and changed into the polar coordinate system and using the Eqs. (45) and (46). Strain and are shown in Figs.13 and 14.

, (45) . (46)

In addition, the distribution of strain for radius direction at on the surface of specimen is shown in Fig.15. It is known that strain for radius direction from singularity point is proportional to , and the parameter indicates order of singularity. In this study, order of singularity is obtained by eigen analysis based on the FEM (See Appendix). [8]

In recent years, this analysis is applied to compute order of singularity at a vertex in transversely isotropic piezoelectric dissimilar material joints, and this is applicable when order of singularity at vertex for bonded structure in three dimension. [9] In this study, the order of singularity at vertex of silicon plate is obtained 0.577. In Fig.15, solid line indicates result of curve fitting for computed strain eighth digit in significant digit. It is seen that the computed strain is good agreement with the result of curve fitting. In addition, it is found that it is necessary to employ the displacement of third digit in significant digit, if computed strain is agreement with the result of curve fitting.

x

z x[!m]

z [ ! m ]

! _r

(a)Region1

(7)

x

z _x[!m]

z[ ! m ]

! _r

(b)Region2

Fig.13 Distribution of strain (In case of eighth digit)

x

z _x[!m]

z [ ! m ]

! _"

"=45°

r O

(a) Region1

! "

"=45°

r O

x z

x[!m]

z [ ! m ]

(b)Region2

Fig.14 Distribution of strain (In case of eighth digit)

Distance from point O[!m]

S tra in "

!

| ! =45° #

:Second digit :Third digit

:Fourth digit

:Fifth digit#

:Eighth digit#

:Curve fittinng using !=0.577

(a)Region1

Distance from point O[!m]

S tra in "

!

| ! =45° #

:Second digit :Third digit

:Fourth digit

:Fifth digit#

:Eighth digit#

:Curve fittinng using !=0.577

(b)Region2

Fig.15 Strain at on surface

Discussion 5 Measurement results of displacement and strain fields

The surface images in region1 and region2 are shown in Figs.16 and 17. These surface images are measured by digital microscope that magnifies 1080 diameters. Light values in surface image are employed for computation of displacement. In Fig.16 and 17, the region surrounded by broken line indicates the region used for computation of displacement and strain.

Re

Si

(b)After loading

S ide of ! fi xe d bounda ry S ide of ! te ns ile l oa di ng

130"m

130 " m

Resolution 150pixel#150pixel

Re

Si

(a)Before loading

Fig.16 Surface image in region 1

(a)Before loading 130!m

130 ! m

Re

Si

Resolution 150pixel"150pixel

Re

Si

(b)After loading

S ide of # fi xe d bounda ry S ide of # te ns ile l oa di ng

Fig.17 Surface image in region 2

The size of target region is set 130×130μm

(150 × 150pixel), and pixel size is approximately 1 × 1 μ m. The

significant digit in target region is third digit by size of target

region and pixel. Therefore, the significant digit in the

(8)

computation of rough search is also third digit. In addition, subset size is 51 ×5 1pixel in displacement computation, and radius of domain influence is given 0.66μm in strain computation.

Distribution of relative displacement with respect to vertex on interface edge obtained by the DICM is shown in Fig.18, broken line denotes interface line. In addition, distribution of displacement on A-A’ line in Fig.18 (a) is shown in Fig.19, and circle mark indicates the displacement component in Fig. 18(a). It is seen that the distribution is step-wise constant, and it can be said that this displacement distribution is not suitable for computation of strain.

Therefore, smoothing technique by the MLSM is applied to the results shown in Fig.18. If the radius of domain influence is given 20μm, results shown in Fig.18 are modified as shown in Fig.20. The distribution of displacement component on A-A’ line is modified as shown in Fig.19 (box mark). Here, strain computation is carried out for the results of Fig. 20. Strain distributions for and are shown in Figs. 21 and 22.

x y

y [ ! m ]

x[!m]

Displacement[!m]

Re Si

A’

A

(a)Region1

y [ ! m ]

x[!m]

Displacement[!m]

Re Si

x y

(b)Region2

Fig.18 Distribution of relative displacement

:Displacement without smoothing!

:Displacement with smoothing

y["m]

D is pl ac em ent v [ " m ]

Fig.19 Gap of displacement on line A-A’

y [ ! m ]

x[!m]

Re Si

x y

Displacement[!m]

(a)Region1

y [ ! m ]

Re Si

x y

x[!m]

Displacement[!m]

(b)Region2

Fig.20 Distribution of relative displacement with

smoothing

(9)

! _x

y [ ! m ]

x[!m]

Re Si

x y

(a)Region1

! _x

y [ ! m ]

Re Si

x y

x[!m]

(b)Region2

Fig.21 Distribution of strain

! _y

y [ ! m ]

x[!m]

Si Re

x y

(a)Region1

! _y

y [ ! m ] Re Si

x y

x[!m]

(b)Region2

Fig.22 Distribution of strain

From these results, it is found that there is area of high strain value, and it appears that crack occurs near interface.

Therefore, measuring surface height variation on line B-B’

shown in Fig. 23, it could be confirmed that there is crack near interface between silicon plate and resin.

x

y ^y[!m]

S urfa ce he ight [ ! m ]

Crack

Interface

Re Si

Re

Si B’

B

Fig.23 Surface image and variation of surface height on line B-B’

Conclusions

In this study, measurement of the displacement and strain fields by image analysis based on the DICM and the MLSM and numerical analysis based on the EFGM were carried out.

Conclusions in this study are written as follows.

l It was found that appropriate strain distribution can be obtained in case that displacement value of third digit in significant digit is employed.

l In this study, displacement distribution obtained by the DICM was step-wise constant. If this kind of displacement distribution is obtained, strain distribution can’t be obtained. Therefore, employing smoothing technique by the MLSM to the displacement distribution by the DICM, appropriate strain distribution could be computed in this study.

l It appears that it is necessary to apply light value of

higher resolution in target region, because size of pixel

in the DICM is larger than value of relative

displacement obtained by the EFGM.

(10)

Appendix

Computation of order of singularity based on the FEM The derivation of characteristic equation shown in the reference of Pegeau et. al. [8] is simply introduced. In this formulation, the computational region is defined by spherical configuration whose radius is , and the spherical coordinate system is introduced (See Fig. 24) . As the final form of the derived equation, the equation on the spherical surface is obtained. Therefore, the surface domain is divided into finite elements, and the computation for the characteristic equation is carried out. In the reference of Pageau et. al. [8], the quadratic isoparametric element is employed as the finite elements. If displacements for each element are expressed by interpolation function shown in Eq.

(47) and the interpolation function is substituted to equation of the principle of virtual work, a characteristic equation is finally derived as shown in Eq. (48).

, (47)

, (48)

where is expressed by , and and represent spherical surface. In addition, indicates the shape function.

In the Eq. (48), indicates characteristic root and vector denotes superposed displacement vector in entire domain, and matrices , and represent the coefficient matrices derived by finite element procedure. Detail of this formulation is shown in reference. [8]. The characteristic root is obtained by solving the Eq. (48) based on eigen analysis. Relation ship between the characteristic root and order of singularity is expressed by . If the parameter is , it is denoted that stress fields has a stress singularity. On the other hand, if the parameter is

, it is denoted that the stress singularity disappears.

y x

z

r ^Q r ₀

!

O

"

#

#=-1

#=1

$

%

$=1

$=-1

%=1

%=-1

1 2 3

4 5 7 6

8 Fig.24 Computational model and isoparametric quadratic element

In this study, characteristic root and order of singularity was obtained as shown in Table. 2. Applying order of singularity to computed strain distribution, the value 0.577 is adopted as order of singularity in this study (See Fig. 15).

Table.2 Computed characteristic root and order of singularity

1 0.367 0.633

2 0.423 0.577

Acknowledgments

This work was supported by "Program for High Reliable Materials Designand Manufacturing in Nagaoka University of Technology", and Grants-in-Aid for scientific Research(B) Grant Number 90143693. We wish to thank you for Mr. Hiroyuki TANAKA at SUMITOMO BAKELITE CO., LTD for providing us specimen employed in this study.

References

1. Koguchi. H, Hoshi. K, Evaluation of joining strength of silicon-resin interface at a vertex in 3D joint structures, InterPACK2011, (2011), pp.1-8.

2. Pan. B, Asundi. A, Xie. H and Gao. J, Digital image correlation using iterative least squares and pointwise least squares for displacement field and strain field measurements, Optics and Lasers in Engineering, vol.47, (2009), pp.865-874.

3. Kirugulige. M.S, Tippur. H.V and Denny. T.S, Measurement of transient deformations using digital image correlation method and high-speed photography:

application to dynamic fracture, Applied Optics, vol.46, (2007), pp. 5083-5096.

4. Kurahashi. T, Yamada. K and Koguchi. H, Measurement of Displacement and Strain Fields near Interface of Bonded Structure Based on Digital Image Correlation Method Using Surface Configuration Data, International Conference on Advanced Technology in Experimental Mechanics (2011), pp.1-10.

5. Kurahashi. T, Ishikawa. A and Koguchi. H, Evaluation of Intensity of Stress Singularity for 3D Dissimilar Material Joints Based on Mesh Free Method,11th International Conference on the Mechanical Behavior of Materials, Procedia Engineering, vol.10, (2011), pp.3095-3100.

6. Belytschko. T, Lu. Y. Y and Gu. L, Element-free Garerkin methods, Intternational. Journal for numerical Methods in Engineering, vol.37, (1994), pp.229-256.

7. Bruck. H.A, McNeill. S.R, Sutton. M.A and Peters. W.H, Digital image correlation using Newton-Raphson method of paritual differential correction, Experimental Mechanics, vol.29, (1989), pp.261-267.

8. Pageau. S.S and Biggers Jr. S.B, Finite element evaluation of free - edge singular stress fields in anisotropic materials, Intternational. Journal for numerical Methods in Engineering, vol.38, (1995), pp.2225-2239.