**Finite Element Analysis of Loading Area Effect on Sandwich Panel Behaviour Beyond the Yield Limit**

Salih Akour and Hussein Maaitah

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/50433

## **1. Introduction**

352 Finite Element Analysis – New Trends and Developments

2011; 67(1):14-24.

28(5),769-776.

Structures, 1999; 72:109-126.

Structural Engineering; 1987.

London, August; 1993.

Symposium on Tubular Structures, Madrid, Spain.

COPPE/UFRJ, Rio de Janeiro, Brazil, pp. 1-309; 2004.

[6] Kuhlmann U, H-P Günther, Saul R, Häderle, M.-U. 2003. Welded circular hollow section (CHS) joints in bridges. ISTS 2003: Proceedings of the 10th International

[7] Leitão FN, Silva JGS da, Vellasco PCG da S, Lima LRO de, Andrade SAL de. Composite (steel-concrete) highway fatigue assessment. Journal of Constructional Steel Research,

[8] Varela WD. Modelo teórico-experimental para análises de vibrações induzidas por pessoas caminhando sobre lajes de edifícios (Theoretical-experimental model to analyse vibrations induced by people walking on floor slabs of buildings), PhD Thesis (in Portuguese), Federal University of Rio de Janeiro, Civil Engineering Department,

[9] Murray TM, Allen DE, Ungar EE. Floor vibrations due to human activity, Steel Design Guide Series, American Institute of Steel Construction, AISC, Chicago, USA; 2003. [10] Pimentel RL, Pavic A, Waldron, P. Evaluation of design requirements for footbridges excited by vertical forces from walking. Canadian Journal of Civil Engineering, 2001;

[11] Chen, Y. Finite element analysis for walking vibration problems for composite precast building floors using ADINA: modelling, simulation and comparison. Computer &

[12] Bachmann, H, Ammann, W. Vibrations in structures induced by man and machines, Structural Engineering Document 3e. International Association for Bridges and

[13] ANSYS. Swanson Analysis Systems, Inc. P.O. Box 65, Johnson Road, Houston, PA,

[14] International Standard Organization / ISO 2631-2. Evaluation of human exposure to whole-body vibration, Part 2: Human Exposure to Continuous and Shock-Induced

[15] Comité Euro-international du Béton. CEB-FIP: Bulletin d'information, N0 209, England,

15342-0065, Version 10.0, Basic analysis procedures, Second edition; 2003.

Vibrations in Buildings (1 to 80Hz), International Standard; 1989.

Research efforts continuously are looking for new, better and efficient construction materials. The main goal of these researches is to improve the structural efficiency, performance and durability. New materials typically bring new challenges to designer who utilizes these new materials. In the past decades various sandwich panels have been implemented in aerospace, marine, architectural and transportation industry. Light-weight, excellent corrosion characteristics and rapid installation capabilities created tremendous opportunities for these sandwich panels in industry. Sandwich panel normally consists of a low-density core material sandwiched between two high modulus face skins to produce a lightweight panel with exceptional stiffness as shown in Figure 1. Face skins act like flanges of an I-beam. These faces are typically bonded to a core to achieve the composite action and to transfer the forces between sandwich panel components.

## **1.1. Main principles of sandwich structures**

Typical sandwich composite construction consists of three main components as illustrated in Figure 1. The sandwich consists of two thin, stiff and strong faces are separated by thick, light and weaker core. Faces and core materials are bonded together with an adhesive to facilitate the load transfer mechanism between the components, therefore effectively utilize all the materials used. The two faces are placed at a distance from each other to increase the moment of inertia, and consequently the flexural rigidity, about the neutral axis of the structure.

In sandwich structure, typically the core material is not rigid compared to face sheets; therefore, the shear deflection within the core is insignificant in most cases. The shear deflection in the faces can be also neglected. The effect of shear rigidity in the core is shown

in Figure 2. Figure 2a shows an ideal sandwich beam using relatively stiff core, therefore the two faces cooperate without sliding relative to each other. Figure 2b shows a sandwich beam using weak core, therefore the faces are no longer coupled together effectively and each face works independently as plates in bending. The use of weak core in shear results in significant loss of the efficiency of the sandwich structures. In a typical sandwich panel the faces carry the tensile and compressive stresses. The local flexural rigidity of each face is typically small and can be ignored. Materials such as steel, stainless steel, aluminum and fiber reinforced polymer materials are often used as materials for the face. The core has several important functions. It has to be stiff enough to maintain the distance between the two faces constant. It should be also rigid to resist the shear forces and to prevent sliding the faces relative to each other. Rigidity of the core forces the two faces to cooperate with each other in composite action. If these conditions are not fulfilled, the faces behave as two independent beams or panels, and the sandwich effect will be totally lost. Furthermore, rigidity of the core should be sufficient to maintain the faces nearly flat, therefore prevent possibility of buckling of the faces under the influence of compressive stress in their plane. The adhesive between the faces and the core must be able to transfer the shear forces between the face and the core.

Finite Element Analysis of Loading Area Effect on Sandwich Panel Behaviour Beyond the Yield Limit 355

low water absorption. Applications include the construction of bulkheads, hulls, decks,

*Transportation Industry:* High strength-to-weight ratios of sandwich composites offer great advantages to the transportation industry. The insulating, sound damping properties and low cost properties make them the choice materials for the constructions of walls, floors,

*Architectural Industry:* The foam offers an excellent thermal and acoustical insulation which makes it ideal choice for the architectural industry. Typical applications include structural

Work on the theoretical description of sandwich structure behaviour began after World War Two. (Plantema, 1966) published the first book about sandwich structures, followed by books by (Allen, 1969), and more recently by (Zenkert, 1995). Although (Triantafillou and Gibson, 1987) developed a method to design for minimum weight, and reported the failure mode map of sandwich construction, without considering the post yield state of the sandwich structure.

The basic sandwich structure theory presented in all these texts is generally called the

columns, portable buildings, office partitions, countertops and building facades.

**Figure 2.** Presentation of the effect of a) rigid and b) week core.

doors, panels and roofs for vans, trucks, trailers and trains.

classical sandwich theory. This theory assumes that:

The face sheets carry the entire bending load.

1. The core and face sheets are elastic.

2. The overall length to thickness ratio is high.

additional load carrying capacity once the core yields.

Core compression is negligible.

The core carries the entire shear load in sandwich beams and plates.

This theory states that the above–mentioned assumptions are true if:

3. The face sheet thickness is small compared to the overall thickness.

4. The ratio of mechanical properties between the face sheet and the core is high. With these assumptions, a sandwich structure is considered to be incapable of acquiring

(Mercado and Sikarskie, 2000) reported that the load carried by sandwich structures continue to increase after core yielding. Knowing that the core could not carry additional load after yield, this increasing load carrying capacity of post yield sandwich structure

transoms and furniture.

**1.3. Literature review** 

**Figure 1.** Schematic of sandwich construction

### **1.2. Applications**

Sandwich construction provides efficient utilization of the materials used for each component to its ultimate limit (Zenkert, 1997). The sandwich structure offers also a very high stiffness-to-weight ratio. It enhances structure flexural rigidity without adding substantial weight and makes it more advantageous as compared to composite materials. Sandwich constructions have superior fatigue strength and exhibit superior acoustical and thermal insulation. Sandwich composites could be used in a wide variety of applications such as:

*Aerospace Industry:* Sandwich composites are increasingly being used in the aerospace industry because of their bending stiffness-to-weight ratio. Floorboards, composite wing, horizontal stabilizer, composite rudder, landing gear door, speed brake, flap segments, aircraft interior and wingspans are typically made of sandwich composites.

*Marine Industry:* Sandwich composites are ideally suited for the marine industries most advanced designs. The foam cores meet the critical requirements of strength, buoyancy and low water absorption. Applications include the construction of bulkheads, hulls, decks, transoms and furniture.

**Figure 2.** Presentation of the effect of a) rigid and b) week core.

*Transportation Industry:* High strength-to-weight ratios of sandwich composites offer great advantages to the transportation industry. The insulating, sound damping properties and low cost properties make them the choice materials for the constructions of walls, floors, doors, panels and roofs for vans, trucks, trailers and trains.

*Architectural Industry:* The foam offers an excellent thermal and acoustical insulation which makes it ideal choice for the architectural industry. Typical applications include structural columns, portable buildings, office partitions, countertops and building facades.

## **1.3. Literature review**

354 Finite Element Analysis – New Trends and Developments

between the face and the core.

**Figure 1.** Schematic of sandwich construction

**1.2. Applications** 

such as:

in Figure 2. Figure 2a shows an ideal sandwich beam using relatively stiff core, therefore the two faces cooperate without sliding relative to each other. Figure 2b shows a sandwich beam using weak core, therefore the faces are no longer coupled together effectively and each face works independently as plates in bending. The use of weak core in shear results in significant loss of the efficiency of the sandwich structures. In a typical sandwich panel the faces carry the tensile and compressive stresses. The local flexural rigidity of each face is typically small and can be ignored. Materials such as steel, stainless steel, aluminum and fiber reinforced polymer materials are often used as materials for the face. The core has several important functions. It has to be stiff enough to maintain the distance between the two faces constant. It should be also rigid to resist the shear forces and to prevent sliding the faces relative to each other. Rigidity of the core forces the two faces to cooperate with each other in composite action. If these conditions are not fulfilled, the faces behave as two independent beams or panels, and the sandwich effect will be totally lost. Furthermore, rigidity of the core should be sufficient to maintain the faces nearly flat, therefore prevent possibility of buckling of the faces under the influence of compressive stress in their plane. The adhesive between the faces and the core must be able to transfer the shear forces

Sandwich construction provides efficient utilization of the materials used for each component to its ultimate limit (Zenkert, 1997). The sandwich structure offers also a very high stiffness-to-weight ratio. It enhances structure flexural rigidity without adding substantial weight and makes it more advantageous as compared to composite materials. Sandwich constructions have superior fatigue strength and exhibit superior acoustical and thermal insulation. Sandwich composites could be used in a wide variety of applications

*Aerospace Industry:* Sandwich composites are increasingly being used in the aerospace industry because of their bending stiffness-to-weight ratio. Floorboards, composite wing, horizontal stabilizer, composite rudder, landing gear door, speed brake, flap segments,

*Marine Industry:* Sandwich composites are ideally suited for the marine industries most advanced designs. The foam cores meet the critical requirements of strength, buoyancy and

aircraft interior and wingspans are typically made of sandwich composites.

Work on the theoretical description of sandwich structure behaviour began after World War Two. (Plantema, 1966) published the first book about sandwich structures, followed by books by (Allen, 1969), and more recently by (Zenkert, 1995). Although (Triantafillou and Gibson, 1987) developed a method to design for minimum weight, and reported the failure mode map of sandwich construction, without considering the post yield state of the sandwich structure.

The basic sandwich structure theory presented in all these texts is generally called the classical sandwich theory. This theory assumes that:

	- 1. The core and face sheets are elastic.
	- 2. The overall length to thickness ratio is high.
	- 3. The face sheet thickness is small compared to the overall thickness.
	- 4. The ratio of mechanical properties between the face sheet and the core is high.

With these assumptions, a sandwich structure is considered to be incapable of acquiring additional load carrying capacity once the core yields.

(Mercado and Sikarskie, 2000) reported that the load carried by sandwich structures continue to increase after core yielding. Knowing that the core could not carry additional load after yield, this increasing load carrying capacity of post yield sandwich structure

initiates the postulation that the additional shear load was transferred to the face sheets. To account for the above-mentioned phenomenon, (Mercado et al, 1999) developed a higher order theory by including a bilinear core material module. This theory yields a fairly accurate prediction on the deflection of a foam cored sandwich structure in four point bending (Mercado et al, 2000). In addition, this theory does not take into account the core compression under localized load, or any geometric non-linearity. The classical sandwich beam theory also assumes that in-plane displacements of the core through its depth are linear. In other words, it was assumed that the core thickness remains constant and crosssections perpendicular to the neutral axis remain plane after deformation. This assumption is generally true for traditional core material such as metallic honeycomb. However, this assumption is not suitable for soft, foam-based cores, especially when the sandwich structure is subjected to a concentrated load (Thomsen, 1995). With a much lower rigidity compared to metallic honeycomb, foam-based cored sandwich structures are susceptible to localized failure. Insufficient support to the face sheets due to core compression near the application points of concentrated loads can lead to failures such as face sheet/ core delamination, face sheet buckling, and face sheet yielding. This localized non-linearity is reported by many researchers such as (Thomsen, 1995), (Thomsen, 1993), (Rothschild 1994), (Caprino, 2000), and (Gdoutos et al, 2001). The shear distribution at localized failure points has not been well defined. (Miers, 2001) investigated the effect of localized strengthening inserts on the overall stiffness of a sandwich structure. This localized strengthening increases the rigidity of the sandwich structure, but the addition of high stiffness inserts complicates the manufacturing process of sandwich structure.

Finite Element Analysis of Loading Area Effect on Sandwich Panel Behaviour Beyond the Yield Limit 357

3. The effect the size of the distributed load area on the behavior of the sandwich panel beyond the core yield limit for different types of materials is investigated. These parameters are the determining factors of the significance of geometric non-linearity

1. Localized core yielding occurs mainly through core compression. Therefore, analysis should be done using material properties determined from compression test. 2. For practical purposes, the assumptions that have been made in developing the

3. The Finite Element Model (FEM) is extended to include the relative dominance of core

4. Localized loads are modeled as load on small partitioned area to better simulate the

This section presents the physical model of the sandwich panel, which includes geometry,

The sandwich panel consists of two face sheets made of metal. The thickness of each face is **t**. Soft core of **c** thickness is sandwiched between those face sheets. The core material is made of foam which is soft compared to the face sheets .The panel is square in shape. The side length is designated by **a.** Figure 3 illustrates the sandwich panel geometry whereas the

1. Post yield behavior of sandwich panel.

and core material nonlinearity

shear failure and face sheet yielding.

actual loading condition.

**2.1. Sandwich panel geometry** 

**2. Physical model** 

2. Effect of geometric non-linearity under distributed loads.

The above investigation is done in view of the following points:

5. Experimental verification is conducted for selected cases.

dimensions of the sandwich panel are shown in Table 1.

**Figure 3.** Illustration sandwich plate geometry.

sandwich panel theory eliminated part of the problem physics.

boundary conditions as well as the materials used in the investigation.

To design an efficient sandwich structure, it is vital to understand the behavior of each layer in the structure. Classical sandwich theory (Zenkert 1995, Plantema 1966, Allen 1969), higher order theory by Mercado (2000) and high order theory developed by Frostig et al. (1992) could predict the sandwich panel behavior fairly accurate in the linear range. However, these theories could not give an accurate prediction of the sandwich structure behavior after core yielding. Large deflection of sandwich structures due to core yielding could vary the direction of the applied load on the structure.

### **1.4. Research objective**

To design an efficient sandwich structure, it is vital to understand the load distribution pattern in each layer of the structure. Most of the previous efforts are made by using classical sandwich theory, and higher order theory, where high order theory predicted the sandwich panel behavior fairly well in the linear range. However, these theories could not give an accurate prediction of the shear distribution in each layer after core yielding. Large deflection of sandwich structures due to core yielding could vary the direction of the applied load on the structure. Change in loading direction would obviously change the shear distribution in the sandwich structure. In order to investigate the exact change of shear distribution due to distributed loads, as well as geometric nonlinearity and localized core failure, finite element analysis is used in this research effort. The main objective of this research is to investigate the following:

1. Post yield behavior of sandwich panel.

356 Finite Element Analysis – New Trends and Developments

complicates the manufacturing process of sandwich structure.

direction of the applied load on the structure.

research is to investigate the following:

**1.4. Research objective** 

initiates the postulation that the additional shear load was transferred to the face sheets. To account for the above-mentioned phenomenon, (Mercado et al, 1999) developed a higher order theory by including a bilinear core material module. This theory yields a fairly accurate prediction on the deflection of a foam cored sandwich structure in four point bending (Mercado et al, 2000). In addition, this theory does not take into account the core compression under localized load, or any geometric non-linearity. The classical sandwich beam theory also assumes that in-plane displacements of the core through its depth are linear. In other words, it was assumed that the core thickness remains constant and crosssections perpendicular to the neutral axis remain plane after deformation. This assumption is generally true for traditional core material such as metallic honeycomb. However, this assumption is not suitable for soft, foam-based cores, especially when the sandwich structure is subjected to a concentrated load (Thomsen, 1995). With a much lower rigidity compared to metallic honeycomb, foam-based cored sandwich structures are susceptible to localized failure. Insufficient support to the face sheets due to core compression near the application points of concentrated loads can lead to failures such as face sheet/ core delamination, face sheet buckling, and face sheet yielding. This localized non-linearity is reported by many researchers such as (Thomsen, 1995), (Thomsen, 1993), (Rothschild 1994), (Caprino, 2000), and (Gdoutos et al, 2001). The shear distribution at localized failure points has not been well defined. (Miers, 2001) investigated the effect of localized strengthening inserts on the overall stiffness of a sandwich structure. This localized strengthening increases the rigidity of the sandwich structure, but the addition of high stiffness inserts

To design an efficient sandwich structure, it is vital to understand the behavior of each layer in the structure. Classical sandwich theory (Zenkert 1995, Plantema 1966, Allen 1969), higher order theory by Mercado (2000) and high order theory developed by Frostig et al. (1992) could predict the sandwich panel behavior fairly accurate in the linear range. However, these theories could not give an accurate prediction of the sandwich structure behavior after core yielding. Large deflection of sandwich structures due to core yielding could vary the

To design an efficient sandwich structure, it is vital to understand the load distribution pattern in each layer of the structure. Most of the previous efforts are made by using classical sandwich theory, and higher order theory, where high order theory predicted the sandwich panel behavior fairly well in the linear range. However, these theories could not give an accurate prediction of the shear distribution in each layer after core yielding. Large deflection of sandwich structures due to core yielding could vary the direction of the applied load on the structure. Change in loading direction would obviously change the shear distribution in the sandwich structure. In order to investigate the exact change of shear distribution due to distributed loads, as well as geometric nonlinearity and localized core failure, finite element analysis is used in this research effort. The main objective of this


The above investigation is done in view of the following points:


## **2. Physical model**

This section presents the physical model of the sandwich panel, which includes geometry, boundary conditions as well as the materials used in the investigation.

## **2.1. Sandwich panel geometry**

The sandwich panel consists of two face sheets made of metal. The thickness of each face is **t**. Soft core of **c** thickness is sandwiched between those face sheets. The core material is made of foam which is soft compared to the face sheets .The panel is square in shape. The side length is designated by **a.** Figure 3 illustrates the sandwich panel geometry whereas the dimensions of the sandwich panel are shown in Table 1.

**Figure 3.** Illustration sandwich plate geometry.


Finite Element Analysis of Loading Area Effect on Sandwich Panel Behaviour Beyond the Yield Limit 359

reaches fracture limit. The distributed load is applied on the top surface of the sandwich panel. The area on which the distributed load is applied (see Figure 5 and 7) is located at the middle of the top face sheet plate. The loading area at the middle top face of sandwich panel is square in shape. This area has been varied from 100X100 mm2 through 200X200 mm2, 400X400 mm2, and 600X600 mm2 so the ratio of these areas relative to the total area of the

sandwich panel is 1/36, 4/36, 16/36 and 36/36 respectively.

**Figure 4.** Sandwich panel boundary condition for a) X-Y plane and b) Y-Z plane.

thickness is selected to be 30mm as shown in Table 1.

In the current research, different materials are used. Their modulus of elasticity is varying from 37.5 MPa through 138.6 MPa, 180 MPa, and 402.6 MPa as shown in Table 2. Core

(b)

(a)

The core of sandwich structure is used to separate the two faces, most often identical in material and thickness, which primarily resist the in plane and bending load. The core is mainly subjected to shear so that the core shear strain produces global deformations and core shear stresses. Thus, core must be chosen such that not to fail under applied transverse load. It should have shear modulus that is high enough to give the required stiffness.

*2.4.2. Core material* 

**2.5. Material properties** 

**Table 1.** The value of the parameters shown in Figure 3.

## **2.2. Assumptions**

This research takes into consideration the geometric non-linearity as well as the material nonlinearity. The following assumptions are made to simplify the model without losing the problem physics:


Due to the significantly higher yield strength and modulus of elasticity of the face sheets compared to the core, face sheets are assumed to remain elastic throughout the loading for simply supported panel. The analysis stops when the face sheets start to yield.

3. Geometric non-linearity has a significant effect: Geometric non-linearity is considered to have significant effect on the load distribution on each layer of the sandwich structure.

## **2.3. Boundary condition**

Due to the symmetry of the sandwich panel (symmetric over X-axis and symmetric over Zaxis), only quarter of it is being modeled. Such symmetric boundary conditions are applied of the X-axis and Z-axis. The two planes of symmetry of the panel have symmetric boundary conditions, (see Fig. 4). A simply supported boundary condition is applied to strip area of the quarter panel as shown in Fig. 5. This simulates the simply supported condition of the panel. The loading area is square in shape, its side length varies in steps of 100, 200, 400 and 600mm for full panel dimension. But when dealing with quarter panel, the side length is 50, 100, 200, and 300mm

## **2.4. Study parameters**

The main parameters that have influence on the performance of the sandwich plate are, the loading area on which the load is distributed and the core material stiffness.

## *2.4.1. Loading*

The load is applied to the sandwich top face sheet as a distributed load which is increased gradually (step by step) till the face sheet stress reaches yield stress or the core material reaches fracture limit. The distributed load is applied on the top surface of the sandwich panel. The area on which the distributed load is applied (see Figure 5 and 7) is located at the middle of the top face sheet plate. The loading area at the middle top face of sandwich panel is square in shape. This area has been varied from 100X100 mm2 through 200X200 mm2, 400X400 mm2, and 600X600 mm2 so the ratio of these areas relative to the total area of the sandwich panel is 1/36, 4/36, 16/36 and 36/36 respectively.

**Figure 4.** Sandwich panel boundary condition for a) X-Y plane and b) Y-Z plane.

### *2.4.2. Core material*

358 Finite Element Analysis – New Trends and Developments

**Table 1.** The value of the parameters shown in Figure 3.

1. Face sheets and core are perfectly bonded.

3. Geometric non-linearity has a significant effect:

on each layer of the sandwich structure.

2. Face sheets remain elastic at all time.

**2.2. Assumptions** 

problem physics:

yield.

**2.3. Boundary condition** 

100, 200, and 300mm

*2.4.1. Loading* 

**2.4. Study parameters** 

Parameter Dimension Note

a 600 mm Constant t 1.0 mm Constant c 30 mm Constant

This research takes into consideration the geometric non-linearity as well as the material nonlinearity. The following assumptions are made to simplify the model without losing the

Due to the significantly higher yield strength and modulus of elasticity of the face sheets compared to the core, face sheets are assumed to remain elastic throughout the loading for simply supported panel. The analysis stops when the face sheets start to

Geometric non-linearity is considered to have significant effect on the load distribution

Due to the symmetry of the sandwich panel (symmetric over X-axis and symmetric over Zaxis), only quarter of it is being modeled. Such symmetric boundary conditions are applied of the X-axis and Z-axis. The two planes of symmetry of the panel have symmetric boundary conditions, (see Fig. 4). A simply supported boundary condition is applied to strip area of the quarter panel as shown in Fig. 5. This simulates the simply supported condition of the panel. The loading area is square in shape, its side length varies in steps of 100, 200, 400 and 600mm for full panel dimension. But when dealing with quarter panel, the side length is 50,

The main parameters that have influence on the performance of the sandwich plate are, the

The load is applied to the sandwich top face sheet as a distributed load which is increased gradually (step by step) till the face sheet stress reaches yield stress or the core material

loading area on which the load is distributed and the core material stiffness.

The FEM model assumes no delamination occur between layers.

In the current research, different materials are used. Their modulus of elasticity is varying from 37.5 MPa through 138.6 MPa, 180 MPa, and 402.6 MPa as shown in Table 2. Core thickness is selected to be 30mm as shown in Table 1.

## **2.5. Material properties**

The core of sandwich structure is used to separate the two faces, most often identical in material and thickness, which primarily resist the in plane and bending load. The core is mainly subjected to shear so that the core shear strain produces global deformations and core shear stresses. Thus, core must be chosen such that not to fail under applied transverse load. It should have shear modulus that is high enough to give the required stiffness. Furthermore, its young's modulus normal to the faces should be high enough to prevent contraction of the core thickness and therefore a rapid decrease in flexural rigidity. The core should have low density in order to add as little as possible to the total weight of sandwich structure. Because of low density requirement, core materials are very different from face sheet materials. A detailed characterization of their mechanical behavior is essential for their efficient use in structural application. Four types of foam H100, H250, AirexR63.50 and Herex C70.200 are investigated.

Finite Element Analysis of Loading Area Effect on Sandwich Panel Behaviour Beyond the Yield Limit 361

Airex R63.50 has high fatigue strength, high three-dimensional formability, and high resistance to dynamic loads. Materials in Airex R63 family are widely used in the production of marine hulls and lightweight cars due to the appreciation of their low density

Material properties of the HerexC70.200 foam core is obtained from (Rao, 2002) work. Herex C70.200 is an isotropic and stiff foam material with high stiffness and strength to weight ratios. The materials in Herex C70 family have excellent chemical resistance and low thermal conductivity and water absorption. The appreciation of these inherent properties of Herex C70 materials makes this material a popular choice for the core materials of structural sandwich structures in marine and railway applications. The stress strain curve of this

In this research a first-order idealized core material property module suggested by (Mercado, Sikarskie, 1999) is used. This first-order idealized model, also called the bi-linear model, describes the material properties of the core with the stress strain curve as shown on

The other material used in this research is linked PVC close called cellular foam (divinycell). The type of divinycell, H100, H250 with densities of 100 and 250 kg/m3, their mechanical properties are stated in Table 2 and their stress strain curves are shown in Figure 6b and 6d respectively.

**Material Property source Young's modulus (MPa) Poisson's ratio Shear modulus (Mpa) Shear strength (Mpa) 0.2% offset yield strength (Mpa) Strain at yield popup (mm/mm)** 

Core A : AirexR63.50 Rao, 2002 37.5 0.335 14.05 0.45 0.637 0.019

C70.200 Rao, 2002 180 0.37 65.69 1.6 2.554 0.0162

Core D: H250 Kuang, 2001 402.6 0.35 117.2 4.5 5 0.014

This section presents the development of finite element models for simply supported sandwich panel. Detailed descriptions of the boundary conditions, element types, and the

Core B: H100 Kuang, 2001 138.6 0.35 47.574 1.2 1.5 0.0108225

Gall 1991 69,000 0.33 25,000 120 145 Not

available

and high strength and stiffness to weight ratio. Airex R63.50 is presented in Table 2.

material is presented in Figure 6.

Figure 6a and 6c.

Face sheet :

Core C: Herex

Aluminum 3003-H14

**3. Finite element model** 

Boyer and

**Table 2.** Compression of sandwich panel material properties

**Figure 5.** Panel span overview of quarter sandwich panel for different loading area

## *2.5.1. Mechanical properties for face sheet*

Material properties for the sandwich plate face sheets are taken from (Boyer and Gall (Eds.), 1991). Aluminum 3003-H14 is a type of aluminum alloy that has high resistance to corrosion and is easy to weld is used in this investigation. The 3003-aluminum family is normally used in the production of cooking utensils, chemical equipment, and pressure vessels. The face sheets are assumed to remain elastic at all times. Therefore only elastic material properties are required for the face sheets and they are presented in Table 2.

## *2.5.2. Mechanical properties for core*

This subsection presents the core material properties used to model the sandwich panel. In all cases, face sheets of the sandwich structures are assumed to remain elastic throughout the analyses. Therefore, only core materials require a good post yield behavior descriptions. The core materials undergo plastic deformation; hence there is a need to obtain a full description of the core materials' behavior upon yield initiation.

Airex R63.50 has high fatigue strength, high three-dimensional formability, and high resistance to dynamic loads. Materials in Airex R63 family are widely used in the production of marine hulls and lightweight cars due to the appreciation of their low density and high strength and stiffness to weight ratio. Airex R63.50 is presented in Table 2.

Material properties of the HerexC70.200 foam core is obtained from (Rao, 2002) work. Herex C70.200 is an isotropic and stiff foam material with high stiffness and strength to weight ratios. The materials in Herex C70 family have excellent chemical resistance and low thermal conductivity and water absorption. The appreciation of these inherent properties of Herex C70 materials makes this material a popular choice for the core materials of structural sandwich structures in marine and railway applications. The stress strain curve of this material is presented in Figure 6.

In this research a first-order idealized core material property module suggested by (Mercado, Sikarskie, 1999) is used. This first-order idealized model, also called the bi-linear model, describes the material properties of the core with the stress strain curve as shown on Figure 6a and 6c.

The other material used in this research is linked PVC close called cellular foam (divinycell). The type of divinycell, H100, H250 with densities of 100 and 250 kg/m3, their mechanical properties are stated in Table 2 and their stress strain curves are shown in Figure 6b and 6d respectively.


**Table 2.** Compression of sandwich panel material properties

## **3. Finite element model**

360 Finite Element Analysis – New Trends and Developments

Herex C70.200 are investigated.

Furthermore, its young's modulus normal to the faces should be high enough to prevent contraction of the core thickness and therefore a rapid decrease in flexural rigidity. The core should have low density in order to add as little as possible to the total weight of sandwich structure. Because of low density requirement, core materials are very different from face sheet materials. A detailed characterization of their mechanical behavior is essential for their efficient use in structural application. Four types of foam H100, H250, AirexR63.50 and

**Figure 5.** Panel span overview of quarter sandwich panel for different loading area

are required for the face sheets and they are presented in Table 2.

description of the core materials' behavior upon yield initiation.

Material properties for the sandwich plate face sheets are taken from (Boyer and Gall (Eds.), 1991). Aluminum 3003-H14 is a type of aluminum alloy that has high resistance to corrosion and is easy to weld is used in this investigation. The 3003-aluminum family is normally used in the production of cooking utensils, chemical equipment, and pressure vessels. The face sheets are assumed to remain elastic at all times. Therefore only elastic material properties

This subsection presents the core material properties used to model the sandwich panel. In all cases, face sheets of the sandwich structures are assumed to remain elastic throughout the analyses. Therefore, only core materials require a good post yield behavior descriptions. The core materials undergo plastic deformation; hence there is a need to obtain a full

*2.5.1. Mechanical properties for face sheet* 

*2.5.2. Mechanical properties for core* 

This section presents the development of finite element models for simply supported sandwich panel. Detailed descriptions of the boundary conditions, element types, and the

loading are presented in the coming subsections. The finite element software used in the development of the finite models is (I-DEAS Master Series 10 1999). The relatively robust and user-friendly solid modeling and finite element meshing interface are the main advantages of this solid modeling and finite element software.

Finite Element Analysis of Loading Area Effect on Sandwich Panel Behaviour Beyond the Yield Limit 363

Due to the high yield strength and high modulus of elasticity of the sandwich face sheets compared to the core, face sheets are assumed to remain elastic throughout the

The loading cases considered are modeled quasi-static instead of dynamic. Incremental loadings are applied slowly during the actual experiments (i.e. simulates exactly the real situation). Therefore, the type of analysis done for this research effort is "static,

Geometric non-linearity is considered to have significant effect on the load distribution on each layer of the sandwich structure. Therefore, all finite element analysis that is done takes into consideration the geometric non-linearity. This is the main difference between the numerical models and the theoretical models. Classical sandwich plate theory and higher order theory do not take shape change of the sandwich structures

5. The panel is simply supported from all sides. It is partitioned into three layers, forming

The symmetric nature of the problem allows only quarter of the whole panel to be meshed. The boundary conditions applied are shown on Figures 4 and 5. The two planes of symmetry of the panel have symmetric boundary conditions, where in-plane displacements and rotation about an axis respective normal to the symmetry plane is allowed. A simply supported boundary condition is applied to the two other sides of the quarter panel. A distributed load is applied on the top surface of the sandwich panel. The area in which the

The panel is loaded with a set of loads that are varying slowly with time, and the analysis is carried out at each load step. The finite element software is set in such a way to solve the model at each load step. This allows all the analysis to be done in a single run of the finite

(a) (b)

2. Face sheets remain elastic all the time:

3. Load scenarios are quasi-static:

non-linear analysis".

into account.

three bonded material layers.

loading for the simply supported panel.

4. Geometric non-linearity has a significant effect:

**3.2. Finite element mesh and boundary conditions** 

distributed load is applied is varying as shown in Figures 5 and 7.

**Figure 7.** The loading area with side length a) 50 mm and b) 200 mm.

**Figure 6.** Stress strain curve for a)material A: AirexR63.50 (Rao, 2002), b) material B: H100 (Kuang, 2001), c)material C: Herex C70.200 (Rao, 2002), d) material D: H250 (Kuang, 2001)

## **3.1. Model assumptions**

All the finite element model analyses done in this research involves the use of non-linear analysis capability of I-DEAS, which includes geometric non-linearity and material nonlinearity. With geometric non-linearity, the software takes the effect of geometry changes into account while calculating the solution. Using material non-linearity option the non-linear behavior of the material response (i.e. post yield material properties) is taken into account.

Below are the assumptions made for the Finite Element Model:

1. Face sheets and core are perfectly bonded:

The numerical model assumes no delamination occur between layers. This assumption is applied by utilizing the partitioning option in the preprocessing module of the software. This option allows the analyst to deal with the whole volume of the structure as one unit also it allows the analyst to assign different material for each partitioned volume.


## **3.2. Finite element mesh and boundary conditions**

362 Finite Element Analysis – New Trends and Developments

**material A**

0 0.05 0.1 0.15 0.2 0.25 **Strain(mm/mm)**

**material C**

advantages of this solid modeling and finite element software.

loading are presented in the coming subsections. The finite element software used in the development of the finite models is (I-DEAS Master Series 10 1999). The relatively robust and user-friendly solid modeling and finite element meshing interface are the main

**material B**

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 **Strain(mm/mm)**

material D

0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 **Strain(mm/mm)**

(d)

(a) (b)

**Stress(MPa)**

**Stress(MPa)**

**Figure 6.** Stress strain curve for a)material A: AirexR63.50 (Rao, 2002), b) material B: H100 (Kuang,

All the finite element model analyses done in this research involves the use of non-linear analysis capability of I-DEAS, which includes geometric non-linearity and material nonlinearity. With geometric non-linearity, the software takes the effect of geometry changes into account while calculating the solution. Using material non-linearity option the non-linear behavior of the material response (i.e. post yield material properties) is taken into

The numerical model assumes no delamination occur between layers. This assumption is applied by utilizing the partitioning option in the preprocessing module of the software. This option allows the analyst to deal with the whole volume of the structure as one unit also it allows the analyst to assign different material for each partitioned

2001), c)material C: Herex C70.200 (Rao, 2002), d) material D: H250 (Kuang, 2001)

Below are the assumptions made for the Finite Element Model:

1. Face sheets and core are perfectly bonded:

(c)

0 0.02 0.04 0.06 0.08 0.1 0.12 **Strain(mm/mm)**

**3.1. Model assumptions** 

account.

0 0.2 0.4 0.6 0.8 1

**stress(MPa)**

**stress(MPa)**

volume.

The symmetric nature of the problem allows only quarter of the whole panel to be meshed. The boundary conditions applied are shown on Figures 4 and 5. The two planes of symmetry of the panel have symmetric boundary conditions, where in-plane displacements and rotation about an axis respective normal to the symmetry plane is allowed. A simply supported boundary condition is applied to the two other sides of the quarter panel. A distributed load is applied on the top surface of the sandwich panel. The area in which the distributed load is applied is varying as shown in Figures 5 and 7.

**Figure 7.** The loading area with side length a) 50 mm and b) 200 mm.

The panel is loaded with a set of loads that are varying slowly with time, and the analysis is carried out at each load step. The finite element software is set in such a way to solve the model at each load step. This allows all the analysis to be done in a single run of the finite

element model. As a result of this, the model would consume less memory space because one single solid model and finite element model can be used for all load steps.

Finite Element Analysis of Loading Area Effect on Sandwich Panel Behaviour Beyond the Yield Limit 365

 One of these challenges is extracting the element force and storing them in a file. This problem is solved by displaying the data on the screen and then copied and stored in a

 Identifying the nodes for a surface of interest so that they can be extracted from the file in which the elemental forces are stored. Since I-DEAS labels the nodes, the nodes corresponding to the surface of interest are copied and stored in a separate node labels

 MATLAB program is developed to extract the elemental forces of the surface of interest from the file in which they are stored by matching the node labels of the surface that are

 Singularity is a serious problem. The post processing analysis for the quality of the elements is utilized to identify the poor elements. The problem is solved by refining the

FEM is capable of capturing the problem details with little approximations compared to

FEM provides solution for many problems like the current case that they do not have

 FEM method is cheap compared to the experimental models. There is no need to produce a prototype or to have high tech facility to conduct the investigation. There is no need for the investigator to be available in a certain place to perform the

The finite element model is verified analytically and experimentally. The analytical verification is based on the classical sandwich panel theory whereas the experimental

Classical sandwich theory has been utilized to obtain close form solution (Zenkret, 1995). The comparison between the numerical and theoretical models in the linear rang are presented in Figure 9. The Figures show very good agreement between theoretical and numerical solution. The classical sandwich plate theory is therefore used to compare and validate the FEM predicted shear distribution of the panel in the linear range. Comparison between the FEM determined shear distribution and the classical sandwich plate theory distribution is performed at all load steps. It is assumed that is the core in the linear range carries the entire shear load. Results obtained from the closed form solution are compared with the total resultant shear load in the global Y direction, *RTOT* (*Yg*), obtained numerically using MATLAB.

The following are some advantages of using FEM over other methods:

separate file for further analysis.

file for further analysis.

stored in node labels file.

the analytical techniques.

investigation is carried out for selected cases.

element size.

*3.2.2. Advantages of FEM* 

analytical solution.

investigation.

**3.3. FEM verification** 

*3.3.1. Analytical verification* 

The numerical model utilizes the map meshing facility in I-DEAS. By controlling the number of nodes along each edge of the solid model, this function provides full control of the mesh size. The element size is chosen by referring to (Miers, 2001) work in mesh refinement. (Mires, 2001) recommended a core element size of 1.5 mm and face element size of 3 mm in order to achieve convergence in the data obtained. For the current case constant mesh density is ensured with the mapped meshing function. This is important because constant mesh density ensures that the data collected from any region in the panel are of the same degree of resolution. Three-dimensional solid brick elements (20 node brick element) are used in this analysis. Second order (parabolic) brick elements are chosen over the first order (linear) brick elements in order to better interpolate the data between nodes. Figure 8 shows the FEM mesh model of the sandwich panel and the brick element utilized in FEM.

**Figure 8.** Illustration of a) Meshed quarter sandwich panel and b) Solid Brick Element (20 node brick element) used in mesh generation.

Since the analysis involves material non-linearity, a yield function or yield criteria needs to be defined for the model. Von Mises yield criteria and its associated flow rule is used in this analysis. Isotropic hardening is also used to describe the change of the yield criterion as a result of plastic straining. Only the core elements are assigned a yield function due to the assumption that only core yielding occurs throughout the loading process. The face sheets are assumed to remain elastic at all time; hence no yield function needs to be assigned to the face sheet elements. However the yield point of the face sheet material is fed to the software to be used as indicator for stopping the analysis.

#### *3.2.1. FEM challenges*

The following challenges are experienced:


## *3.2.2. Advantages of FEM*

364 Finite Element Analysis – New Trends and Developments

element utilized in FEM.

element) used in mesh generation.

*3.2.1. FEM challenges* 

to be used as indicator for stopping the analysis.

The following challenges are experienced:

element model. As a result of this, the model would consume less memory space because

The numerical model utilizes the map meshing facility in I-DEAS. By controlling the number of nodes along each edge of the solid model, this function provides full control of the mesh size. The element size is chosen by referring to (Miers, 2001) work in mesh refinement. (Mires, 2001) recommended a core element size of 1.5 mm and face element size of 3 mm in order to achieve convergence in the data obtained. For the current case constant mesh density is ensured with the mapped meshing function. This is important because constant mesh density ensures that the data collected from any region in the panel are of the same degree of resolution. Three-dimensional solid brick elements (20 node brick element) are used in this analysis. Second order (parabolic) brick elements are chosen over the first order (linear) brick elements in order to better interpolate the data between nodes. Figure 8 shows the FEM mesh model of the sandwich panel and the brick

**Figure 8.** Illustration of a) Meshed quarter sandwich panel and b) Solid Brick Element (20 node brick

(a) (b)

Since the analysis involves material non-linearity, a yield function or yield criteria needs to be defined for the model. Von Mises yield criteria and its associated flow rule is used in this analysis. Isotropic hardening is also used to describe the change of the yield criterion as a result of plastic straining. Only the core elements are assigned a yield function due to the assumption that only core yielding occurs throughout the loading process. The face sheets are assumed to remain elastic at all time; hence no yield function needs to be assigned to the face sheet elements. However the yield point of the face sheet material is fed to the software

one single solid model and finite element model can be used for all load steps.

The following are some advantages of using FEM over other methods:


## **3.3. FEM verification**

The finite element model is verified analytically and experimentally. The analytical verification is based on the classical sandwich panel theory whereas the experimental investigation is carried out for selected cases.

## *3.3.1. Analytical verification*

Classical sandwich theory has been utilized to obtain close form solution (Zenkret, 1995). The comparison between the numerical and theoretical models in the linear rang are presented in Figure 9. The Figures show very good agreement between theoretical and numerical solution. The classical sandwich plate theory is therefore used to compare and validate the FEM predicted shear distribution of the panel in the linear range. Comparison between the FEM determined shear distribution and the classical sandwich plate theory distribution is performed at all load steps. It is assumed that is the core in the linear range carries the entire shear load. Results obtained from the closed form solution are compared with the total resultant shear load in the global Y direction, *RTOT* (*Yg*), obtained numerically using MATLAB. Sample of the total shear resultant comparisons between the numerical and theoretical models in the linear rang are shown Figure 9a and 9b for load of 17.2 kPa and 51.7 kPa respectively.

Finite Element Analysis of Loading Area Effect on Sandwich Panel Behaviour Beyond the Yield Limit 367

2. Fixture for applying simply supported boundary condition is produced. Figure 10

4. Distributed load is applied to the specimen by adaptors manufactured for this purpose.

**Figure 10.** Pictures of the fixture that is produced for applying simply supported boundary condition,

(a) (b)

**Figure 11.** Uniaxial testing machine a) with specimen and b) without specimen

(a) (b)

3. The test is performed on a uniaxial testing machine that is shown in Figure 11.

Figure 12 illustrates the adapters used in experimental setup.

shows two different views of the fixture.

a) top view and b) 3D view.

**Figure 9.** Total plate shear distribution comparison along X-axis at a) 17.2 kPa and b) 51.7 kPa.

## *3.3.2. Experimental verification*

To assure accuracy and validity of the results some selected cases are investigated experimentally. The results obtained from the FEM are compared against those obtained experimentally. Both results show excellent agreement.

### *3.3.2.1. Test setup*

Here is a description of the experimental setup used in the study and consists of the following:

1. The core of the sandwich panel is made of polyurethane foam. Top and bottom sheets of the sandwich panel are made of steel. The dimension of the panels used in the investigation is 250X250 mm2. Mechanical properties of the sheet metal are obtained experimentally.



**Global Y Force (N)**

100

250

**Global Y Force (N)**

400

550

*3.3.2. Experimental verification* 

*3.3.2.1. Test setup* 

experimentally.

following:

experimentally. Both results show excellent agreement.

Sample of the total shear resultant comparisons between the numerical and theoretical models in the linear rang are shown Figure 9a and 9b for load of 17.2 kPa and 51.7 kPa respectively.

(a)

0 50 100 150 200 250 300 350 **Distance from left support (mm)**

Numerical Theoretical

Numerical Teoretical

**Figure 9.** Total plate shear distribution comparison along X-axis at a) 17.2 kPa and b) 51.7 kPa.

To assure accuracy and validity of the results some selected cases are investigated experimentally. The results obtained from the FEM are compared against those obtained

(b)

0 100 200 300 400 **Distance from left support (mm)**

Here is a description of the experimental setup used in the study and consists of the

1. The core of the sandwich panel is made of polyurethane foam. Top and bottom sheets of the sandwich panel are made of steel. The dimension of the panels used in the investigation is 250X250 mm2. Mechanical properties of the sheet metal are obtained 4. Distributed load is applied to the specimen by adaptors manufactured for this purpose. Figure 12 illustrates the adapters used in experimental setup.

**Figure 10.** Pictures of the fixture that is produced for applying simply supported boundary condition, a) top view and b) 3D view.

**Figure 11.** Uniaxial testing machine a) with specimen and b) without specimen

Finite Element Analysis of Loading Area Effect on Sandwich Panel Behaviour Beyond the Yield Limit 369

FEM EXP

**Figure 14.** Comparison of load versus center deflection for core thickness=71 mm, sheet Thickness = 0.5

0 2 4 6 810 **Center Panel Deflection(mm)**

The main advantage of these results over the sandwich panel theory is that both geometric and material nonlinearities are considered without approximation. Usually these approximations eliminate part of the problem physics. By utilizing "I-DEAS' post processing module, stress and its all components, strain and it is all components including

Figures 15a and 15b present Von Mises stress contours for both panel and core respectively whereas Figures 16a and 16b present the plastic strain for both panel and core respectively. It is clear from Figure 16a and 16b that the plastic deformation occurs close to the panel

The criterion, which is adopted by this investigation at what load step the FEM should stop the analysis, is when any of face sheets starts to yield or core material reaches fracture limit. This criterion fulfills the need of the designer; in general design engineer tries to avoid panel face sheets permanent distortion. As soon as the face sheet metal starts to yield, this means that permanent deformation is taking place. So all results produced neither exceed the

Figure 17a present the effect of loading area (area on which the load is applied) for core material A. It is obvious as the loading area increases the stress decreases for the same amount of loading. Same thing can be said for the bottom face sheet in Figure 17b. The core

The effect of loading area at sandwich panels of cores A, B, C and D (see Table 2) are presented 18 through 21. The maximum shear stress of each core in these graphs is

mm, applied load area = 150 X 150 mm2.

the plastic strain, and deformations are obtained.

support (close to the area where boundary conditions are applied).

loading that could cause face - sheet yielding nor exceed core fracture limit.

material (Figure 17a) reaches yield at low loads when the loading area is small.

**4. Results** 

**Load (KN)**

**Figure 12.** The adapters used in the experiments for applying distributed load on specimen of side length a) 200mm, b) 150mm and c) 100mm.

#### *3.3.2.2. Mechanical properties of the specimen*

The sandwich panel is made of polyurethane foam and steel sheets. The mechanical properties are obtained experimentally for both the sheets and the core. ASTM Designation: C 365 – 00 used for testing the core material whereas ASTM Designation: D 638 – 00 used for testing the sheets.

#### *3.3.2.3. Analysis*

The relation between the applied load and the deflection of the specimen center point are shown in Figures 13 and Figure 14 that present a comparison between the experimental results and FEM results. It may be seen that the results are in very good agreement.

To assure accuracy of the experimental results, the experiment is performed many times and the average values are plotted. The variation in the experimental results dose not exceeds 7% of the average value.

**Figure 13.** Comparison of load versus center deflection for core thickness = 49 mm, Sheet Thickness = 0.5 mm, applied load area = 200 X 200 mm2.

**Figure 14.** Comparison of load versus center deflection for core thickness=71 mm, sheet Thickness = 0.5 mm, applied load area = 150 X 150 mm2.

## **4. Results**

FEM EXP

368 Finite Element Analysis – New Trends and Developments

length a) 200mm, b) 150mm and c) 100mm.

testing the sheets.

7% of the average value.

0.5 mm, applied load area = 200 X 200 mm2.

*3.3.2.3. Analysis* 

**Load (KN)**

*3.3.2.2. Mechanical properties of the specimen* 

**Figure 12.** The adapters used in the experiments for applying distributed load on specimen of side

(a) (b) (c)

The sandwich panel is made of polyurethane foam and steel sheets. The mechanical properties are obtained experimentally for both the sheets and the core. ASTM Designation: C 365 – 00 used for testing the core material whereas ASTM Designation: D 638 – 00 used for

The relation between the applied load and the deflection of the specimen center point are shown in Figures 13 and Figure 14 that present a comparison between the experimental

To assure accuracy of the experimental results, the experiment is performed many times and the average values are plotted. The variation in the experimental results dose not exceeds

**Figure 13.** Comparison of load versus center deflection for core thickness = 49 mm, Sheet Thickness =

012345678 **Center Panel Deflection (mm)**

results and FEM results. It may be seen that the results are in very good agreement.

The main advantage of these results over the sandwich panel theory is that both geometric and material nonlinearities are considered without approximation. Usually these approximations eliminate part of the problem physics. By utilizing "I-DEAS' post processing module, stress and its all components, strain and it is all components including the plastic strain, and deformations are obtained.

Figures 15a and 15b present Von Mises stress contours for both panel and core respectively whereas Figures 16a and 16b present the plastic strain for both panel and core respectively. It is clear from Figure 16a and 16b that the plastic deformation occurs close to the panel support (close to the area where boundary conditions are applied).

The criterion, which is adopted by this investigation at what load step the FEM should stop the analysis, is when any of face sheets starts to yield or core material reaches fracture limit. This criterion fulfills the need of the designer; in general design engineer tries to avoid panel face sheets permanent distortion. As soon as the face sheet metal starts to yield, this means that permanent deformation is taking place. So all results produced neither exceed the loading that could cause face - sheet yielding nor exceed core fracture limit.

Figure 17a present the effect of loading area (area on which the load is applied) for core material A. It is obvious as the loading area increases the stress decreases for the same amount of loading. Same thing can be said for the bottom face sheet in Figure 17b. The core material (Figure 17a) reaches yield at low loads when the loading area is small.

The effect of loading area at sandwich panels of cores A, B, C and D (see Table 2) are presented 18 through 21. The maximum shear stress of each core in these graphs is

normalized by the maximum shear yield strength of its corresponding core material and the loading area is normalized by the total surface area of the panel. Also the shear stress of the face sheets is normalized by the corresponding shear yield strength of the face sheets. It can be seen from Table 2, the core materials are labeled from A to D in ascending order according to their stiffness. It is obvious from Figures 18 through 21 that the load carrying capacity of sandwich panel increases by increasing core stiffness. It is observed through all the results that the lower face sheet reaches yield limit before the top face sheet so in the Figures 18 through 21 the lower face sheet is presented. The results of this work are generated according to the univariate search optimization technique (Chapra and Canal, 2006).

Finite Element Analysis of Loading Area Effect on Sandwich Panel Behaviour Beyond the Yield Limit 371

**Figure 16.** Demonstration of the plastic deformations contour for panel A of loading area 4/36 at load

(b)

(a)

As illustrated in Figure 16, the face sheet material starts to yield (entering the plastic range) close to the support (where the boundary conditions are applied). This is physically true, the distributed load over the loading area becomes concentrated reaction force on the strip area on which the boundary conditions (simply supported boundary condition) are applied, i.e., distributed load is converted to concentrated load. So the area where the boundary conditions are applied reaches the yield stress range before any other part of the panel.

As the loading area decreases the load is getting closer to the concentrated load, this is why in Figure 17 panel A of area ratio 1/36 reaches yield (plastic range) at lower load, than the

step 145 kPa for a) the whole panel and b) the core of the panel.

**5. Discussion** 

**Figure 15.** Von Mises stress contour (in MPa) for panel A of loading area 4/36 at load step 145kPa for a) the whole panel and b) the core of the panel.

Finite Element Analysis of Loading Area Effect on Sandwich Panel Behaviour Beyond the Yield Limit 371

**Figure 16.** Demonstration of the plastic deformations contour for panel A of loading area 4/36 at load step 145 kPa for a) the whole panel and b) the core of the panel.

## **5. Discussion**

370 Finite Element Analysis – New Trends and Developments

2006).

normalized by the maximum shear yield strength of its corresponding core material and the loading area is normalized by the total surface area of the panel. Also the shear stress of the face sheets is normalized by the corresponding shear yield strength of the face sheets. It can be seen from Table 2, the core materials are labeled from A to D in ascending order according to their stiffness. It is obvious from Figures 18 through 21 that the load carrying capacity of sandwich panel increases by increasing core stiffness. It is observed through all the results that the lower face sheet reaches yield limit before the top face sheet so in the Figures 18 through 21 the lower face sheet is presented. The results of this work are generated according to the univariate search optimization technique (Chapra and Canal,

**Figure 15.** Von Mises stress contour (in MPa) for panel A of loading area 4/36 at load step 145kPa for a)

(b)

(a)

the whole panel and b) the core of the panel.

As illustrated in Figure 16, the face sheet material starts to yield (entering the plastic range) close to the support (where the boundary conditions are applied). This is physically true, the distributed load over the loading area becomes concentrated reaction force on the strip area on which the boundary conditions (simply supported boundary condition) are applied, i.e., distributed load is converted to concentrated load. So the area where the boundary conditions are applied reaches the yield stress range before any other part of the panel.

As the loading area decreases the load is getting closer to the concentrated load, this is why in Figure 17 panel A of area ratio 1/36 reaches yield (plastic range) at lower load, than the

other panels presented in the Figure. Increasing the loading area increases the load carrying capacity of the panel.

Finite Element Analysis of Loading Area Effect on Sandwich Panel Behaviour Beyond the Yield Limit 373

**Figure 20.** Presentation of maximum shear stress versus loading for A, B, C, and D core material panels

(a) (b)

**Figure 21.** Presentation of maximum shear stress versus loading for A, B, C, and D core material panels

(a) (b)

Figure 19 through 21 present that the lower face sheet for core material B, C and D reaches yield limit before their corresponding core material. This can be referred to the high stiffness

It is obvious from Figure 17 through 21 that panel carrying capacity increases beyond core yield limit. In yield range the core material keeps deforming while the stress is constant (see Figure 22). This deformation works as a mechanism for transferring the excess load to the face sheets. For example in Figure 20, the shear stress of core material A after 100kPa load does not change whereas the shear stress of the corresponding lower face sheet keeps

To replace the core material with same material of the top and bottom sheets, core's width should be shrunk according to the ratio of the modulus of elasticity of the core to that of the metal. The materials B, C and D are relatively stiff in comparison with A. Equivalent crosssection of core material (see Figure 23) has the same height for all cases and the width is increasing according to the modulus of elasticity ratios. For a rectangle the second moment of area (wh3/12) is varying linearly with the width (equivalent width). The effect of the

of its core material, i.e., the panel gets closer in its behavior to isotropic plate.

of load area ratio 16/36 and core thickness 30mm for a) Core and b) Lower Sheet.

of load area ratio 36/36 for a) Core and b) Lower Sheet.

increasing.

**Figure 17.** Presentation of panel A maximum shear stress versus loading for different load area ratio for a) Core and b) Lower Sheet.

**Figure 18.** Presentation of maximum shear stress versus loading for A, B, C, and D core material panels of load area ratio 1/36 for a) Core and b) Lower Sheet.

**Figure 19.** Presentation of maximum shear stress versus loading for A, B, C, and D core material panels of load area ratio 4/36 for a) Core and b) Lower Sheet.

capacity of the panel.

a) Core and b) Lower Sheet.

of load area ratio 1/36 for a) Core and b) Lower Sheet.

of load area ratio 4/36 for a) Core and b) Lower Sheet.

other panels presented in the Figure. Increasing the loading area increases the load carrying

**Figure 17.** Presentation of panel A maximum shear stress versus loading for different load area ratio for

(a) (b)

**Figure 18.** Presentation of maximum shear stress versus loading for A, B, C, and D core material panels

(a) (b)

**Figure 19.** Presentation of maximum shear stress versus loading for A, B, C, and D core material panels

(a) (b)

**Figure 20.** Presentation of maximum shear stress versus loading for A, B, C, and D core material panels of load area ratio 16/36 and core thickness 30mm for a) Core and b) Lower Sheet.

**Figure 21.** Presentation of maximum shear stress versus loading for A, B, C, and D core material panels of load area ratio 36/36 for a) Core and b) Lower Sheet.

Figure 19 through 21 present that the lower face sheet for core material B, C and D reaches yield limit before their corresponding core material. This can be referred to the high stiffness of its core material, i.e., the panel gets closer in its behavior to isotropic plate.

It is obvious from Figure 17 through 21 that panel carrying capacity increases beyond core yield limit. In yield range the core material keeps deforming while the stress is constant (see Figure 22). This deformation works as a mechanism for transferring the excess load to the face sheets. For example in Figure 20, the shear stress of core material A after 100kPa load does not change whereas the shear stress of the corresponding lower face sheet keeps increasing.

To replace the core material with same material of the top and bottom sheets, core's width should be shrunk according to the ratio of the modulus of elasticity of the core to that of the metal. The materials B, C and D are relatively stiff in comparison with A. Equivalent crosssection of core material (see Figure 23) has the same height for all cases and the width is increasing according to the modulus of elasticity ratios. For a rectangle the second moment of area (wh3/12) is varying linearly with the width (equivalent width). The effect of the

difference between the materials B, C, and D is relatively small. So the stress curves for these panels are close to each other and the differences are small as it can be seen in Figures 17 through 21.

Finite Element Analysis of Loading Area Effect on Sandwich Panel Behaviour Beyond the Yield Limit 375

The following are recommendations for further extension of the FEM analysis: Investigate the bonding between the face sheets and the core after yielding.

*Sustainable and Renewable Energy Engineering Program, College of Engineering,* 

*Jordan Aeronautical-Systems Company (JAC), Amman Civil Airport, Marka, Amman, Jordan* 

"C 393-94 (1995). Standard Test Method for Flexural Properties of Sandwich Constructions,

"D 6419-99 (1999). Standard Test Method for Two-Dimensional Flexural Properties of Simply Supported Composite Sandwich Plates Subjected to a Distributed Load", *Annual* 

Airex R63 – (2008), *Typical Mechanical Properties*, Baltek Corporation, www.baltek.com. Allen, H.G., (1969). *Analysis and Design of Structural Sandwich Panels*, Pergamon Press, ISBN

Boyer, H.E. and Gall, T.L. (Eds.). (1991). *Metals Handbook*, American Society for Metals. Caprino, G., Langelan, A. (2000). Study of a Three-Point Bending Specimen for Shear Characterization of Sandwich Cores, *Journal of Composite Materials*, Vol. 34, No. 9, pp.

Chapra, Steve; Canale, Raymond (2005). *Numerical Methods for Engineers*, fifth edition pp.358,

Frostig, Y., Baruch, M., Vilnai, O., Sheinman, I. (1993). Higher-Order Theory for Sandwich Beam Behavior with Transversely Flexible Core, *Journal of Engineering Mechanics*, Vol.

 Modeling face sheets other than metal face sheets such as fiber composite materials Extending the FEM to include the bonding strength between the face sheets and the core so the relative dominance of core shear failure, face sheet yielding, or face sheet

**7. Recommendations** 

**Author details** 

Hussein Maaitah

**8. References** 

791-814.

Corresponding Author

 \*

Salih Akour\*

delamination could be determined.

*University of Sharjah, Sharjah, United Arab Emirates* 

*Mechanical Engineering Department, Amman, Jordan* 

*Annual Book of ASTM Standards*, Vol. 15, No. 3.

*Book of ASTM Standards*, Vol. 15, No. 3.

008012870X 9780080128702, Oxford.

McGraw-Hill, ISBN-10: 0071244298.

119, No. 5, pp. 955-972.

*The University of Jordan, Faculty of Engineering and Technology,* 

**Figure 22.** Schematic drawing of the shear stress for both face sheets and the core within plastic range.

**Figure 23.** Equivalent cross-section of core material with the same height

## **6. Conclusions**


## **7. Recommendations**

The following are recommendations for further extension of the FEM analysis:


## **Author details**

## Salih Akour\*

374 Finite Element Analysis – New Trends and Developments

through 21.

**6. Conclusions** 

experimental one.

face sheet.

capacity.

difference between the materials B, C, and D is relatively small. So the stress curves for these panels are close to each other and the differences are small as it can be seen in Figures 17

**Figure 22.** Schematic drawing of the shear stress for both face sheets and the core within plastic range.

 Investigation of sandwich panel behavior beyond core material yield is carried out. The investigation is accomplished in sight of the core material nonlinearity and the geometric nonlinearity of the whole panel. High tech software 'I-DEAS' (Integrated

 Finite element model is generated using 'I-DEAS' software. This model is validated against experimental and analytical cases available in the literature. To assure model accuracy experimental investigation for selected cases is carried out and compared with FEM. The model shows very good agreement with the analytical as well as the

 It is proved that the load carrying capacity of sandwich panel can be improved by loading the panel beyond the core yield limit. This load is going to be transmitted to the

 Increasing the stiffness of the core material to a certain extent leads to face sheet yielding before the core material. It is proved that increasing core stiffness increases the

 Loading area plays good roll in the load carrying capacity of sandwich panel. Distributing loads over large area of panel surface leads to higher load carrying

Design Engineer Analysis software) is utilized to carry out the investigation.

**Figure 23.** Equivalent cross-section of core material with the same height

load carrying capacity of the sandwich panel.

*Sustainable and Renewable Energy Engineering Program, College of Engineering, University of Sharjah, Sharjah, United Arab Emirates The University of Jordan, Faculty of Engineering and Technology, Mechanical Engineering Department, Amman, Jordan* 

Hussein Maaitah *Jordan Aeronautical-Systems Company (JAC), Amman Civil Airport, Marka, Amman, Jordan* 

## **8. References**


<sup>\*</sup> Corresponding Author

Gdoutos, E.E., Daniel, I.M., Wang, K.-A., Abot, J.L. (June 2001). Non-linear Behavior of Composite Sandwich Beams in Three-point Bending, *Experimental Mechanics*, Vol. 41, No. 2, pp. 182–189.

**Chapter 17** 

© 2012 Weide-Zaage, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

© 2012 Weide-Zaage, licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

**The Finite Element Analysis of Weak** 

**Spots in Interconnects and Packages** 

Concerning the mechanical stress and the electrical-mechanical behavior ULSI multilevel metallization systems are more and more sensitive against influences of geometrical and material changes. The mechanical and electrical reliability of these metallization systems is influenced by this. The reliability of such metallization systems is investigated by thermal and thermal-electrical accelerated stress tests under high temperature load. This leads to degradation due to electro-, thermo- and stress migration. This is one major concern in reliability investigations. Generally measurements are time consuming, expensive and the time-to-market cycle is in the focus of interest too. The prediction of local weak spots in interconnects, vias and solder bumps by finite element simulations are a helpful procedure. Beside this the modern 3-d integration leads to more complex material compositions in the systems concerning the coefficient of thermal expansion (CTE) and other material properties. Higher applied currents on the interconnects and bumps result in Joule heating, high temperature gradients and mechanical stress gradients in the bump and metallization systems. The temperature gradients are much higher compared to systems with wide interconnect lines and bumps with large diameters like in conventional packages. Due to this in ball grid array (BGA) bumps as well as µ-bumps and small through silicon via (TSV) connections current induced migration effects as electromigration (EM) and as result of the high temperature gradients thermomigration (TM) can occur. In interconnects consisting of copper caused by stress gradients due to the different material properties under high

In this chapter the degradation phenomena in dual damascene copper metallization structures as well as degradation in bumps, µ-bumps and TSV are presented. The degradation is current, temperature and mechanical stress induced under a high applied current and temperature load. The finite element analyses and the mass flux divergence

Additional information is available at the end of the chapter

temperature load also stress migration (SM) occurs.

Kirsten Weide-Zaage

http://dx.doi.org/10.5772/50779

**1. Introduction** 


## **Chapter 17**

## **The Finite Element Analysis of Weak Spots in Interconnects and Packages**

Kirsten Weide-Zaage

376 Finite Element Analysis – New Trends and Developments

Thesis, Michigan Technological University.

Ltd., ISBN 0947817778, London.

Services, ISBN 0947817964, United Kingdom.

Science Thesis, Michigan Technological University,.

No. 2, pp. 182–189.

New York.

511-520.

Northwestern University.

Gdoutos, E.E., Daniel, I.M., Wang, K.-A., Abot, J.L. (June 2001). Non-linear Behavior of Composite Sandwich Beams in Three-point Bending, *Experimental Mechanics*, Vol. 41,

Kuang-An, W. (2001). *Failure Analysis of Sandwich Beams*, Doctor of Philosophy Dissertation,

Mercado LL, Sikarskie DL, Miskioglu I. (September 2000). Higher order theory for sandwich beams with yielded core. *Proceedings of ICSS-5 Conference*, Zurich. pp. 141–53. Mercado, L.L., Sikarskie, D.L. (1999). On Response of a Sandwich Panel with a Bilinear Core,

Miers, S.A. (2001). *Analysis and Design of Edge Inserts in Sandwich Beams*, Master of Science

Plantema, F.J. (1966). *Sandwich Construction*, John Wiley and Sons, ISBN-10: 0471691062,

Rao, T. (2002). *Study of Core Compression Using Digital Image Correlation (DIC)*, Master of

Thomsen OT. (1993). Analysis of local bending effects in sandwich plates with orthotropic face layers subjected to localized loads. *Journal of Composite Structure*, vol. 25, no. 1-4, pp.

Thomsen, O.T. (1995). Theoretical and Experimental Investigation of Local Bending Effects

Triantafillou, T.C., Gibson. L.J. (1987). Minimum Weight of Foam Core Sandwich panels for a Given Strength, *Material Science and Engineering*, Vol. 95, November 1987, pp. 55–62. Zenkert, D. (1995). *An Introduction to Sandwich Construction*, (EMAS), The Chameleon Press

Zenkert, D. (1997). *The Handbook of Sandwich Construction*, Engineering Materials Advisory

*Mechanics of Composite Materials and Structures*, Vol. 6, No. 1, pp.57–67.

in Sandwich Plates, *Composite Structures*, Vol. 30, No. 1. pp. 85-101.

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/50779

## **1. Introduction**

Concerning the mechanical stress and the electrical-mechanical behavior ULSI multilevel metallization systems are more and more sensitive against influences of geometrical and material changes. The mechanical and electrical reliability of these metallization systems is influenced by this. The reliability of such metallization systems is investigated by thermal and thermal-electrical accelerated stress tests under high temperature load. This leads to degradation due to electro-, thermo- and stress migration. This is one major concern in reliability investigations. Generally measurements are time consuming, expensive and the time-to-market cycle is in the focus of interest too. The prediction of local weak spots in interconnects, vias and solder bumps by finite element simulations are a helpful procedure. Beside this the modern 3-d integration leads to more complex material compositions in the systems concerning the coefficient of thermal expansion (CTE) and other material properties. Higher applied currents on the interconnects and bumps result in Joule heating, high temperature gradients and mechanical stress gradients in the bump and metallization systems. The temperature gradients are much higher compared to systems with wide interconnect lines and bumps with large diameters like in conventional packages. Due to this in ball grid array (BGA) bumps as well as µ-bumps and small through silicon via (TSV) connections current induced migration effects as electromigration (EM) and as result of the high temperature gradients thermomigration (TM) can occur. In interconnects consisting of copper caused by stress gradients due to the different material properties under high temperature load also stress migration (SM) occurs.

In this chapter the degradation phenomena in dual damascene copper metallization structures as well as degradation in bumps, µ-bumps and TSV are presented. The degradation is current, temperature and mechanical stress induced under a high applied current and temperature load. The finite element analyses and the mass flux divergence

© 2012 Weide-Zaage, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2012 Weide-Zaage, licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

calculation of these phenomena will show the suitability of the method in comparison with experimental results. For a prediction of the weakest spot the suitability of the finite element mesh as well as the modeling concerning the structure shape has to be investigated. Especially edges in the model can influence the quality of the results. The geometrical data of the different metallization or package structures can be taken from the layout as well as optical, scanning microscope or other analytical techniques. Especially for the prediction of the electromigration induced weakest spot in the system, the location of the maximum current density is a major indicator for the fault location. Due to this the maximum current density at localizations of structure inhomogeneity must be checked up. Out of this the maximum current density is calculated by conformal mapping to predict an optimized modeling concerning the shape of the edges and the use of a radius instead of edges. Geometrical variations like the thickness of the first and second metallization and a comparison of the different migration mechanisms will be presented. Concerning the mechanical stress of the DD-Cu metallization the process induced stress will be considered with different processing temperatures of the copper metallization. The influence of different dielectrics on the mechanical stress is also determined. Compared to DDmetallizations and the traces, the bumps are only exposed by electro- and thermomigration. The thermal-electrical-mechanical behavior of µ-bumps and TSV will be shown for a Waferon-Wafer (WoW) structure.

The Finite Element Analysis of Weak Spots in Interconnects and Packages 379

*Hk TB D De* (1)

The diffusion coefficient with the activation enthalpy of ∆H is defined as follows:

[Heumann, Philibert] applies:

/

(2)

0

The change of place takes place at cubic (fcc) metals such as Cu, Au and Ag across the holes. The activation enthalpy is additively composed of the enthalpy of formation and the activation energy for the hole formation. Both are nearly the same size. With equation (1) and D0 as temperature independent diffusion constant as well as the melting temperature Tm

> 6 4 3 <sup>0</sup> 5 10 5 10 [ ² / ], 1.5 10 [ ] *D m s dH T eV <sup>m</sup>*

> > Parameter Al Cu Ag Tm 660 1083 960 ∆V/Ω 0.71-1.3 0.9 0.7-0.9 E [eV] 1.234 2.15 1.921

In the case of multiple components, for example solder bumps, all main materials must be involved in the material transport. In the example of a material of two components A and B, and the chemical or interdiffusion coefficients, there are the partial diffusion coefficients DA

> *A B B A AB <sup>N</sup> D D with D D x D x*

(3)

Due to this a diffusion coefficient of D(Tm)≈10-8 cm2/s [Heumann] is resulting.

**Table 1.** Melting Temperature, Active Volume and Activation Energy [Philibert].

thermotransport moves the Pb in Sn to the cold part of the sample.

predominantly interface or surface migration is occurring.

determined.

and DB [Wilkenson]. Using the Darken equation, the diffusion results as follows:

<sup>~</sup> *<sup>A</sup>*

*x* 

From the migration velocity υ, the chemical diffusion coefficient ܦ ෩and the mole fraction, xA and xB of the considered components the partial diffusion coefficients DA and DB can be

The direction of electrical transport in alloys such as leaded solder materials based on Sn, goes, depending on the proportion of Pb to the anode or the cathode. In terms of

There are different diffusion paths in the material. These are the grain boundaries, surfaces to adjacent materials and the volume or bulk. In the aluminum as a metallization material primarily grain boundary diffusion occurs. In contrary in copper as metallization material

In the calculation of grain boundaries, the surface diffusion as well as the intermediate phase diffusion must be considered in thin layers. Atomistic transport along a grain boundary or phase boundary has a low activation energy and is therefore by orders of magnitude faster than in the crystal itself. The analytical models for the grain boundary diffusion δ DGB are only valid for the self diffusion in pure metals [Kaur]. The part of

## **2. Calculation of the migration mechanism**

To figure out the possibilities of determining the effects of migration in solder bumps, interconnects, via and conductive pathways shall be shown here. Migration, particularly in copper or aluminum metallization or migration of solder bumps of a flip-chip package, are presented in [Banas, Hou, Liu 2007, Liu2008, Tan, Wang]. Simulation algorithms are based on analogies and allow only partly- or no material or sizing variations.

Basically, migration processes are described in the metallization on the consideration of diffusion processes. These diffusion processes can lead to a change in dimension or change of geometries of the metallization [Wever]. On one side material loss or the hole formation at these locations leads to tensile stress, while it comes in places of material accumulation to a compressive stress. This leads to a reflux, also called the back stress in the material. There is a critical length existing, which compensates the electromigration related flux and the reflux. The migration of atoms leads to a change of the expansion in the metal, which then in turn changes the chemical potential. Since the chemical potential is very often influenced by a stress term, this leads to stress gradients and a reflux. In microelectronic applications the current flux pathways are subject to large thermal strain before any current load applies. This thermal strain, along with the strain caused by the electromigration can lead to a nonlinear reflux. Also, in the interconnects local self-heating take place, providing a contribution to the migration of the thermal gradients due to the local temperature increases. The components of the electrical, thermal and stress migration flux are added by superposition.

The diffusion coefficient with the activation enthalpy of ∆H is defined as follows:

378 Finite Element Analysis – New Trends and Developments

on-Wafer (WoW) structure.

superposition.

**2. Calculation of the migration mechanism** 

on analogies and allow only partly- or no material or sizing variations.

calculation of these phenomena will show the suitability of the method in comparison with experimental results. For a prediction of the weakest spot the suitability of the finite element mesh as well as the modeling concerning the structure shape has to be investigated. Especially edges in the model can influence the quality of the results. The geometrical data of the different metallization or package structures can be taken from the layout as well as optical, scanning microscope or other analytical techniques. Especially for the prediction of the electromigration induced weakest spot in the system, the location of the maximum current density is a major indicator for the fault location. Due to this the maximum current density at localizations of structure inhomogeneity must be checked up. Out of this the maximum current density is calculated by conformal mapping to predict an optimized modeling concerning the shape of the edges and the use of a radius instead of edges. Geometrical variations like the thickness of the first and second metallization and a comparison of the different migration mechanisms will be presented. Concerning the mechanical stress of the DD-Cu metallization the process induced stress will be considered with different processing temperatures of the copper metallization. The influence of different dielectrics on the mechanical stress is also determined. Compared to DDmetallizations and the traces, the bumps are only exposed by electro- and thermomigration. The thermal-electrical-mechanical behavior of µ-bumps and TSV will be shown for a Wafer-

To figure out the possibilities of determining the effects of migration in solder bumps, interconnects, via and conductive pathways shall be shown here. Migration, particularly in copper or aluminum metallization or migration of solder bumps of a flip-chip package, are presented in [Banas, Hou, Liu 2007, Liu2008, Tan, Wang]. Simulation algorithms are based

Basically, migration processes are described in the metallization on the consideration of diffusion processes. These diffusion processes can lead to a change in dimension or change of geometries of the metallization [Wever]. On one side material loss or the hole formation at these locations leads to tensile stress, while it comes in places of material accumulation to a compressive stress. This leads to a reflux, also called the back stress in the material. There is a critical length existing, which compensates the electromigration related flux and the reflux. The migration of atoms leads to a change of the expansion in the metal, which then in turn changes the chemical potential. Since the chemical potential is very often influenced by a stress term, this leads to stress gradients and a reflux. In microelectronic applications the current flux pathways are subject to large thermal strain before any current load applies. This thermal strain, along with the strain caused by the electromigration can lead to a nonlinear reflux. Also, in the interconnects local self-heating take place, providing a contribution to the migration of the thermal gradients due to the local temperature increases. The components of the electrical, thermal and stress migration flux are added by

$$\begin{array}{rcl}D&=&D\_0\,e^{-\Lambda H/k\_B T} \end{array} \tag{1}$$

The change of place takes place at cubic (fcc) metals such as Cu, Au and Ag across the holes. The activation enthalpy is additively composed of the enthalpy of formation and the activation energy for the hole formation. Both are nearly the same size. With equation (1) and D0 as temperature independent diffusion constant as well as the melting temperature Tm [Heumann, Philibert] applies:

$$5 \cdot 10^{-6} < D\_0 < 5 \cdot 10^{-4} \text{ [}m^2/\text{s}\text{]} \quad dH \equiv 1.5 \cdot 10^{-3} \text{ }T\_m \text{ [}eV\text{]} \tag{2}$$

Due to this a diffusion coefficient of D(Tm)≈10-8 cm2/s [Heumann] is resulting.


**Table 1.** Melting Temperature, Active Volume and Activation Energy [Philibert].

In the case of multiple components, for example solder bumps, all main materials must be involved in the material transport. In the example of a material of two components A and B, and the chemical or interdiffusion coefficients, there are the partial diffusion coefficients DA and DB [Wilkenson]. Using the Darken equation, the diffusion results as follows:

$$
\delta\nu = \left(D\_A - D\_B\right) \frac{\delta N\_A}{\delta \mathbf{x}} \quad \text{with} \quad \tilde{D} = D\_B \mathbf{x}\_A + D\_A \mathbf{x}\_B \tag{3}
$$

From the migration velocity υ, the chemical diffusion coefficient ܦ ෩and the mole fraction, xA and xB of the considered components the partial diffusion coefficients DA and DB can be determined.

The direction of electrical transport in alloys such as leaded solder materials based on Sn, goes, depending on the proportion of Pb to the anode or the cathode. In terms of thermotransport moves the Pb in Sn to the cold part of the sample.

There are different diffusion paths in the material. These are the grain boundaries, surfaces to adjacent materials and the volume or bulk. In the aluminum as a metallization material primarily grain boundary diffusion occurs. In contrary in copper as metallization material predominantly interface or surface migration is occurring.

In the calculation of grain boundaries, the surface diffusion as well as the intermediate phase diffusion must be considered in thin layers. Atomistic transport along a grain boundary or phase boundary has a low activation energy and is therefore by orders of magnitude faster than in the crystal itself. The analytical models for the grain boundary diffusion δ DGB are only valid for the self diffusion in pure metals [Kaur]. The part of

electromigration at grain boundaries depends also on the effective width of the grain δ (d) concerning the mass transportation in relation to the average grain size.

A composition of the interconnect from almost poly crystalline parts and parts where bamboo structured grains occur, leads to the predominant grain diffusion on one hand and on the other to volume diffusion and thus to high gradients in the mass flux. If the mass flux is blocked by a large grain, then the material accumulates relating to direction of the electron or material flux, in front of the blocking grain, while it comes behind the blocking grain to a material loss. For grain boundaries in the area of 280nm the diffusion coefficient of the grain boundary diffusion is in the range of the bulk material DGB ≅ DBulk [Kaur].

Only above a defined threshold current density ��th, resulting from the Blech effect and denoted as 'short-length effect', it comes to effective place change and thus to the mass transportation or material flux [Blech]. The mass flux is dependent on the atomic particle density N, the Boltzmann constant kB, the local temperature T, the current density ��, the specific resistance ρ, the diffusion process relating activation energy EA, the diffusion coefficient D0, and the effective charge eZ.

$$\overrightarrow{J\_{EM}} = \frac{N}{k\_B T} e Z^\* \left(\overrightarrow{j} - \overrightarrow{j}\_{th}\right) \rho \,\,\, D\_0 \exp\left(-\frac{E\_A}{k\_B T}\right) \tag{4}$$

The Finite Element Analysis of Weak Spots in Interconnects and Packages 381

(6)

<sup>0</sup> exp *<sup>A</sup> SM H B B <sup>N</sup> <sup>E</sup> <sup>J</sup> <sup>D</sup> grad kT k T*

The Ω represents the atomic volume and σH the hydrostatic stress. Out of the components of the electrical, thermal, and stress migration, the total mass flux in the metallization structure

The formation of material accumulation or void formation requires a divergence in the flux of the total mass flux. Only divergences lead to a change in the density of the material. The calculation of the mass flux and mass flux divergence in metallization, traces and bumps in the past were done using analogies between the electrical and thermal and the thermal mechanical behavior or a variation of the materials or dimension were not possible [Banas, Hou, Liu 2007, Liu2008, Tan, Wang]. This lead to an overestimation of the temperature gradients and current density like proposed in [Ogurtani]. Without any applied current and

With the new program code the divgrad of T and σ are calculated directly based on the simulation results. The calculation of the mass flux is done for each element under consideration of the neighbor elements. Out of this the stress gradients are calculated. The divgrad terms are calculated on base of the super elements under consideration of an Ansatzfunction. The calculated values for the different migration mechanisms are reloaded into ANSYS for graphic display. Out of this the stress migration can be calculated under SM stress test conditions without any applied current or temperature gradients. The simplified equations neglecting concentration gradients for the different migration mechanism are:

emerges through superposition from the equations (4, 5 and 6).

������� ������� � � ��

������� ������� � � ��

�������������� � � ��

**3. Modeling and simulation improvements** 

temperature gradients the calculation of the stress migration is not possible.

���� <sup>+</sup> ��

���� <sup>−</sup> � �

���� <sup>−</sup> � �

���������) <sup>−</sup> �

In equation (7-9) N is the atomic concentration, the Boltzmann constant kB, the local temperature T, the resistivity ρ, the atomic volume Ω, Q\* the heat of transport, the activation energy EA (taken from grain boundary and interface migration as strongest influence) and σH is the hydrostatic stress. The simulation and calculation sequence is shown in figure 1.

The numbers of elements, which are determining the mesh and due to this the density of the nodes have a strong influence on the accuracy of the simulation. Out of this the mesh of the investigated metallization or metallic material as main interesting point in the simulation plays a major rule. In the metallization itself the potential as well as the temperature is calculated. The metallization is surrounded by dielectric material and the traces by FR-4 or PCB. The bumps can be surrounded by underfillers made of different plastic materials.

�

����� ������� � ������ − ���

������������ � ������ − ���

���

����� ������� � ������ (7)

���� � �� (8)

� ��� (9)

The thermotransport denoted as Soret-effect is a phenomenon of overlay of diffusion and heat conduction. They closely resemble the electrotransport with the difference that the considered system is not isothermal [Wever]. There is a flux of soluted atoms and the heat flux. These fluxes are described about the chemical potential and the thermal gradients. Without the occurrence of concentration gradients the equation for thermal flux is:

$$\overrightarrow{J\_{TM}} = -\frac{N\mathcal{Q}^\*}{k\_B T^2} \,^D D\_0 \exp\left(-\frac{E\_A}{k\_B T}\right) \text{grad}\, T \tag{5}$$

Q\* represents here the transported energy at a constant temperature and is commonly referred to as reduced heat of transport or transfer. The ratio Q\* / kBT² is referred as Soret coefficient. The heat of transport is the heat flux per unit of material without temperature gradients. Is the value of Q\* > 0 a heat flux is generated to keep the soluted atoms isothermal, which takes place towards the dissolved flux. Is Q\* < 0 the flux of dissolved particles and the heat flux are counter set. It follows that in an isothermal system, a density gradient produces a thermal flux and vice versa a temperature gradient leads to a material flux [Shewman]. The heat of transport is approximately equal to the activation energy for the material flux [Jaffe].

The differences of the coefficient of thermal expansion (CTE) between the metallization material and surrounding materials produce, depending on the ambient temperature a mechanical stress. This in turn leads to a material flux within the metallization. Under strain the enthalpy for example of tensile-stress formation in a grain boundary is reduced, which increases the concentration of holes [Heumann].

The Finite Element Analysis of Weak Spots in Interconnects and Packages 381

$$\overrightarrow{J\_{SM}} = -\frac{N\Omega}{k\_B T} \, D\_0 \exp\left(-\frac{E\_A}{k\_B T}\right) \text{grad}\,\sigma\_H \tag{6}$$

The Ω represents the atomic volume and σH the hydrostatic stress. Out of the components of the electrical, thermal, and stress migration, the total mass flux in the metallization structure emerges through superposition from the equations (4, 5 and 6).

380 Finite Element Analysis – New Trends and Developments

coefficient D0, and the effective charge eZ.

the material flux [Jaffe].

electromigration at grain boundaries depends also on the effective width of the grain δ (d)

A composition of the interconnect from almost poly crystalline parts and parts where bamboo structured grains occur, leads to the predominant grain diffusion on one hand and on the other to volume diffusion and thus to high gradients in the mass flux. If the mass flux is blocked by a large grain, then the material accumulates relating to direction of the electron or material flux, in front of the blocking grain, while it comes behind the blocking grain to a material loss. For grain boundaries in the area of 280nm the diffusion coefficient of the grain

Only above a defined threshold current density ��th, resulting from the Blech effect and denoted as 'short-length effect', it comes to effective place change and thus to the mass transportation or material flux [Blech]. The mass flux is dependent on the atomic particle density N, the Boltzmann constant kB, the local temperature T, the current density ��, the specific resistance ρ, the diffusion process relating activation energy EA, the diffusion

<sup>0</sup> \* exp *<sup>A</sup>*

(4)

(5)

*B B <sup>N</sup> <sup>E</sup> J eZ j j D k T k T* 

\* exp *<sup>A</sup>*

*B B NQ <sup>E</sup> <sup>J</sup> <sup>D</sup> gradT k T k T* 

Q\* represents here the transported energy at a constant temperature and is commonly referred to as reduced heat of transport or transfer. The ratio Q\* / kBT² is referred as Soret coefficient. The heat of transport is the heat flux per unit of material without temperature gradients. Is the value of Q\* > 0 a heat flux is generated to keep the soluted atoms isothermal, which takes place towards the dissolved flux. Is Q\* < 0 the flux of dissolved particles and the heat flux are counter set. It follows that in an isothermal system, a density gradient produces a thermal flux and vice versa a temperature gradient leads to a material flux [Shewman]. The heat of transport is approximately equal to the activation energy for

The differences of the coefficient of thermal expansion (CTE) between the metallization material and surrounding materials produce, depending on the ambient temperature a mechanical stress. This in turn leads to a material flux within the metallization. Under strain the enthalpy for example of tensile-stress formation in a grain boundary is reduced, which

The thermotransport denoted as Soret-effect is a phenomenon of overlay of diffusion and heat conduction. They closely resemble the electrotransport with the difference that the considered system is not isothermal [Wever]. There is a flux of soluted atoms and the heat flux. These fluxes are described about the chemical potential and the thermal gradients.

Without the occurrence of concentration gradients the equation for thermal flux is:

2 0

concerning the mass transportation in relation to the average grain size.

boundary diffusion is in the range of the bulk material DGB ≅ DBulk [Kaur].

*EM th*

*TM*

increases the concentration of holes [Heumann].

The formation of material accumulation or void formation requires a divergence in the flux of the total mass flux. Only divergences lead to a change in the density of the material. The calculation of the mass flux and mass flux divergence in metallization, traces and bumps in the past were done using analogies between the electrical and thermal and the thermal mechanical behavior or a variation of the materials or dimension were not possible [Banas, Hou, Liu 2007, Liu2008, Tan, Wang]. This lead to an overestimation of the temperature gradients and current density like proposed in [Ogurtani]. Without any applied current and temperature gradients the calculation of the stress migration is not possible.

With the new program code the divgrad of T and σ are calculated directly based on the simulation results. The calculation of the mass flux is done for each element under consideration of the neighbor elements. Out of this the stress gradients are calculated. The divgrad terms are calculated on base of the super elements under consideration of an Ansatzfunction. The calculated values for the different migration mechanisms are reloaded into ANSYS for graphic display. Out of this the stress migration can be calculated under SM stress test conditions without any applied current or temperature gradients. The simplified equations neglecting concentration gradients for the different migration mechanism are:

$$\overrightarrow{div}\,\overrightarrow{\nu\_{EM}} = \left(\frac{E\_A}{k\_B T^2} + \frac{a\_\rho}{1 + a\_\rho (T - T\_0)} - \frac{1}{T}\right) \cdot \overrightarrow{f\_{EM}} \cdot \overrightarrow{grad}\,T\tag{7}$$

$$\overrightarrow{\nabla}\,d\upsilon\,\overrightarrow{J\_{TM}} = \left(\frac{E\_A}{k\_B T^2} - \frac{2}{T}\right) \cdot \overrightarrow{J\_{TM}} \cdot \overrightarrow{grad\ T} - \frac{\text{QND}}{k\_B T^2} \cdot \Delta T \tag{8}$$

$$\overrightarrow{div}\,\overrightarrow{\upsilon\_{SM}} = \left(\frac{E\_A}{k\_B T^2} - \frac{1}{T}\right) \cdot \overrightarrow{f\_{SM}} \cdot \overrightarrow{grad}\,\,T - \frac{\Omega ND}{k\_B T} \cdot \Delta\sigma\_H \tag{9}$$

In equation (7-9) N is the atomic concentration, the Boltzmann constant kB, the local temperature T, the resistivity ρ, the atomic volume Ω, Q\* the heat of transport, the activation energy EA (taken from grain boundary and interface migration as strongest influence) and σH is the hydrostatic stress. The simulation and calculation sequence is shown in figure 1.

#### **3. Modeling and simulation improvements**

The numbers of elements, which are determining the mesh and due to this the density of the nodes have a strong influence on the accuracy of the simulation. Out of this the mesh of the investigated metallization or metallic material as main interesting point in the simulation plays a major rule. In the metallization itself the potential as well as the temperature is calculated. The metallization is surrounded by dielectric material and the traces by FR-4 or PCB. The bumps can be surrounded by underfillers made of different plastic materials.

The Finite Element Analysis of Weak Spots in Interconnects and Packages 383

like ANSYS®, ABAQUS or COMSOL using the finite element analyses. In the following

In the thermal electrical calculation the degree of freedom is only the temperature. The mesh of the surrounding material depends in the investigated case of the metallization and due to this of the metallization mesh. A coarse mesh of the metallization can lead to an insufficient determination of the temperature in the isolation materials. Out of this the mesh of the

The mesh can be easily refined by varying the count of the nodes. If a convergence concerning the change of the potential voltage is found the mesh is sufficient. A predication concerning the convergence of the current density at inhomogeneities like vias or edges in the metallization cannot be done. At these positions a strong dependence of the maximum current density due to current crowding in relation to the model of this inhomogeneity occurs. The layout from the EDA tools give at such positions rectangular edges or edges with a defined angle. Taking into account such an angle in the model may lead to insufficient calculation of the local current density at this position. Due to this the relation of

Considering the two-dimensional case, choosing a right angle at the edge of a metallization a singularity will occur. In this case the current density can be calculated analytical with an infinity value [Betz]. The result of the numerical analysis of the current density at this position converges to a defined value due to the fact that the neighbor elements are taken into account for the calculation. If the right angle is replaced by a radius, the calculated current density depends on the selected shape of the radius, the element size and the element shape itself. For a determination of the real current density in the structure and the optimized radius, the edge of a simplified metallization consisting of aluminum, is investigated for the 2-dimensional and 3-dimensional case in sub-models. The boundaries in

In the 2-dimensional as well as in the 3-dimensional case the current crowding jmax/jin increases with increasing element count cubically and the relation of the maximum current density to the applied current density jmax/jin increases with decreasing element size exponentially. Out of this the numerical solution converges in the case of infinite small elements to the theoretical expected value. In the case of a curvature of the edges, the current crowding jmax/jin has for increasing element counts as well as decreasing element size for every taken radius a constant value. Due to this the element count used in the simulations depends on the element size as well as the used radius in the model. A decreasing radius down to the range of an rectangular angle with constant element count

Not in every case optical or scanning electron microscopy pictures are available for a determination of the structure shapes after the production process. Due to this for a

metallization should be investigated with regard of an optimization.

the current crowding and angle of the edge will be investigated here.

the calculations were taken of the coarse model.

leads to an increase of jmax/jin.

investigation ANSYS® is used as simulation tool.

**3.2. Mesh density and singularities** 

**Figure 1.** Simulation and calculation sequence for one cycle (static) and more cycle (dynamic).

Concerning the mechanical stress calculation the interfaces play also an important rule. Due to the different materials and their material properties the stress in the structures can be high and can therefore have a strong influence on the reliability of the structures. The calculation of the stress gradients and stress migration mass flux divergence is strongly affected by the accuracy of the mesh.

## **3.1. Advantages of FEM compared with other applicable methods**

Simulation methods like finite difference, finite element and Monte Carlo method are common in the electrical engineering. Monte Carlo method is described for the calculation of time to failure (TTF) distributions induced by electromigration [Huang]. Also in the scope of semiconductor simulation it is widely used [Jacoboni]. The complexity of the metallization and bump structures in combination with the question of multiphysic investigations makes this approach inconvenient in use. The application of the finite difference method is described in investigations of the intermetallic phase growth for instance [Chao]. A model for electromigration calculated by concentration gradients is described in [Joo] and the void evolution and motion was investigated by [Averbuch]. The mass transport along interface caused by electromigration using the level-set method is described in [Li]. All these approaches consider only a part of the thermal-electrical and mechanical behavior.

Calculations in the scope of the coupled thermal-electrical-mechanical behavior of metallization structures and bumps can be sufficiently carried out by commercial programs like ANSYS®, ABAQUS or COMSOL using the finite element analyses. In the following investigation ANSYS® is used as simulation tool.

## **3.2. Mesh density and singularities**

382 Finite Element Analysis – New Trends and Developments

affected by the accuracy of the mesh.

mechanical behavior.

**Figure 1.** Simulation and calculation sequence for one cycle (static) and more cycle (dynamic).

**3.1. Advantages of FEM compared with other applicable methods** 

Concerning the mechanical stress calculation the interfaces play also an important rule. Due to the different materials and their material properties the stress in the structures can be high and can therefore have a strong influence on the reliability of the structures. The calculation of the stress gradients and stress migration mass flux divergence is strongly

Simulation methods like finite difference, finite element and Monte Carlo method are common in the electrical engineering. Monte Carlo method is described for the calculation of time to failure (TTF) distributions induced by electromigration [Huang]. Also in the scope of semiconductor simulation it is widely used [Jacoboni]. The complexity of the metallization and bump structures in combination with the question of multiphysic investigations makes this approach inconvenient in use. The application of the finite difference method is described in investigations of the intermetallic phase growth for instance [Chao]. A model for electromigration calculated by concentration gradients is described in [Joo] and the void evolution and motion was investigated by [Averbuch]. The mass transport along interface caused by electromigration using the level-set method is described in [Li]. All these approaches consider only a part of the thermal-electrical and

Calculations in the scope of the coupled thermal-electrical-mechanical behavior of metallization structures and bumps can be sufficiently carried out by commercial programs In the thermal electrical calculation the degree of freedom is only the temperature. The mesh of the surrounding material depends in the investigated case of the metallization and due to this of the metallization mesh. A coarse mesh of the metallization can lead to an insufficient determination of the temperature in the isolation materials. Out of this the mesh of the metallization should be investigated with regard of an optimization.

The mesh can be easily refined by varying the count of the nodes. If a convergence concerning the change of the potential voltage is found the mesh is sufficient. A predication concerning the convergence of the current density at inhomogeneities like vias or edges in the metallization cannot be done. At these positions a strong dependence of the maximum current density due to current crowding in relation to the model of this inhomogeneity occurs. The layout from the EDA tools give at such positions rectangular edges or edges with a defined angle. Taking into account such an angle in the model may lead to insufficient calculation of the local current density at this position. Due to this the relation of the current crowding and angle of the edge will be investigated here.

Considering the two-dimensional case, choosing a right angle at the edge of a metallization a singularity will occur. In this case the current density can be calculated analytical with an infinity value [Betz]. The result of the numerical analysis of the current density at this position converges to a defined value due to the fact that the neighbor elements are taken into account for the calculation. If the right angle is replaced by a radius, the calculated current density depends on the selected shape of the radius, the element size and the element shape itself. For a determination of the real current density in the structure and the optimized radius, the edge of a simplified metallization consisting of aluminum, is investigated for the 2-dimensional and 3-dimensional case in sub-models. The boundaries in the calculations were taken of the coarse model.

In the 2-dimensional as well as in the 3-dimensional case the current crowding jmax/jin increases with increasing element count cubically and the relation of the maximum current density to the applied current density jmax/jin increases with decreasing element size exponentially. Out of this the numerical solution converges in the case of infinite small elements to the theoretical expected value. In the case of a curvature of the edges, the current crowding jmax/jin has for increasing element counts as well as decreasing element size for every taken radius a constant value. Due to this the element count used in the simulations depends on the element size as well as the used radius in the model. A decreasing radius down to the range of an rectangular angle with constant element count leads to an increase of jmax/jin.

Not in every case optical or scanning electron microscopy pictures are available for a determination of the structure shapes after the production process. Due to this for a determination of the optimized radius two dimensional simulations can be compared with results from the calculation of the maximum current density out of conformal mapping. The maximum current density and details about the homogeneity of the current density and resistance behavior at the investigated places can be achieved by the conformal mapping.

The calculation of the maximum current density is done by the following approach. In general for the analytical calculation of the current density distribution at the edge of an interconnection, a bent metallization with different width g and h and a radius r can be converted by two times conformal mapping into a flat conducting band with width π in the Z' plane. In this flat conducting band a homogenous current density distribution can be assumed (figure 2). With the back transformation of this homogenous current density distribution into the Z plane the current density distribution along the path A-B-C can be achieved and due to this the maximum current density. With the parameters S=g/h as ratio of the metallization width and P=r/h as ratio of the radius r to the smaller metallization width h and the applied current density jin in the vertical part of the metallization shall apply under the condition the P << 1 [Hagedorn]:

$$\dot{j}\_{\text{max}} = \mathbf{1}, \mathbf{04} \begin{bmatrix} \frac{\mathbf{s}^2 + \mathbf{1}}{\mathbf{P} \cdot \mathbf{s}^2} \end{bmatrix}^{\frac{1}{3}} \dot{\mathbf{j}}\_{\text{in}} \tag{10}$$

The Finite Element Analysis of Weak Spots in Interconnects and Packages 385

a good correspondence between simulation and analytical solution is found. For a radius above 30nm the simulated values are caused by the increased element size above the analytical values. A radius below 30nm has an insufficient mesh and due to this the simulated values are beneath the analytical values. To get compliance between simulation and analytical solution the mesh has to be refined. This leads to an increase of the

**Figure 3.** Max. current density depending on the radius for the interconnect/via edge. Analytical

In the finite element method simulation of the mechanical stress, the calculated unknown parameters are the displacement of the nodes. The different forces are calculated by the derivation of the displacements. The form function is continuously but not smoothed between the elements. Due to this the values at the nodes have to be averaged. Due to this and caused by the simplified form function and the limited element expansion the FEM is an approximation method. This can lead to failures in the calculation of the dilatations. When the results of an inaccurate mesh are derived the inaccuracy concerning a second derivation will be intensified. The correlation of the accuracy and mesh refinement is shown in figure 4 [Eichelseder]. From this knowledge it is clear that the distribution of the dilatation of the components S1 to S3 is more accurate than the strain distribution resp. the Von Mises Stress (VMS). This fact also implies that a simplified mesh or a coarse mesh can deliver accurate results concerning the dilatation. The derivate strain shows accurate values only after a refinement. After an additional refinement the gradient of the strain and out of this the hydrostatic stress (HS) can be calculated with a good accuracy. The calculation of the mass flux divergence due to stress migration includes also a derivation of the stress gradients, which will affect the accuracy additional. The accuracy of the finite element mesh has to be investigated before the mass flux divergence calculation is done. And the results have to be

calculation time and is therefore not useful.

solution done by conformal mapping and simulation.

**3.3. Mechanical stress and stress gradients** 

handled with care if the mesh is not adequate tested.

With equation 10 the current crowding jmax/jin can be determined depending on the metallization width and the radius.

**Figure 2.** Original figure of a two dimensional edge of an metallization in the Z-plane (left). Metallization after two times of conformal mapping in the Z' plane (right).

In the 2-dimensional case a good correlation between the analytical solution and the maximum current density determined by the simulations can be found. In the 3-dimensional case the analytical solution is valid under the assumption of a homogenous current density distribution in the via, which means that the current density distribution can be mapped on the 2-dimensional case.

The maximum current density depending on the radius for a fixed element count is compared to the simulated values. The results are shown in figure 3. For a radius of 0.03µm a good correspondence between simulation and analytical solution is found. For a radius above 30nm the simulated values are caused by the increased element size above the analytical values. A radius below 30nm has an insufficient mesh and due to this the simulated values are beneath the analytical values. To get compliance between simulation and analytical solution the mesh has to be refined. This leads to an increase of the calculation time and is therefore not useful.

**Figure 3.** Max. current density depending on the radius for the interconnect/via edge. Analytical solution done by conformal mapping and simulation.

#### **3.3. Mechanical stress and stress gradients**

384 Finite Element Analysis – New Trends and Developments

apply under the condition the P << 1 [Hagedorn]:

metallization width and the radius.

the 2-dimensional case.

determination of the optimized radius two dimensional simulations can be compared with results from the calculation of the maximum current density out of conformal mapping. The maximum current density and details about the homogeneity of the current density and resistance behavior at the investigated places can be achieved by the conformal mapping.

The calculation of the maximum current density is done by the following approach. In general for the analytical calculation of the current density distribution at the edge of an interconnection, a bent metallization with different width g and h and a radius r can be converted by two times conformal mapping into a flat conducting band with width π in the Z' plane. In this flat conducting band a homogenous current density distribution can be assumed (figure 2). With the back transformation of this homogenous current density distribution into the Z plane the current density distribution along the path A-B-C can be achieved and due to this the maximum current density. With the parameters S=g/h as ratio of the metallization width and P=r/h as ratio of the radius r to the smaller metallization width h and the applied current density jin in the vertical part of the metallization shall

(10)

max in 1,04 j

With equation 10 the current crowding jmax/jin can be determined depending on the

*P S j*

**Figure 2.** Original figure of a two dimensional edge of an metallization in the Z-plane (left).

In the 2-dimensional case a good correlation between the analytical solution and the maximum current density determined by the simulations can be found. In the 3-dimensional case the analytical solution is valid under the assumption of a homogenous current density distribution in the via, which means that the current density distribution can be mapped on

The maximum current density depending on the radius for a fixed element count is compared to the simulated values. The results are shown in figure 3. For a radius of 0.03µm

Metallization after two times of conformal mapping in the Z' plane (right).

*S*

 

> In the finite element method simulation of the mechanical stress, the calculated unknown parameters are the displacement of the nodes. The different forces are calculated by the derivation of the displacements. The form function is continuously but not smoothed between the elements. Due to this the values at the nodes have to be averaged. Due to this and caused by the simplified form function and the limited element expansion the FEM is an approximation method. This can lead to failures in the calculation of the dilatations. When the results of an inaccurate mesh are derived the inaccuracy concerning a second derivation will be intensified. The correlation of the accuracy and mesh refinement is shown in figure 4 [Eichelseder]. From this knowledge it is clear that the distribution of the dilatation of the components S1 to S3 is more accurate than the strain distribution resp. the Von Mises Stress (VMS). This fact also implies that a simplified mesh or a coarse mesh can deliver accurate results concerning the dilatation. The derivate strain shows accurate values only after a refinement. After an additional refinement the gradient of the strain and out of this the hydrostatic stress (HS) can be calculated with a good accuracy. The calculation of the mass flux divergence due to stress migration includes also a derivation of the stress gradients, which will affect the accuracy additional. The accuracy of the finite element mesh has to be investigated before the mass flux divergence calculation is done. And the results have to be handled with care if the mesh is not adequate tested.

The Finite Element Analysis of Weak Spots in Interconnects and Packages 387

the metallization height here was varied. In the investigated cases for electromigration the current flux direction was calculated downstream from the second metallization M2 to the first metallization M1 and in the opposite direction upstream. The mass flux divergence of the variation is shown in figure 5 (right). It was found that the variation of the first metallization M1 and the variation of the second metallization M2 show a big influence on the mass flux divergence and out of this on the reliability of the structure. The weakest

For technology nodes in the range of 500nm and below SiO2 or combinations of USG (undoped-silicate-glass), PSG (phosphosilicated-glass) or FSG (fluorine-doped-silicate-glass) as a customary dielectric (IMD) is used. For products produced in nowadays low-κ dielectrics, with a low dielectric constant such as Silk ™, SiCOH, Black Diamond, or MSQ ™

Figure 6 shows the maximum hydrostatic stress as a function of the applied current density for the different dielectric materials. The greatest mechanical stress occurs at a SiCOH

A second model is based on dimensions of the 65nm technology node with a wide line and different via bottom geometries. The model is shown in figure 7 based on SEM pictures from the literature [Delsol, Lee]. Tantalum and TaN were chosen as barrier material, SiCN as capping material and SiCOH was chosen as dielectric material. The structures were investigated under test conditions with an applied current density of 1.5MA/cm² taken from [Lin2006a, Lin2006b]. The thermo-mechanical investigations were carried out considering

The stress state in the copper metallization is strongly related to the process temperature of the copper. A high process temperature leads to a higher nearly stress free state of the whole structure compared to structures with lower temperatures [Matsuyama]. The nearly stress

location found here is in and beneath the via.

**Figure 6.** Maximum hydrostatic Stress for BDII™, Silk und SiCOH.

**4.2. Separated migration mechanism with different via shape**

the process-induced stress. The process temperatures are given in table 2.

II are used.

dielectric and the lowest at Silk ™.

**Figure 4.** Correlation of the accuracy and mesh refinement [Eichelseder].

## **4. Metallization**

## **4.1. Metallization variation in the 90nm node**

The investigated DD copper structure based on the dimension of the 90nm scaled down to the 65nm process technology. The model consists of SiCN as cap layer, Ta/TaN as barrier layer and different dielectrics like Silk™, Black Diamond II™ and SiCOH. In the mechanical calculations the process temperature using the birth and die algorithm in ANSYS are included, caused by the fact that the use of a reference temperature for the stress free state is not sufficient [Weide-Zaage 2008].

**Figure 5.** Mesh of the investigated structure (left) and mass flux divergence vs metallization height for a variation of M1 and M2 (right).

The mesh of the investigated structure is shown in figure 5 (left). For simplification one via with two metal links right and left hand was generated. As one example for this proposal the metallization height here was varied. In the investigated cases for electromigration the current flux direction was calculated downstream from the second metallization M2 to the first metallization M1 and in the opposite direction upstream. The mass flux divergence of the variation is shown in figure 5 (right). It was found that the variation of the first metallization M1 and the variation of the second metallization M2 show a big influence on the mass flux divergence and out of this on the reliability of the structure. The weakest location found here is in and beneath the via.

**Figure 6.** Maximum hydrostatic Stress for BDII™, Silk und SiCOH.

386 Finite Element Analysis – New Trends and Developments

**4. Metallization** 

not sufficient [Weide-Zaage 2008].

variation of M1 and M2 (right).

**Figure 4.** Correlation of the accuracy and mesh refinement [Eichelseder].

The investigated DD copper structure based on the dimension of the 90nm scaled down to the 65nm process technology. The model consists of SiCN as cap layer, Ta/TaN as barrier layer and different dielectrics like Silk™, Black Diamond II™ and SiCOH. In the mechanical calculations the process temperature using the birth and die algorithm in ANSYS are included, caused by the fact that the use of a reference temperature for the stress free state is

**Figure 5.** Mesh of the investigated structure (left) and mass flux divergence vs metallization height for a

The mesh of the investigated structure is shown in figure 5 (left). For simplification one via with two metal links right and left hand was generated. As one example for this proposal

**4.1. Metallization variation in the 90nm node** 

For technology nodes in the range of 500nm and below SiO2 or combinations of USG (undoped-silicate-glass), PSG (phosphosilicated-glass) or FSG (fluorine-doped-silicate-glass) as a customary dielectric (IMD) is used. For products produced in nowadays low-κ dielectrics, with a low dielectric constant such as Silk ™, SiCOH, Black Diamond, or MSQ ™ II are used.

Figure 6 shows the maximum hydrostatic stress as a function of the applied current density for the different dielectric materials. The greatest mechanical stress occurs at a SiCOH dielectric and the lowest at Silk ™.

## **4.2. Separated migration mechanism with different via shape**

A second model is based on dimensions of the 65nm technology node with a wide line and different via bottom geometries. The model is shown in figure 7 based on SEM pictures from the literature [Delsol, Lee]. Tantalum and TaN were chosen as barrier material, SiCN as capping material and SiCOH was chosen as dielectric material. The structures were investigated under test conditions with an applied current density of 1.5MA/cm² taken from [Lin2006a, Lin2006b]. The thermo-mechanical investigations were carried out considering the process-induced stress. The process temperatures are given in table 2.

The stress state in the copper metallization is strongly related to the process temperature of the copper. A high process temperature leads to a higher nearly stress free state of the whole structure compared to structures with lower temperatures [Matsuyama]. The nearly stress

free state determined by investigating of the process induced stress is approx. 200°C. The calculation of the stress migration in the structures is described in [Weide-Zaage 2010].

The Finite Element Analysis of Weak Spots in Interconnects and Packages 389

The failure locations due to the void formation after the electromigration test in the literature are found downstream at the bottom of the via and upstream in and over the via at the interface metallization cap layer. A schematically overview of the different migration effects like electro-, thermo- and stress migration mass flux electron flux upstream and downstream and supposed void formation (yellow) is drawn in the via structure shown in figure 6. The electromigration mass flux is depending on the current flux direction and temperature gradients, the thermomigration mass flux depend on the temperature gradients and the stress migration on the stress gradients, which occur at tensile regions also under electromigration test conditions. An arbitrary assumed grain distribution is grey colored in this graphic. Thermomigration mass flux (blue) in the investigated structures proceeds independent of the current flux direction from first as well as the second metallization into the via. The electromigration mass flux (green) is at the upper part of the via downstream in the same direction and at the via bottom in opposite direction of the temperature gradients and thermomigration mass flux. In the upstream case it is vice versa. Due to the process induced stress distribution in the metallization, under electromigration test conditions, the metallization is mostly compressive with some tensile regions. The stress gradients in the tensile regions may lead to a current direction independent migration. High stress migration (orange) is found at the bottom and beneath the via above and below the barrier. The stress gradients show above the barrier upwards and below the barrier downwards. At the interface of the cap layer and the metallization they show down. Only the occurrence of a migration pathway leads to voiding. For interface migration the stress gradients have to have a component into the direction of the interface. Also the existence of a grain boundary pathway supports the voiding. Out of this foot voiding as well as voids in the via itself, both

**Figure 8.** Schematically electro-, thermo- and stress migration mass flux electron flux upstream and

The thermomigration (Soret-Effect) was investigated under electromigration test conditions of 325°C. The calculation of the thermomigration shows high, values for the mass flux divergence. With temperature gradients of 50K/µm the gradients are high but not high

downstream and supposed void formation (yellow) in the via structure EM>SM>TM.

indicated in yellow can be explained by this (figure 8).

*4.2.2. Thermomigration stress temperature 325°C*

**Figure 7.** SEM pictures from [Lin2006, Lin2006(2)] (left) and Finite Element Mesh of the via region of the four different bottom geometries (right).


**Table 2.** Mechanical properties and process temperatures used in the simulations

## *4.2.1. Electromigration test temperature 325°C*

The electromigration behavior of the four models was investigated with current (flux in both directions) and an activation energy of 0.9eV. The different via bottom geometries were verified by comparison with investigations from literature. The calculated maximum mass flux divergences as well as the reciprocal values which are related to the MTF are given in table 3**.** Comparing a Gouging and a Flat via and a Cone Shape and a V-Gouging via the Gouging via is more reliable than the Flat via and the Cone shape via more reliable than the V-Gouging via. The simulated models show that the Cone Shape via has the best results concerning the electromigration behavior.


**Table 3.** Mass flux divergence and reciprocal mass flux divergence.

The failure locations due to the void formation after the electromigration test in the literature are found downstream at the bottom of the via and upstream in and over the via at the interface metallization cap layer. A schematically overview of the different migration effects like electro-, thermo- and stress migration mass flux electron flux upstream and downstream and supposed void formation (yellow) is drawn in the via structure shown in figure 6. The electromigration mass flux is depending on the current flux direction and temperature gradients, the thermomigration mass flux depend on the temperature gradients and the stress migration on the stress gradients, which occur at tensile regions also under electromigration test conditions. An arbitrary assumed grain distribution is grey colored in this graphic. Thermomigration mass flux (blue) in the investigated structures proceeds independent of the current flux direction from first as well as the second metallization into the via. The electromigration mass flux (green) is at the upper part of the via downstream in the same direction and at the via bottom in opposite direction of the temperature gradients and thermomigration mass flux. In the upstream case it is vice versa. Due to the process induced stress distribution in the metallization, under electromigration test conditions, the metallization is mostly compressive with some tensile regions. The stress gradients in the tensile regions may lead to a current direction independent migration. High stress migration (orange) is found at the bottom and beneath the via above and below the barrier. The stress gradients show above the barrier upwards and below the barrier downwards. At the interface of the cap layer and the metallization they show down. Only the occurrence of a migration pathway leads to voiding. For interface migration the stress gradients have to have a component into the direction of the interface. Also the existence of a grain boundary pathway supports the voiding. Out of this foot voiding as well as voids in the via itself, both indicated in yellow can be explained by this (figure 8).

**Figure 8.** Schematically electro-, thermo- and stress migration mass flux electron flux upstream and downstream and supposed void formation (yellow) in the via structure EM>SM>TM.

#### *4.2.2. Thermomigration stress temperature 325°C*

388 Finite Element Analysis – New Trends and Developments

the four different bottom geometries (right).

*4.2.1. Electromigration test temperature 325°C*

concerning the electromigration behavior.

**Table 3.** Mass flux divergence and reciprocal mass flux divergence.

free state determined by investigating of the process induced stress is approx. 200°C. The calculation of the stress migration in the structures is described in [Weide-Zaage 2010].

**Figure 7.** SEM pictures from [Lin2006, Lin2006(2)] (left) and Finite Element Mesh of the via region of

Cu 125 0.34 16.7 10-6 200 TaN 185 0.33 6.6 10-6 40 SiCN 100 0.17 3 10-6 300 SiCOH 15 0.3 11.6 10-6 350 Si 98 0.45 2.64 10-6 25

**Table 2.** Mechanical properties and process temperatures used in the simulations

Young Modul (GPa) Poisson CTE (300K) Process Temperature (°C)

The electromigration behavior of the four models was investigated with current (flux in both directions) and an activation energy of 0.9eV. The different via bottom geometries were verified by comparison with investigations from literature. The calculated maximum mass flux divergences as well as the reciprocal values which are related to the MTF are given in table 3**.** Comparing a Gouging and a Flat via and a Cone Shape and a V-Gouging via the Gouging via is more reliable than the Flat via and the Cone shape via more reliable than the V-Gouging via. The simulated models show that the Cone Shape via has the best results

> Via-Shape **div EM (div EM)-1**  Gouging Via 1.143 10-3 875 Flat Via 0.917 10-3 1091 V-Gouging Via 0.939 10-3 1065 Cone Via 0.672 10-3 1488

The thermomigration (Soret-Effect) was investigated under electromigration test conditions of 325°C. The calculation of the thermomigration shows high, values for the mass flux divergence. With temperature gradients of 50K/µm the gradients are high but not high enough to induce thermomigration. Gradients of 100K/µm are supposed to be an acceleration factor for the thermomigration. In figure 9 the mass flux divergence distribution due to thermomigration is shown. High values in the Gouging model are found at the via bottom in the via and at the corner of the metallization, barrier- and cap-layer. For the Cone Shape via high values are found only at the corner of the corner the metallization, barrierand cap-layer. The activation energy for the thermomigration calculation is supposed to be too small and should be >1.2eV. Measurements of thermomigration activation energies of copper in copper are not known until now to verify this.

The Finite Element Analysis of Weak Spots in Interconnects and Packages 391

**Figure 10.** Model of the µ-bump (left). The different materials are indicated by colors. The mass flux

In figure 10 (right) the mass flux divergence depending on the applied current for the SAC PoP bump and the µ-bump is shown. Under the same applied current the mass flux divergence of µ-bumps is compared to PoP bumps about four orders of magnitude higher. Due to this fact the reliability of the µ-bump has a high electromigration risk. The weakest point of the µ-bump was found at the top of the bump due to current crowding combined

**Figure 11.** Temperature distribution in a BGA PoP Package with ICs (left) and maximum temperature

[Meinshausen 2010(b), Meinshausen 2012] investigated in previous simulations the heat flux density of IC stacks with a constant surface load. The heat loss of the ICs is one reason for the existence of an inhomogeneous temperature distribution in packages. An applied power of 2-8W for all ICs in the IC stack was investigated as boundary condition for the simulations. The temperature distribution in the package for a power loss of 8W is shown in figure 11, left side. The maximum temperature occurs in the IC stack. The bump temperature is much lower (figure11, right). This leads to a strong inhomogeneous temperature distribution. On the other hand the temperature gradients in the bumps are not as high compared to the EM simulation. Due to this in this investigated case

divergence vs. the applied current for SAC bump and µ-bump (right).

with high temperature gradients.

in the bumps vs applied IC's power loss (right).

thermomigration will not occur.

**Figure 9.** Thermomigration mass flux divergence distribution for the Gouging and the Cone Shape model.

## **5. µ-Bump, CuSn-pillar and TSV**

## **5.1. µ-Bump in comparison with BGA-PoP**

A variation of the applied current in a Package-on-Package (PoP) bumps [Meinshausen 2010 (a)] and µ-bump [Meinshausen 2011] was carried out and the mass flux divergence distribution was determined. The bumps in the FE model of the PoP device with one bump chain consist of SAC305 (SnAgCu). For the under-bump metallization (UBM) and the surface finishes 6µm thick Ni layers were used. As mold compound (MC) around the upper contact surfaces "Stycast 1090" was chosen. The simulations were carried out with anisotropic and temperature depending material parameters. In the FR-4 substrates five layers of Cu traces and Cu vias between the top and the bottom package are placed to connect the top and the bottom bumps. The height of the upper, the middle and the lower copper layers is 20µm, 18µm and 36µm. The width of all layers is 250µm. The model of a µbump between two ICs is shown in figure 10 (left). The dimensions of the µ-bumps are similar to the test structures used in [Labie]. The diameter of the µ- bump is 25µm and the height is 10µm. Over and under the µ-bump a 100µm silicon layer resp. a 50µm thick silicon layer is representing the ICs of a CoC (Chip-on-Chip) structure. The ICs are covered with a 1µm thick Si3N4 passivation layer. The copper traces at the upper and the lower contact surface have a height of 0.5µm and a width of 32µm. The pitch is 40µm.

**5. µ-Bump, CuSn-pillar and TSV** 

**5.1. µ-Bump in comparison with BGA-PoP** 

model.

copper in copper are not known until now to verify this.

enough to induce thermomigration. Gradients of 100K/µm are supposed to be an acceleration factor for the thermomigration. In figure 9 the mass flux divergence distribution due to thermomigration is shown. High values in the Gouging model are found at the via bottom in the via and at the corner of the metallization, barrier- and cap-layer. For the Cone Shape via high values are found only at the corner of the corner the metallization, barrierand cap-layer. The activation energy for the thermomigration calculation is supposed to be too small and should be >1.2eV. Measurements of thermomigration activation energies of

**Figure 9.** Thermomigration mass flux divergence distribution for the Gouging and the Cone Shape

A variation of the applied current in a Package-on-Package (PoP) bumps [Meinshausen 2010 (a)] and µ-bump [Meinshausen 2011] was carried out and the mass flux divergence distribution was determined. The bumps in the FE model of the PoP device with one bump chain consist of SAC305 (SnAgCu). For the under-bump metallization (UBM) and the surface finishes 6µm thick Ni layers were used. As mold compound (MC) around the upper contact surfaces "Stycast 1090" was chosen. The simulations were carried out with anisotropic and temperature depending material parameters. In the FR-4 substrates five layers of Cu traces and Cu vias between the top and the bottom package are placed to connect the top and the bottom bumps. The height of the upper, the middle and the lower copper layers is 20µm, 18µm and 36µm. The width of all layers is 250µm. The model of a µbump between two ICs is shown in figure 10 (left). The dimensions of the µ-bumps are similar to the test structures used in [Labie]. The diameter of the µ- bump is 25µm and the height is 10µm. Over and under the µ-bump a 100µm silicon layer resp. a 50µm thick silicon layer is representing the ICs of a CoC (Chip-on-Chip) structure. The ICs are covered with a 1µm thick Si3N4 passivation layer. The copper traces at the upper and the lower contact

surface have a height of 0.5µm and a width of 32µm. The pitch is 40µm.

**Figure 10.** Model of the µ-bump (left). The different materials are indicated by colors. The mass flux divergence vs. the applied current for SAC bump and µ-bump (right).

In figure 10 (right) the mass flux divergence depending on the applied current for the SAC PoP bump and the µ-bump is shown. Under the same applied current the mass flux divergence of µ-bumps is compared to PoP bumps about four orders of magnitude higher. Due to this fact the reliability of the µ-bump has a high electromigration risk. The weakest point of the µ-bump was found at the top of the bump due to current crowding combined with high temperature gradients.

**Figure 11.** Temperature distribution in a BGA PoP Package with ICs (left) and maximum temperature in the bumps vs applied IC's power loss (right).

[Meinshausen 2010(b), Meinshausen 2012] investigated in previous simulations the heat flux density of IC stacks with a constant surface load. The heat loss of the ICs is one reason for the existence of an inhomogeneous temperature distribution in packages. An applied power of 2-8W for all ICs in the IC stack was investigated as boundary condition for the simulations. The temperature distribution in the package for a power loss of 8W is shown in figure 11, left side. The maximum temperature occurs in the IC stack. The bump temperature is much lower (figure11, right). This leads to a strong inhomogeneous temperature distribution. On the other hand the temperature gradients in the bumps are not as high compared to the EM simulation. Due to this in this investigated case thermomigration will not occur.

## **5.2. TSV with µ-bump**

Through silicon vias (TSV) have a wide range of applications in the modern packaging. Three different kinds of TSV can be identified which are common under the topic of 3D IC packaging, 3D IC integration and 3D Si integration. Due to this the processing as well as the size of the TSV is completely different. The application can be divided into three topics.

The Finite Element Analysis of Weak Spots in Interconnects and Packages 393

On top of the TSV high VMS was found caused by the fact that the copper in the TSV is under pressure and the material is pressed out of the TSV. Concerning the EM high mass fluxes were found at the bottom of the TSV to the metal trace. At this position also high thermal mass flux occurs. High hydrostatic stress occurs at the bottom of the TSV to the metal trace. The stress increases about 10% for an increase of the applied current from 0.1 to

For a simplified assembly, a smaller pitch between the bumps and a low interconnect inductance CuSn-pillars can be used. Copper pillars are copper bumps with a thin layer of Sn on the top. These layer thicknesses can vary [Huffman, Syed]. The CuSn-pillars can be formed under pressure and a temperature load. This leads to the formation of Cu3Sn and Cu6Sn5 phases. Due to this CuSn pillars with different Sn thickness and location in the bump were investigated. The bumps had a diameter of 24µm and a complete high of 45µm. Above the bumps a Cu trace and below the bumps a TSV (figure 12) were placed in the model. The applied current was set to 175 mA, the substrate temperature was varied from -50 to 150°C and the stress free temperature was set to 135°C. In figure 13 a part of the mesh and the

**Figure 13.** Mesh of the CuSn Pillar with Cu in magenta and Sn in yellow (above) and temperature

gradient distribution for the different models (below).

temperature gradient distribution in the different bumps is shown.

0.3A.

**5.3. CuSn-pillar** 


In the case of the 3D IC integration the TSVs are necessary. Due to the fact that the data width is limited, the use of TSV with small sizes in the range of 5-10µm with a pitch in the range of 20-40µm a much wider I/O width is possible. The TSV shortens the way of the connections compared to the common wire bonding. Using the wire bonding the chip size has to be increased and due to this the costs increase. TSV as solution will solve this problem in future. For future 3D integration Cu-TSV as well as tungsten W-TSV will be used [Murugesan] caused by the fact that W is a capable material for sub µ-vias but not usable for power and ground applications due to its high specific resistance. Therefore Cu with a Ta barrier will be used.

**Figure 12.** Schematic picture of a TSV and µ-bump (left) and FE-Model (right).

As an example in figure 12 (right) the mesh of a TSV with a SAC pillar is shown. The right side of figure 12 shows schematically a TSV with a µ-bump [Leduc]. The geometrical data for the simulations are taken from [Kitada, Lo]. The diameter of the TSV was set to 9µm and the height was set to 80µm. The copper barrier in the TSV consists of SiN and TiN. The thickness of the interconnections is 7µm. The passivation consists of SiN.

On top of the TSV high VMS was found caused by the fact that the copper in the TSV is under pressure and the material is pressed out of the TSV. Concerning the EM high mass fluxes were found at the bottom of the TSV to the metal trace. At this position also high thermal mass flux occurs. High hydrostatic stress occurs at the bottom of the TSV to the metal trace. The stress increases about 10% for an increase of the applied current from 0.1 to 0.3A.

## **5.3. CuSn-pillar**

392 Finite Element Analysis – New Trends and Developments

the chips; active or passive interposer.

Through silicon vias (TSV) have a wide range of applications in the modern packaging. Three different kinds of TSV can be identified which are common under the topic of 3D IC packaging, 3D IC integration and 3D Si integration. Due to this the processing as well as the size of the TSV is completely different. The application can be divided into three topics.

b. Chip to Chip, Chip to Wafer or Wafer to Wafer stacking with SiO2 to Sio2 and Cu to Cu

c. Wafer to Wafer stacking; Memory chip stacking TSV are used as connection between

In the case of the 3D IC integration the TSVs are necessary. Due to the fact that the data width is limited, the use of TSV with small sizes in the range of 5-10µm with a pitch in the range of 20-40µm a much wider I/O width is possible. The TSV shortens the way of the connections compared to the common wire bonding. Using the wire bonding the chip size has to be increased and due to this the costs increase. TSV as solution will solve this problem in future. For future 3D integration Cu-TSV as well as tungsten W-TSV will be used [Murugesan] caused by the fact that W is a capable material for sub µ-vias but not usable for power and ground applications due to its high specific resistance. Therefore Cu with a Ta

a. Package on Package stacking with TSV as connections between them

**Figure 12.** Schematic picture of a TSV and µ-bump (left) and FE-Model (right).

thickness of the interconnections is 7µm. The passivation consists of SiN.

As an example in figure 12 (right) the mesh of a TSV with a SAC pillar is shown. The right side of figure 12 shows schematically a TSV with a µ-bump [Leduc]. The geometrical data for the simulations are taken from [Kitada, Lo]. The diameter of the TSV was set to 9µm and the height was set to 80µm. The copper barrier in the TSV consists of SiN and TiN. The

**5.2. TSV with µ-bump** 

bonding.

barrier will be used.

For a simplified assembly, a smaller pitch between the bumps and a low interconnect inductance CuSn-pillars can be used. Copper pillars are copper bumps with a thin layer of Sn on the top. These layer thicknesses can vary [Huffman, Syed]. The CuSn-pillars can be formed under pressure and a temperature load. This leads to the formation of Cu3Sn and Cu6Sn5 phases. Due to this CuSn pillars with different Sn thickness and location in the bump were investigated. The bumps had a diameter of 24µm and a complete high of 45µm. Above the bumps a Cu trace and below the bumps a TSV (figure 12) were placed in the model. The applied current was set to 175 mA, the substrate temperature was varied from -50 to 150°C and the stress free temperature was set to 135°C. In figure 13 a part of the mesh and the temperature gradient distribution in the different bumps is shown.

**Figure 13.** Mesh of the CuSn Pillar with Cu in magenta and Sn in yellow (above) and temperature gradient distribution for the different models (below).

The homogenous temperature distribution the bumps can be achieved by a placement of the Sn in the middle (model A). In the case of model B-C high temperature gradients occur above the Sn near the copper trace. At this position also current crowding occurs. Both can lead to a weak link at this position. Depending on the current flow direction the flux will be increased or decreased. In figure 14 the hydrostatic stress depending of the substrate temperature is shown for model A the highest stress occurs for temperatures about -50°C near the reference temperature of the stress free state the stress is the lowest. The processing temperatures should be included in the model.

The Finite Element Analysis of Weak Spots in Interconnects and Packages 395

The finite element analysis of metallization structures or packages and bumps provide an insight into local temperatures, current density and stress gradient distributions with the possibility of mass flux divergence calculation. A prediction of weak links in the structures helps to increase the reliability during the design phase. Also the costs will be decreased and

In future phase separation in the CuSn bumps by IMC growth as well as the influence of special placed grain boundaries and concentration gradients should be included in the model. The process induced stress should be also included in the package

Averbuch, A.; Israeli, M.; Ravve, I.: "Electromigration of intergranular voids in metal films for microelectronic interconnects", Journal of Computational Physics 186, 2, 2003, pp.

Banas, L.; Nürnberg, R.: "Finite Element Approximation of a Three Dimensional Phase Field Model for Void Electromigration", Springer, J. Sci. Computation, DOI 10.1007/s10915-

Blech, I.A.: "Electromigration in thin aluminum films on titanium nitride", J. Appl. Phys. 47,

Delsol, J. R.; Jacquemin, J.-P.; et. al.: "Improved electrical and reliability performance of 65 nm interconnects with new barrier integration schemes", Microelectronic Engineering

Chao, B; Chae, S.H.; et.al: "Electromigration enhanced intermetallic growth and void

Eichelseder, W.; Unger, B.; Dannbauer, H.: "Neue Aspekte der Submodelltechnik in der finite Elemente Methode zur Beschleunigung des Entwicklungsprozesses", FEM-

Hagedorn, F.B.; Hall, P.M.: "Right-Angle Bends in Thin Film Conductors", Journal Appl.

Heumann, T.: " Diffusion in Metallen", Springer Verlag, WFT 10, ISBN 3-540-55379-7, 1992.

*Gottfried Wilhelm Leibniz Universität Hannover, Information Technology Laboratory,* 

Betz, A.: "Konforme Abbildungen", Berlin, Springer Verlag, 1964, pp.22-33.

formation in Pb-free solder joints", J. App. Phys. 100, 084909 (2006)

**6. Conclusion** 

redesigns avoided.

**Author details** 

**7. References** 

281–304.

008-9203-y.

1976, p. 1203.

83 (2006) pp. 2377–2380.

Konferenz Baden Baden, 1998.

Phys., Vol. 34, No.1, 1963, pp. 128-131.

Kirsten Weide-Zaage

modeling.

*Germany* 


**Table 4.** Hydrostatic Stress, Maximum Temperature and Maximum Mass flux Divergence.

**Figure 14.** Hydrostatic stress depending of the substrate temperature (model A).

## **6. Conclusion**

394 Finite Element Analysis – New Trends and Developments

temperatures should be included in the model.

The homogenous temperature distribution the bumps can be achieved by a placement of the Sn in the middle (model A). In the case of model B-C high temperature gradients occur above the Sn near the copper trace. At this position also current crowding occurs. Both can lead to a weak link at this position. Depending on the current flow direction the flux will be increased or decreased. In figure 14 the hydrostatic stress depending of the substrate temperature is shown for model A the highest stress occurs for temperatures about -50°C near the reference temperature of the stress free state the stress is the lowest. The processing

> B 174 396,86 0,486e12 C 175 396,54 0,459e12 D 172 397,02 0,512e12

**Table 4.** Hydrostatic Stress, Maximum Temperature and Maximum Mass flux Divergence.

**Figure 14.** Hydrostatic stress depending of the substrate temperature (model A).

HS [MPa] Tmax [K] Div Jmax [a.u.]

The finite element analysis of metallization structures or packages and bumps provide an insight into local temperatures, current density and stress gradient distributions with the possibility of mass flux divergence calculation. A prediction of weak links in the structures helps to increase the reliability during the design phase. Also the costs will be decreased and redesigns avoided.

In future phase separation in the CuSn bumps by IMC growth as well as the influence of special placed grain boundaries and concentration gradients should be included in the model. The process induced stress should be also included in the package modeling.

## **Author details**

Kirsten Weide-Zaage

*Gottfried Wilhelm Leibniz Universität Hannover, Information Technology Laboratory, Germany* 

## **7. References**


Hou, Y; Tan, C.M.: "Comparison of stress-induced voiding phenomena in copper line-via structures with different dielectric materials", Semiconductor Sci. Technology, 2009, pp.1-8.

The Finite Element Analysis of Weak Spots in Interconnects and Packages 397

Lo, W.-C. ; Chen, Y.-H. ; et.al: TSV and 3D Wafer Bonding Technologies For Advanced Stacking System and Application at ITRI, Symp. on VLSI Tech. Digest of Technical

Matsuyama, H.; et al.: "Investigation of stress-induced voiding inside and under vias in copper interconnects with "wing" pattern", Proc. Conf. International Reliability Physics

Meinshausen, L.(a); Weide-Zaage, K.; et.al: ";Virtual prototyping of PoP interconnections regarding electrically activated mechanisms" IEEE Thermal, Mechanical and Multi-Physics Simulation and Experiments in Microelectronics and Micro-Systems, EuroSimE

Meinshausen, L. (b); Weide-Zaage, K.: "Exploration of Migration and Stress Effects in PoPs Considering Inhomogeneous Temperature Distribution", International Wafer-Level

Meinshausen, L.; Weide-Zaage, K.; et.al: "Electro- and Thermomigration in Microbump Interconnects for 3D Integration", IEEE Electronic Components and Technology Conf.,

Meinshausen, L.; Weide-Zaage, K.; et.al:"Thermal Management for stackable package with

Murugesan, M; Kino, H.; et.al: High Density §D LSI Technology using W/CU Hybrid TSVs,

Ogurtani, T.O.; Akyildiz, O.: "Morphological evolution of voids by surface drift diffusion driven by capillary, electromigration, and thermal-stress gradients induced by steadystate heat flow in passivated metallic thin films and flip chip solder joints. I. Theory,

Philibert, J.: "Atom Movements Diffusion and Mass Transport in Solids", Monographies de

Shewman, P. G.: "Diffusion in Solids", Mac Graw-Hill Series in Materials Science and

Syed, A.; Dhandapani, K. et.al: "Cu Pillar and µ-bump Electromigration Reliability and Comparison with High Pb, SnPb, and SnAg bumps, IEEE Electronic Components and

Tan, C.M.; Roy, A.; et.al: "Current Crowding Effect on Copper Dual Damascene Via Bottom Failure for ULSI Applications" IEEE Transactions Dev. Materials Reliability 2005; Vol.

Wang, S.N.; Liang, L; et.al.;" Solder Joint Reliability under Electromigration and Thermal-

Wever, H.: "Stofftransport im metallischen Festkörper", Angewandte Chemie, Nr.7, 1963,

Weide-Zaage, K.; Zhao, J.; et.al: "Determination of Migration Effects in Cu-Via Structures with Respect to Process Induced Stress", Microelectronic Reliability 48, 2008, pp. 1393–

Mechanical Load ", Elec. Comp. & Tec. Conf, 2007 pp1074-1083.

Papers, 2009, pp. 70-71.

June 2011, pp. 1444-1451.

IEEE/IEDM, 2011

Physique, 1991.

5,198.

1397.

pp.309-352.

Engineering, 1963/1989.

Technology Conf. 2011, pp. 332-339.

stacked ICs", EuroSimE 2012, p. 1-4.

Journal of applied Physics 104, 023521, 2008

2010, p. 1-4.

Symposium, IEEE/IRPS 2008, pp. 683-684.

Packaging Conference (IWLPC), October 2010, pp. 137-145.


Lo, W.-C. ; Chen, Y.-H. ; et.al: TSV and 3D Wafer Bonding Technologies For Advanced Stacking System and Application at ITRI, Symp. on VLSI Tech. Digest of Technical Papers, 2009, pp. 70-71.

396 Finite Element Analysis – New Trends and Developments

pp.1-8.

210.

978-3-211-82110-7.

Tokyo 2010.

1999, pp. 281–304.

2006, pp. 671-672.

700–709.

pp. 811-824.

Hou, Y; Tan, C.M.: "Comparison of stress-induced voiding phenomena in copper line-via structures with different dielectric materials", Semiconductor Sci. Technology, 2009,

Huang, J.S.; Oates, T.S.: "Monte-Carlo simulation of electromigration failure distributions of submicron contacts and vias: a new extrapolation methodology for reliability estimate", Proceedings of the IEEE Interconnect Technology Conference 2000, pp. 208-

Huffman, A.; Lueck, M.; et.al: "Effects of Assembly Process Parameters on the Structure and Thermal Stability of Sn-Capped Cu Bump Bonds", ECTC Conf. 2007, pp. 1589-1596. Jacoboni, C; Lugli, P.: "The Monte Carlo Method for Semiconductor Device Simulation", Computional Microelectronics", Edited by S. Selberherr, Vienna, Springer 1989, ISBN:

Jaffe, D.; Shewman, P.G.: "Termal Diffusion of Substitutional Impurities in Copper, Gold

Joo, Y.C.; Yang, T.Y.; et.al.: "Driving forces of mass transport in phase change materials and

Kaur, I.; Mishin,Y.; Gust, W.;: "Fundamentals of Grain and Interphase Boundary and

Kitada, H.; Maeda, N.; et.al.: Development of Low Temperature Dielectrics down to 150°C for Multiple TSVs Structure with Wafer-on Wafer (WOW) Technology, IEEE, CMPT

Labie, R.; Limaye, P.; et.al: "Reliability testing of Cu-Sn intermetallic micro-bump interconnections for 3D-device stacking", IEEE/Electronics-System-Integration-

Leduc, P.; Assous, M.; et.al: "First integration of Cu TSV using die-to-wafer direct bonding and planarization", IEEE Int. Conf. on 3D System Integration, 3DIC, 2009, pp. 1-5. Lee, K.-D.; Park, Y.-J.; Kim, T.: "Via Processing Effects on Electromigration in 65nm

Li, Z; Zhao, H.; Gao, H.: "A Numerical Study of Electro-migration Voiding by Evolving Level Set Functions on a fixed Cartesian Grid", Journal of Computational Physics 152, 2,

Lin, M.H.(a); Chang, K.P.; et.al: 'Effects of Width Scaling, Length Scaling, and Layout Variation on Electromigration in Dual Damascene Copper Interconnects', IEEE RPSP,

Lin, M.H.(b); Lin, M.T.; et.al: "Copper Interconnect Electromigration Behavior in Various Structures and Precise Bimodal Fitting", Jap. J. Appl. Phys. Vol. 45, No. 2A, 2006, pp.

Liu, Y.; Irving, S.; et.al: "3-D Modelling of Elec-tromigration Combined with Thermal-Mechanical Effect for IC Device Package", 8th Int. Conf. EuroSimE 2007, pp.1-13. Liu, Y.; Liang, L.; et al.: "3D Modeling of electromigration combined with thermalmechanical effect for IC device and package", Microelectronics Reliability, Vol.48 (2008),

Technology, IEEE 06CH37728 44th Annual IRPS, 2006, pp.:103-106.

and Silver", Acta Metallurgica, Vol.12, 1964, pp.515-527.

their effect on device failures", EPCOS 2011.

Diffusion", John Wiley & Sons Ltd., 1995.

Conference (ESTC), Berlin, September 2010.

	- Weide-Zaage, K.; Ciptokusumo, J.; Aubel, O.: "Influence of the Activation Energy of the Different Migration Effects on Failure locations in Metallizations", AIP Conf. Proc. -- November 24, 2010 -- Volume 1300, pp. 85-90.
	- Wilkenson, D.S.: "Mass Transport in Solids and Fluids", Cambridge University Press, ISBN 0521624096, 2000.

0521624096, 2000.

November 24, 2010 -- Volume 1300, pp. 85-90.

Weide-Zaage, K.; Ciptokusumo, J.; Aubel, O.: "Influence of the Activation Energy of the Different Migration Effects on Failure locations in Metallizations", AIP Conf. Proc. --

Wilkenson, D.S.: "Mass Transport in Solids and Fluids", Cambridge University Press, ISBN

## *Edited by Farzad Ebrahimi*

In the past few decades, the Finite Element Method (FEM) has been developed into a key indispensable technology in the modeling and simulation of various engineering systems. The present book reports on the state of the art research and development findings on this very broad matter through original and innovative research studies exhibiting various investigation directions of FEM in electrical, civil, materials and biomedical engineering. This book is a result of contributions of experts from international scientific community working in different aspects of FEM. The text is addressed not only to researchers, but also to professional engineers, students and other experts in a variety of disciplines, both academic and industrial seeking to gain a better understanding of what has been done in the field recently, and what kind of open problems are in this area.

Photo by FotoMak / iStock

Finite Element Analysis - New Trends and Developments

Finite Element Analysis

New Trends and Developments