Tips & Tricks
5. Simulating Stress Effects in Back-End-of-Line Structures

5.1 Effective Material Approach
5.2 Converting Elastic Material Parameters
5.3 Avoiding Zero Plastic Hardening Modulus
5.4 Taking Advantage of Symmetry
5.5 Mechanical Coupling Boundary Conditions
5.6 Structure Generation for Cylinders and Spheres in BEOL
5.7 Structure Generation for Unconventional BEOL Process Steps
5.8 Defining Contact and Crack Locations in the Unified Coordinate System
5.9 Linearizing Dielectric Materials at BEOL Process Temperatures
5.10 Time-Step Control for BEOL Stress Simulation With Only Linear Models
5.11 Time-Step Control for BEOL Stress Simulation With Nonlinear Models
5.12 Modeling Deformation With Moving Mesh
5.13 Modeling Crack Opening on Top Surface
5.14 J-Integral Calculations With Multiple Contours
5.15 General Guidelines for Submodel Technique
5.16 Meshing Parameter Recommendations for Large Structures With Curved Surfaces
5.17 Guidelines for Cohesive Zone Modeling
5.18 Defining Surface Contacts to Apply Mechanical Boundary Conditions

Objectives

5.1 Effective Material Approach

When analyzing stress for a large structure, for example, the thermal mismatch stress between the silicon chip and package substrate, it is not practical to consider fine features such as metal wires in the BEOL layers. However, the contribution of the stressed BEOL stack to the stress equilibrium of the chip-package structure cannot be simply ignored.

One approach to effectively consider the BEOL contribution is to use the concept of "smear" material. With this approach, the BEOL layers are treated as a single layer of "smear" material. The "smear" material properties are the volume average of the metal wires and dielectrics. Typically, three elastic model parameters must be defined: Young's modulus, Poisson ratio, and coefficient of thermal expansion (CTE). For Young's modulus, the "smear" material parameter is calculated as:

Esmear = Emetal x Volume percentage of metal + Edielectrics x Volume percentage of dielectrics

Similar expressions can be written for the Poisson ratio and CTE. These new material properties are introduced in the simulation flow with the following statements:

mater add name= Smear new.like= Nitride
pdbSet Smear Mechanics BulkModulus 33.0e10
pdbSet Smear Mechanics ShearModulus 28.2e10
pdbSet Smear Mechanics ThExpCoeff 15.3e-6

5.2 Converting Elastic Material Parameters

Sentaurus Interconnect uses the bulk modulus and the shear modulus to define material elastic properties, while many times, available BEOL material properties are given in terms of Young's modulus and the Poisson ratio. The following statements convert user data in terms of Young's modulus and the Poisson ratio to the bulk modulus and the shear modulus required in the simulations:

pdbSet Copper Mechanics BulkModulus  [Enu2K 111.5e10 0.343]
pdbSet Copper Mechanics ShearModulus [Enu2G 111.5e10 0.343]

5.3 Avoiding Zero Plastic Hardening Modulus

When modeling plasticity with bilinear hardening laws, the onset of plastic deformation is controlled by the material yield stress, and the material hardening is measured by the hardening modulus. For numeric reasons, it is recommended to avoid using the zero plastic hardening modulus.

Even for rigid plastic materials, a small nonzero value of the hardening modulus can often improve numeric convergence. The following statements can adjust the hardening modulus:

pdbSetDouble <material> Mechanics Hardening.Modulus.Isotropic <n>
pdbSetDouble <material> Mechanics Hardening.Modulus.Kinematic <n>

The parameter Hardening.Modulus.Isotropic corresponds to the isotropic hardening modulus, and the parameter Hardening.Modulus.Kinematic corresponds to the kinematic hardening modulus.

5.4 Taking Advantage of Symmetry

Symmetry is commonly observed in BEOL structures. Using symmetric boundary conditions can substantially reduce the computational task. For example, stress analysis for a through-silicon via (TSV) array can be performed using a quarter of a TSV with symmetric boundary conditions applied. The mechanical symmetric boundary conditions can be specified as:

pdbSet Mechanics Left  BoundaryCondition Dirichlet
pdbSet Mechanics Right BoundaryCondition Dirichlet
pdbSet Mechanics Front BoundaryCondition Dirichlet
pdbSet Mechanics Back  BoundaryCondition Dirichlet

Similarly, the symmetric boundary conditions can be applied with stressdata commands:

stressdata bc.location=Left  bc.value= { dy=0.0 }
stressdata bc.location=Right bc.value= { dy=0.0 }
stressdata bc.location=Front bc.value= { dz=0.0 }
stressdata bc.location=Back  bc.value= { dz=0.0 }

After the simulation, the structure with solution fields can be reflected multiple times to generate a large structure for visualization purposes.

5.5 Mechanical Coupling Boundary Conditions

When simulating material thermal mismatch of a unit cell inside a large structure with finite layer thicknesses, the deformation of all layers often needs to be modeled explicitly. To approximate the constraints from surrounding materials, mechanical coupling boundary conditions on the unit cell can be applied as:

pdbSet Mechanics <Left | Right | Back | Front> <Cx | Cy | Cz> 1

Here, Cx, Cy, or Cz stands for the coupling of the displacement component in the x-, y-, or z-direction.

All other boundary conditions defined with the pdbSet method become invalid with the application of the coupling boundary condition. These boundary conditions must be redefined using the stressdata command.

5.6 Structure Generation for Cylinders and Spheres in BEOL

Cylinder and sphere objects are common in BEOL structures, for example, interconnect vias, TSVs, microbumps, and solder bumps. These objects cannot easily be generated by deposition and etching with standard masks.

There are several practical approaches to generate these objects. One approach is to generate these objects in Sentaurus Structure Editor and then insert them in the simulation. The following example shows the statements to generate a solder ball:

sde external { (sdegeo:create-sphere (position $x $y $z) 
  $radius "material" "region") }
polyhedron name= sphere1 external.sde
sde off
insert polyhedron= sphere1 replace.materials= { replace_material }

Another approach for generating cylinders is to create circular masks either in IC WorkBench Edit/View Plus (ICWBEV Plus) or in simulation input, and then perform etching and deposition with circular marks. When applying this approach, you must be aware that the number of points defined around a circle will affect the final mesh. Usually, 36 points provide good enough resolution. Too many points can lead to an unnecessarily dense mesh. You also can use the shape library for some shapes, for example, PolygonWaferMask.

5.7 Structure Generation for Unconventional BEOL Process Steps

BEOL structure generation can involve unconventional process steps such as die stacking, bumping, and underfilling. These process steps result in structures that cannot be generated by conventional deposition and etching operations. Sometimes, the structure must be flipped and the location must be translated.

A practical approach for generating such structures is to build structure components individually and then to assemble them. This assembly process can be achieved using statements such as polyhedron tdr=... to insert one structure component into the target structure, transform flip to flip a component before assembly, and transform translate ... to align components in the required locations.

5.8 Defining Contact and Crack Locations in the Unified Coordinate System

BEOL stress simulation and failure analysis often require the definition of contact and crack locations. Although the unified coordinate system (UCS) is recommended, the DF–ISE coordinate system also is supported for structure generation in Sentaurus Structure Editor. The coexistence of the UCS and the DF–ISE coordinate system can lead to confusion when defining contact and crack locations.

Be aware that all geometry location definitions must be specified in UCS coordinates. For example:

contact box xlo=-0.01 ylo=-0.46 xhi=0.1 yhi=-0.16 name=source Aluminum
crack name=Si_Ox mat.1=Silicon mat.2=Oxide CZM=none tolerance=1.0e-6 \ 
  x.min=-0.1 x.max=-0.05 y.min=-0.1 y.max=-0.05 z.min=0.0 z.max=0.4 \
  normal= { 0.0 0.0 1.0 } angle.tolerance=1.0e-3

The specified coordinates are all in UCS. You must manually convert the coordinates if the original structures are in DF–ISE coordinates.

5.9 Linearizing Dielectric Materials at BEOL Process Temperatures

Dielectric materials often exhibit viscoelastic behavior at elevated temperatures but elastic behavior at BEOL process temperatures. By default, TCAD tools use viscoelastic models for dielectric materials. This default setting often leads to unnecessary computation during BEOL stress simulations. It is recommended to linearize dielectric-material mechanical properties at BEOL temperatures with the following statements:

pdbSet <material> Mechanics Vcrit0 1e-30
pdbSet <material> Mechanics VcritW 0.0
pdbSet <material> Mechanics Viscosity0 1e40
pdbSet <material> Mechanics ViscosityW 0.0

5.10 Time-Step Control for BEOL Stress Simulation With Only Linear Models

For elastic materials in BEOL structures, the constitutive relations are linear. When performing a thermal ramp from the initial temperature to the final temperature, the solve command must be constructed in such a way that only one solve step is performed. For example:

solve temperature=170.0<C> t.final=25.0<C> delNT= 145.0 \
  time=0.1<s> init= 0.1<s> maxstep= 0.1<s>

In this example, the total time, the initial time step, and the maximum time step are all set to 0.1 s. The maximum temperature step during the temperature ramp-down from 170°C to 25°C is defined as 145°C. To avoid large ramp rate–induced numeric issues, the time step should not be too small.

5.11 Time-Step Control for BEOL Stress Simulation With Nonlinear Models

For some materials in BEOL structures such as copper, plastic deformation can occur when the material stress level exceeds the yield strength. For some other materials, creep and viscoplastic deformation can occur. In addition, cohesive surface separation is a nonlinear process. The nonlinear constitutive relation requires small time steps to avoid numeric convergence problems.

It is recommended to limit both the initial time step and the maximum time step for each loading step. The following command sets the initial time step to 1% of the thermal loading and the maximum time step to 10% of the loading:

solve temperature= 170.0<C> t.final= 25.0<C> time= 0.1<s> \
  init= 0.001<s> maxstep= 0.01<s>

In this example, the initial loading is set to (170 – 25) x 0.001°C, and the maximum time step is limited to (170 – 25) x 0.01°C.

5.12 Modeling Deformation With Moving Mesh

When a moving mesh is used to resolve structural deformation and displacement under either a point force or a contact loading, it is important to ensure that only small deformations and small displacements occur under each loading/moving step.

For point-force loading, the following parameter must be set to update the point-force location during deformation:

pdbSetBoolean Mechanics Point.Force.Snap 1

For contact loading, the new contact location after mesh moving must be within the original contact box.

5.13 Modeling Crack Opening on Top Surface

When modeling a crack opening on the top surface, the top gas region must be removed in the simulation to allow proper crack surface separation and crack opening visualization. This requirement might be removed in a future release.

5.14 J-Integral Calculations With Multiple Contours

When calculating J-integral at the crack tip, the mesh must be refined around the crack tip front. The J-integral calculation should always be performed with multiple contours. The following example selects nine contours for the calculation:

j_integral crack=edge x=2.0 y=3.0 z=1.0 number=9

The calculated J-integral values should be examined, and the stable contour value away from the crack tip front should be selected.

5.15 General Guidelines for Submodel Technique

When applying the submodel technique to obtain accurate resolution in some local regions of interest, several guidelines need to be observed:

5.16 Meshing Parameter Recommendations for Large Structures With Curved Surfaces

The default meshing parameters, which are selected to resolve geometry details in small structures with flat surfaces, can lead to highly anisotropic meshes with thin layers of elements at material interfaces and overrefinements on curved surfaces. To generate isotropic-like meshes with regular equilateral elements and to prevent overrefinement on curved surfaces, it is recommended to use meshing parameters that are more suitable for stress analysis in large structures. These meshing parameters can be invoked by the following command:

SetMechanicsMeshMode

After this command, regular mesh-generation commands, such as grid remesh, are needed to start mesh-generation operations. Since geometries can be very different in real applications, some manual adjustments of these parameters might still be needed to achieve optimal mesh qualities for stress analysis. The important meshing parameters include:

pdbSet Grid SnMesh Apply.Brep.DelPSC 1               ;# default 0
pdbSet Grid SnMesh Apply.Brep.DelPSC.Resolution 0.1  ;# default 0.01 [um] 
pdbSet Grid SnMesh Apply.Brep.DelPSC.Accuracy 0.005  ;# default 0.0001 [um]
pdbSet Grid SnMesh ImprintAccuracy 0.1               ;# default 1e-5 [um]
pdbSet Grid SnMesh SliverDistance 0.1                ;# default 1e-3 [um]

where:

5.17 Guidelines for Cohesive Zone Modeling

Cohesive zone modeling requires defining (1) a cohesive zone model (CZM) law that governs the traction and separation behavior on cohesive surfaces, and (2) a cohesive path where cohesive surfaces can separate to form a crack.

The following commands define a CZM law named law1 and a cohesive path named crk01:

CZM name=law1 law=exponential \
   normal.stress.max=1.0e9 shear.stress.max=1.0e9 \
   normal.toughness=1e5 shear.toughness=1e5 mat.1=Mat1 mat.2=Mat2

crack name=crk01 CZM=law1 \
   x.min=0 y.min=0.0 x.max=40.0 y.max=3.0 z.min=0.0 z.max=3.0 \
   tolerance=1.0e-3 mat.1=Mat1 mat.2=Mat2 angle.tolerance=1.0e-4 \
   normal= "0 1 0" \
   CZM.box= {xmin= 5.0 xmax= 40.0 ymin= 0.0 ymax= 3.0 zmin= 0.0 zmax= 3.0}

In the CZM law definition, the traction and separation relation is chosen to be an exponential function. The maximum material strength in the normal and shear directions are both 109 dyn/cm2. The material work of separation (that is, material fracture toughness) in the normal and shear directions are both 105 dyn/cm. The CZM law applies to a cohesive path between the first material Mat1 and the second material Mat2.

In the crack and cohesive path definition, the entire cracking region contains an initial crack and a cohesive path ahead of the initial crack. The cracking region is named crk01 and the CZM law to be applied is law1.

The larger cube that encloses the initial crack and the cohesive path extends from the lower-front corner (0, 0, 0) to the upper-back corner (40 µm, 3 µm, 3 µm) with a snapping distance of 10–3 µm.

The smaller cube that encloses only the cohesive path extends from the lower-front corner (5 µm, 0 µm, 0 µm) to the upper-back corner (40 µm, 3 µm, 3 µm).

The initial crack and the cohesive path are between the first material Mat1 and the second material Mat2 with a surface normal (0 1 0) and a tolerance of 10–4. All of the coordinates must be defined in the unified coordinate system (UCS).

When cohesive surfaces separate along the cohesive path, the cohesive traction first increases, reaches the material strength, then decreases to zero after the total separation. To capture this entire separation event, the mesh must be relatively uniform around the cohesive path, and the local element size must be comparable to the ratio of the material fracture toughness divided by the material strength. To allow 3D region separation in the mesh and vertex separation at a point on the cohesive path, the following settings must be applied:

pdbSet Mechanics Model CZM
pdbSet Grid No3DMerge 1
pdbSet Mechanics Merge.Vertices 0

5.18 Defining Surface Contacts to Apply Mechanical Boundary Conditions

When defining surface contacts to apply mechanical boundary conditions, it is important to note the following:

main menu    |   module menu    |   << previous section    |   next section >>