Microelectronic devices are complex systems fabricated through the manipulation and integration of many different materials. In order to be able to understand and improve such technology a multidisciplinary approach is necessary.
As explained in
The finite element method originated from the needs to solve mechanically related problems. In this work a FEM tool is used to calculate the resulting stress/strain in microelectronic components under several loading conditions and therefore an overview of FEM is presented.
The continuum mechanics approach ignores the discrete nature of matter, considering materials as uniformly and continuously distributed throughout regions of space. This simplification leads to the definition of material properties such as conductivity, permittivity, yield strength, coefficient of thermal expansion and simulation quantities such as density, displacement, velocity, etc., as continuous functions of position.
Two main types of equations are used in continuum mechanics [35]:
In this section the kinematic and the stress equations are described. Subsequently, a brief description of the constitutive equations concerning elasticity is provided.
This section, which is based on [35, 36], deals with the kinematic relations of a continuous medium. To define the basic deformation and strain concepts we must analyze material kinematics, which refers to the mathematics of motion without considering the masses of objects or the forces which may have caused the motion. Therefore, a kinematic analysis permits to study motion and displacement within a material, without regard to the associated stresses. This permits to analyze the system essentially as a geometric problem. Kinematics is widely applied to solve mechanical problems [36].
In order to describe the motion of a body the definition of a reference coordinate system is
necessary (cf.
At time
 (2.1) 
which corresponds to
 (2.2) 
In this equation
The deformation of the body can be described by only considering two configurations of the body, an initial and a final configuration where the motion of the body can be regarded as a singleparameter sequence of deformations. In (2.1) the mapping has the unique inverse
 (2.3) 
or
 (2.4) 
From (2.3) the existence of the Jacobian determinant is ensured
 (2.5) 
which is positive for all points
The coordinates
Lagrangian or Eulerian coordinates can be applied to describe the displacement vector
 (2.6) 
or by using (2.3) as a reference
 (2.7) 
The material displacement gradient is obtained by a partial differentiation of the displacement vector (2.6) with respect to the Lagrangian coordinate
 (2.8) 
where
 (2.9) 
In a similar way the material deformation gradient is defined by
 (2.10) 
and the spatial deformation gradient is defined by
 (2.11) 
From (2.1) and (2.10) the distance differential
 (2.12) 
where
 (2.13) 
An important value is the difference
 (2.14) 
In the equation above, the secondorder tensor with component
 (2.15) 
where the tensor
 (2.16) 
is called right CauchyGreen tensor.
The Lagrange strain tensor (2.15), described in terms of displacement (2.6), takes the form
 (2.17) 
if
In addition the strain measure can be formulated by employing the spatial coordinate
 (2.18) 
where
 (2.19) 
If each displacement gradient component
 (2.20) 
Implementing the same methodology for
 (2.21) 
which is know as the classical strain tensor.
The Eulerian infinitesimal strain tensor is an approximation which is only valid for small changes and can be expressed as [37]
 (2.22) 
or using the engineering notation
 (2.23) 
The components
 (2.24) 
however we will stick to the notation from (2.23) in this work. The formulation in (2.21) forms a
system of six equations for the three displacement components. The strain components must
satisfy the compatibility conditions. Because the displacement field
 (2.25) 
When a body is deformed internal forces known as stress arise inside it. The stress generates to restore the undeformed state of the body. At mechanical equilibrium the sum of all acting and internal forces must equal to zero. This section, based on [35, 38], provides the relationships between force and stress.
In order to study how the body exterior interacts with the body interior we imagine a closed
surface
We consider a unit vector
 (2.26) 
The components of the stress vector are normal stress and shear stress. Normal stress is defined
as
By considering three specific sections corresponding to the coordinate system
 (2.27) 
where
 (2.28) 
The stress component
The stress components
The stress tensor (also called Cauchy stress tensor)
Considering a 3D system, the stress tensor is described through nine stress components [39]
 (2.29) 
where
If the stress tensor is symmetric the number of stress components reduces to six because
 (2.30) 
By applying Newton’s first law to a solid subjected to body force components
represented by
 (2.31) 
where
 (2.32) 
and then equation (2.31) becomes
 (2.33) 
where
 (2.34) 
The Cauchy stress tensor
In engineering the concept of the scalar Von Mises effective Stress is commonly used to predict yielding of materials under a loading conditions as well as other physical phenomena [36]
 (2.35) 
The equations described above are applicable for all continuous materials. However additional equations, called constitutive equations, are necessary to describe the mechanical behavior of any particular material.
Four theories for material behavior can be distinguished [35]:
Creep behavior exists in viscoelasticity and viscoplasticity theories. Creep is the phenomenon where a solid material moves slowly or deforms permanently under stresses.
Materials are additionally described as brittle or ductile [37]. In brittle materials rupture occurs without any noticeable prior deformation, behaving elastically before rupture. For ductile material the rupture arrives after a plastic deformation.
The different material behaviors listed above can be analyzed using stressstrain curves. For
each material a unique relationship between stress and strain can be found. Stressstrain curves
are obtained by recording the strain at distinct intervals of applied stress.
The stressstrain curve describes the elastic modulus
The behavior of materials and their properties are correlated to the arrangement of atoms in
the material body. Singlecrystal is a body where atoms are ordered in a pattern which repeats
itself; therefore, only a single crystalline phase is present. Polycrystalline materials are composed
of more than one crystal (grains) and often combine different crystalline phases (cf.
At the atomic scale elastic and plastic deformations are significantly different (cf.
It is therefore clear that the plastic mechanism is influenced by the microstructural features such as slip plane, lattice imperfections, alloying, and impurity elements. The plasticity is also influenced by the movement, multiplication, and interaction between dislocations [43]. Therefore, the microstructures of materials define properties such as elasticity, plasticity, and creep.
An elasticplastic material is usually studied assuming that the strains
 (2.36) 
We take a closer look on the elastic component of a linear stressstrain relationship in the next subsection. For the plastic component a different constitutive law, based on plastic strain increments or total plastic strains, are in use [35, 37].
For a linear elastic material, assuming small (infinitesimal) strain (
 (2.37) 
where
 (2.38) 
which reduces the number of independent components from 81 to 36. Equation (2.37) can be inverted and the elasticity law becomes
 (2.39) 
where
For an isotropic material,
 (2.40) 
where
 (2.41) 
where
 (2.42) 
Many 3D solid mechanics problems can be treated as 2D problems. Plane strain and plane stress
are two idealized stress configurations, which simplify problem analysis [37]. The
two approximations are described by using the following engineering notation with
coordinates
The plane strain is defined to be a state which assumes that
 (2.43) 
the kinematic relations to
 (2.44) 
and compatibility conditions become
 (2.45) 
This assumption permits to reduce the stressstrain relation to
 (2.46) 
On the other hand, in the plane stress approximation,
 (2.47) 
In microelectronic devices, thin films can have residual stresses, which arise due to two main sources: intrinsic and thermal stress [45]. High values of residual stress can affect the quality and performance of thin films and subsequently of the entire devices. Therefore, the investigation of the origin of the residual stress is a fundamental aspect to examine.
Intrinsic stress emerges during the deposition processes (e.g. sputtering, spraying, painting, spin
coating, CVD, PVD, and electrodeposition). The magnitude of the intrinsic stress depends on
the process, deposition temperature, and deposited material. During the deposition process, the
microstructure evolves from single atom islands to a uniform film, influencing the magnitude of
intrinsic stress [45]. In
In absence of mechanical loads a body can deform due to temperature changes. This deformation is called thermal deformation, which gives rise to so called thermal strain inside the body [39]. For small temperature changes the thermal strain is linearly proportional to the temperature change and it is described by the Coefficient of Thermal Expansion (CTE). CTE is a material property which describes how a material expands upon heating or contracts upon cooling. For changes in volume the volumetric CTE can be expressed as
 (2.48) 
where
 (2.49) 
where
 (2.50) 
where the change from
Thermal stress is the stress generated in a constrained structure due to a thermal expansion or contraction. Two different set of constraints are possible: external constraint and internal constraint. External constraint can be every external object which obstruct the free expansion or contraction of the system, when temperature changes occur. Internal constraint usually appear in composite systems, formed by different materials, or when a nonuniform temperature distribution is present. Due to the different CTEs every material behaves differently under thermal variation. The thermal strain generated in a material can affect and then generate thermal stress, thereby impacting the system.
The relationships between stress and strain for an isotropic material, including the thermal strain are given as
 (2.51) 
The principle of virtual work is applied to study the mechanics of deformable bodies. The term virtual refers to the implementation of infinitesimal displacements, called virtual displacement, which do not represent a real displacement acting on the system. By considering the equilibrium equations of a system, where all the forces acting on a system in static equilibrium must be equal zero, virtual displacement can be applied and the resulting deformation of the body can be obtained [36, 37].
We consider a body of volume
 (2.52) 
To derive the concept of the virtual work it is necessary to recall the equilibrium condition and the classic strain tensor:
 (2.53) 
 (2.54) 
The equilibrium condition is multiplied by the virtual displacement
 (2.55) 
By using the product rule
 (2.56) 
and the symmetry of (2.54) equation (2.55) can be rewritten as
 (2.57) 
which, when is integrated over the volume of
 (2.58) 
Gauss’s divergence theorem is applied at the first volume integral, becoming a surface integral
 (2.59) 
By using the boundary condition the equation above can be rewritten in a simplified form
 (2.60) 
where the integrals on the left represent the virtual work done by external forces. On the right hand side the integral represents the virtual work done by internal forces. The system is in equilibrium if the entire virtual work of the imposed virtual displacements vanishes, therefore (2.60) is equivalent to (2.34).
Fracture is a degradation process which leads to the generation of new surfaces. Fracture develops through cracking, which is the development of displacement discontinuity surfaces within the solid. Previously a continuummechanical analysis was presented, however a microscopic approach is necessary to understand the role of the microstructure in fracture mechanics [37]. At the microscopic level, the fracture process consists of breaking the bonds between elements (atoms, ions, molecules, etc.). The bonding force between two elements is described by
 (2.61) 
where
 (2.62) 
where
In a polycrystalline material the separation of atomic planes (microcracks) is often combined with pronounced microplastic processes (dislocations). The dislocation pileup at obstacles is an important phenomenon during microcrack formation.
It can cause a high stress concentrations leading to bond breaking along favored lattices places:
These processes describe fracture in brittle materials. In plastic deformation the load is not enough to break the atomic bonds but it is large enough to produce permanent deformation by shearing atomic planes (dislocations); therefore, the failure is the plastic regime.
Cracks are classified from different phenomenological viewpoints. A stationary crack is a
crack which does not change in length. If the crack starts to propagate and it becomes
nonstationary a specified load or deformation is called crack initiation. Crack propagation
can further be distinguished as stable, when the increase in crack length requires an
increases external load, or unstable, when the crack advances without any increase in the
external load. In
In
A physical problem can often be described by the following general differential equation
 (2.63) 
where
For a welldefined problem set, boundary conditions are required. There are three types of
boundary conditions: The first is called Dirichlet boundary condition where a domain
 (2.64) 
where
 (2.65) 
where
 (2.66) 
where
In this section we describe the general method of the weighted residual and the derivation of the solution using different approaches.
Onedimensional physical problems in a domain
 (2.67) 
In (2.67)
 (2.68) 
When an exact analytical solution cannot be found, an approximated solution of the problem is obtained by employing numerical methods. The goal is to calculate a useful approximation of the exact solution
 (2.69) 
In the case of the weighted residual method [47, 49] an approximation of the solution of the following form is chosen
 (2.70) 
where
A local error know as residual
 (2.71) 
Performing an integration over the entire domain
 (2.72) 
which represents an inner product. By choosing the appropriate
 (2.73) 
where
 (2.74) 
where the Dirac delta function is:
 (2.75) 
By applying the test function (2.74) in the inner product (2.72) and considering the property of the Dirac delta function,
 (2.76) 
 (2.77) 
and choosing
 (2.78) 
Because
 (2.79) 
which can be written in matrix form
 (2.80) 
The results can be obtained by inverting the matrix:
 (2.81) 
By obtaining the parameters
 (2.82) 
where the indicator function is defined by
 (2.83) 
Therefore the test function is defined as
 (2.84) 
The integral of the error becomes zero over the different domains:
 (2.85) 
 (2.86) 
The parameters
 (2.87) 
where
In the method of least squares the weighted integral is written as
 (2.88) 
where, employing the chain rule, the following form can be derived:
 (2.89) 
This finally leads to the minimization of the square residual:
 (2.90) 
In detail, the employed test function is
 (2.91) 
thus this equation can be inserted in (2.88)
 (2.92) 
The equation can then be rewritten as
 (2.93) 
and
 (2.94) 
From this linear system it is possible to obtain the
 (2.95) 
Equation (2.95) can be plugged into (2.72)
 (2.96) 
The parameters
 (2.97) 
which can be written as system of linear equations:
 (2.98) 
from which is possible to obtain the
 (2.99) 
where
The definition (2.99) can be plugged into (2.72), obtaining
 (2.100) 
By inserting the approximation function (2.70) in the above equation we obtain
 (2.101) 
This integral can be separated
 (2.102) 
and due to the linearity of the differential operator the summation can be extracted
 (2.103) 
The weighting functions can be arbitrary chosen and therefore we consider only one
 (2.104) 
Using this method for all
 (2.105) 
which can be expressed in matrix form by
 (2.106) 
with the following definitions:
 (2.107) 
where
A solution to (2.63), which is called the strong formulation, is very difficult or impossible for typical engineering problems. It requires continuity on the dependent field variables and that the approximation function be differentiable up to the order of the differential equations. For such problems, the weak formulation can be used, which mathematically corresponds to the strong formulation. Usually the weak formulation has an integral form and it requires a weaker continuity on the field variables. It also allows for a reduction in the necessity for the differentiability of the approximation function by one order and permits a native incorporation of the Neumann boundary condition [50]. Here the transition from the strong formulation to the weak formulation is described [48].
In a 3D system a differential equation in domain
 (2.108) 
where
 (2.109) 
where
The weak formulation is obtained by multiplying (2.108) by (2.73) and integrating over the
entire domain
 (2.110) 
Subsequently, the
 (2.111) 
Therefore, the weak formulation for the initial boundary problem becomes:
 (2.112) 
The Neumann boundary conditions are incorporated in (2.112), while the Dirichlet boundary
condition
The approximated solution can be obtained by plugging the equations (2.73) and (2.70) into
(2.112) and by choosing
 (2.113) 
where for
 (2.114) 
whit values
 (2.115) 
The
FEM is based on splitting complex domains or subdomains into geometrically simple shapes called finite elements. These elements, forming a mesh, are defined in 1D by intervals, in 2D usually by triangles or quadrilaterals, and in 3D usually by tetrahedrals or hexahedras. The size of these elements is a crucial parameter, as a finer mesh results in a more accurate approximation of the solution, while a courser mesh reduces the computational time.
All the single elements
In each subdomain, due to the meshing, a number of nodal points are defined. These nodal
points are the supporting points of the basis functions
As an example, we consider a 1D domain divided into elements of small intervals
 (2.116) 
 (2.117) 
 (2.118) 
In order to simplify the example we use an equidistant discretization

where the mesh size
 (2.119) 
The approximation function
 (2.120) 
then, on the nodal points, the chosen basis functions fulfill the conditions
 (2.121) 
and the approximate function
 (2.122) 
This allows the calculation of the unknown parameters
The chosen polynomial functions have to be continuous over the entire domain.
A typical application for FEM are solid mechanics problems. Here, a general derivation and implementation for a linear elastic solid is presented [36]. The description follows a general framework used to find solutions for problems in solid mechanics. Different approaches to construct the weak form, as well as expressions for the strain and the constitutive law, can be found in [36, 49].
The system depicted in
The coefficients
In
If
 (2.123) 
for all virtual displacement fields and virtual strains, the equation of stress equilibrium
To better explain the FEM procedure the first integral in (2.123) can be rewritten. In fact,
by employing the symmetry of the stress tensor
 (2.124) 
Furthermore for linear elastic solids we recall
 (2.125) 
The expression (2.125) can be plugged into the first term of (2.123) and the virtual work can be rewritten as
 (2.126) 
The displacement field
The strain components can be obtained from
The integral form of (2.126) can be approximated by discretizing the displacement field,
which can be calculated as a set of
The FEM procedure specifies the displacement field at an arbitrary point within the solid by interpolating between nodal values. The interpolation function can be defined as
 (2.127) 
The same procedure can be used for the virtual displacement, which represents the test function
 (2.128) 
Now (2.127) and (2.128) are substituted in (2.126). By choosing
 (2.129) 
The system of equations in (2.129) can be rewritten in matrix form as
 (2.130) 
where the notations indicate
 (2.131) 
Therefore, the