# Modeling the Boron-Doping Silicon Beam by a Multilayer Model

Authors: Huijun Yu; Wu Zhou; Bei Peng; Xiaoping He; Xiaohong Hao; Zhi Zeng
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The boron-doping silicon beam commonly used in microdevices exhibits a nonuniform material property along its thickness or width because of the gradient of boron concentration induced by diffusion process. The constant of rigidity, one of the most important parameters of microbeam, needs to be accurately calculated and designed in the development of high precise sensors and actuators. Current design methods, mainly depending on the analytical solutions derived under the assumption of a uniform material property or some commercial software for a varied property, are not adequate and time consuming to calculate the constant of rigidity of boron-doping silicon beam. A multilayer model is proposed in this paper to replace the continuous solid model by dividing the beam into separated layers glued together. The finite element lamination method is utilized to acquire the equivalent Young modulus and moment of inertia of cross section of multilayer model. The equivalent values are calculated from double-layer structures to multilayer ones based on the small deformation theory and the material mechanics theory. The proposed method provides an effective method to design the stiffness or frequency of microdevice and its results are validated by COMSOL simulation.


---

## Body

## 1. Introduction

Boron diffusion is the process commonly used in silicon based MEMS for acquiring a certain electrical property in microelectronics [1] or a process requirement in bulk silicon micromachining [2, 3]. It is one of the most important steps in bulk silicon dissolved wafer process (BSDWP) which is extensively used in fabrication of microaccelerometers [4, 5], microgyroscopes [6, 7], microswitches [8, 9], and microgears [10]. The microstructures fabricated by BSDWP, however, exhibit a varied property from top surface to bottom surface. When the movable or deformable parts are involved in devices, the boron-doping silicon is always chosen as the structural material for the flexible components, like beams [11], flexure hinges [12, 13], and so forth. These components play the role of springs in microsensors to connect the movable parts to the substrate. Their stiffness or rigidity, therefore, has significant importance to the design of movable structures in MEMS, because it directly determines the sensitivity, stability, natural frequency, and measurement range of microdevices. Microbeam design always focuses on the dependence of the stiffness on its structure dimensions and material property. The analytical formula of stiffness of folded beam which is used for flexible supporting of mass of microaccelerometers has been given in details by Senturia [14]. The formula was derived by replacing the folded beam by a fixed-fixed beam with the load on its middle position. Wong et al. have given out the theoretical analysis of stiffness constant of round-folded beams in MEMS accelerometers [15]. The mentioned formulae are based on the assumption the beams have a uniform material property or Young’s modulus at least. However, the microstructures always exhibit a varied property after the silicon undergoes a series of doping and deposition processes, specifically in the self-stop etching process where a heavy boron-doping step exists [3].The doping process including boron ion implantation and diffusion will lead to a changed property of doping silicon from the surface to the interior. This paper proposes a finite element method to study the stiffness of boron-doped microfolded beams of in-plane capacitive microaccelerometers. This method utilizes an approximated multilayer model to replace the continuous solid beam. Each layer of the model has a uniform property and Young’s modulus of one layer is different from another, and then mechanics equations and boundary conditions are established analytically to analyze the stiffness of the approximated model under the transverse force vertical to stratification plane. In Section2, the folded beam form in microaccelerometer is introduced, and its material distribution of boron-doped silicon is recognized as an extended Gaussian function along the thickness. Young’s modulus corresponding to the concentration of boron is identified by the characteristics of self-stop etching of silicon. Section 3 gives out the equivalent model by dividing the cross section of beam into a series of rectangles, and its corresponding effective modulus, effective moment of inertia, and bending stiffness are derived by material mechanics equations. The followed simulation and experiments indicate a good agreement between proposed model and simulation. The application condition of this model and conclusion can been found in Section 4.


## 2. Material Property of Beam

The microfabrication, mainly consisting of surface silicon and bulk silicon technologies, always involves diffusion, deposition, ion implantation and etching, and so forth. Those processing steps generally concern the rearrangement or redistribution of materials, especially in the boron-doping silicon of BSDWP [3]. The studied folded beams in this paper are fabricated by BSDWP and applied in microaccelerometers for supporting the movable mass block which can transfer the input acceleration into the shift of differential capacitance. The whole sensor of accelerometer is given in Figure 1(a), and one of its supporting beams is shown in Figure 1(b). The sensitive axis of sensor is parallel to X-axis, and the deformation plane of the folded beams is parallel to X
Y plane under the symmetrical layout of mass block and supporting beams. Therefore, the load along with F is called the in-plane load, while the load perpendicular to X
Y plane is called the out-of-plane load. The folded beam can be replaced by a fixed-fixed one with the load applied on the center according to [14]. Figure 1(c) gives its cross section whose material property is represented by gradually changing colors. The bottom part of the beam has the highest concentration of 4 boron and the top part has the lightest according to [3, 16] in which the authors described the whole process flow of BSDWP in detail.Figure 1
Beam model in the sensor of microaccelerometer.The boron concentration in silicon is assumed to be an extended half Gaussian function of the depth according to the solid diffusion mechanism [17]. The maximum concentration occurs at surface and then gradually decreases. Therefore, only two points of concentration and diffusion depth are needed to predict the boron profile in silicon solid. This condition can be achieved by two separate self-stopped etching processes using different etchants [18]. One etchant is ψ (HF : HNO3 : CH3COOH = 1 : 3 : 8); the stop point is at a boron atom concentration of 1
×
10
17 cm−3. The other etchant is 30%–40% KOH; the stop point is at the concentration of 5
×
10
19 cm−3. Therefore, the Gaussian functions (f
(
z
)) of the boron profile with the boundary conditions are given as f
(
z
s
)
=
5
×
10
19 cm−3 and f
(
z
b
)
=
1
×
10
17 cm−3, where z is the coordinate along depth and z
s and z
b are the distances of stop points from the surface, respectively.The relative concentration, defined as the ratio of boron atoms to silicon atoms per unit volume, can be expressed as(1)
C
(
z
)
=
f
(
z
)
f
Si
(
z
)
=
f
(
z
)
5
×
1
0
22
with the boundary conditions C
(
z
s
)
=
1
×
10
-
3 and C
(
z
b
)
=
2
×
10
-
6.Compared with undoped silicon, heavily boron-doped silicon has a 0.8% higher thermal expansion coefficient [19] and a 20%–30% higher Young’s modulus [20, 21]. Since it is related to the concentration of boron, Young’s modulus can be assumed as
(2)
E
(
z
)
=
E
Si
+
C
(
z
)
E
max
⁡
-
E
Si
C
max
⁡
,
where C
max
⁡ (=8
×
10
-
3) is the maximum relative concentration determined by the solid solubility of boron in silicon, E
max
⁡ (=200–220 GPa) is Young’s modulus to C
max
⁡, and E
Si (=169 GPa) is Young’s modulus of silicon. The values of E
(
z
), therefore, can be calculated under the condition that Young’s modulus is proportional to the boron concentration.


## 3. Equivalent Model and Simulation

The constant of rigidity under the out-of-plane load is considered in this paper in order to evaluate the cross sensitivity of microaccelerometers. The cross section of beam is depicted in Figure2(a); the density of red represents the concentration of boron in silicon, which decreases with the depth from the top surface according to the Fick’s diffusion laws. Figure 2(b) is a simplified model of boron-doping silicon with four separated layers divided by three contour lines. When the beam is divided into a series of layers by the method applied to cross section, a multilayer beam model is gained to replace the continuous structure by such a finite element lamination (FEL) procedure.Figure 2
Multilayer model.According to the material mechanics theory [22], the double-layer beam in Figure 3(a) with different Young’s modulus and same width can by replaced by the one in Figure 3(b) with the same Young’s module but a different width under the transferring condition that the total strains in both beams are equal to each other under the same load. It is concluded that “to replace the original material structure by another equivalent material structure the area must be multiplied by the modular ratio E
2
/
E
1”; therefore, the equivalent width in Figure 3(b) with the same thickness is
(3)
w
2
=
E
2
E
1
w
.Figure 3
Cross section of double-layer beam.The equivalent Young modulus of double-layer beam can be acquired by the same condition of material replacement. The relationship can be expressed as(4)
E
e
w
(
h
1
+
h
2
)
=
E
1
(
w
h
1
+
w
2
h
2
)
,
where E
e is the equivalent Young modulus which equals
(5)
E
e
=
E
1
h
1
+
E
2
h
2
h
1
+
h
2
.When a series of layers are considered, the double-layer FEL procedure should be extended to a multilayers model of which the equivalent Young modulus is(6)
E
e
=
∑
i
=
1
n
E
i
h
i
∑
i
=
1
n
h
i
.Another important parameter involved in the evaluation of out-of-plane stiffness is the moment of inertia of cross section which can be calculated by the following expression:(7)
I
e
=
∑
i
=
1
n
I
i
=
∑
i
=
1
n
(
I
c
i
+
A
i
(
y
c
-
y
c
i
)
2
)
,
where I
e is moment of inertia of cross section, I
i is the moment of inertia of the ith section relative to y
c, I
c
i is the moment of inertia of the ith section relative to its own neutral axis, A
c is the area of the ith section, y
c
i is the coordinate of neutral axis of the ith section, and y
c is the coordinate of neutral axis in the coordinate system defined in Figure 2; its position has been determined in [3]:
(8)
y
c
=
(
1
/
2
)
∑
i
=
1
n
w
i
h
i
2
+
∑
i
=
2
n
(
w
i
h
i
∑
j
=
1
i
-
1
h
j
)
∑
i
=
1
n
w
i
h
i
=
(
1
/
2
)
∑
i
=
1
n
E
i
h
i
2
+
∑
i
=
2
n
(
E
i
h
i
∑
j
=
1
i
-
1
h
j
)
∑
i
=
1
n
E
i
h
i
,
where w
i and h
i are the width and thickness of the ith layer, respectively.Replacing Young’s modulus of the 2nd tonth layer by that of the first layer, the corresponding width of the 2nd to nth layer, w
i, becomes E
i
w
/
E
1. Thus the parameters mentioned above can be expressed as
(9)
I
c
i
=
1
12
w
i
h
i
3
=
E
i
w
h
i
3
12
E
1
,
(10)
A
i
=
w
i
h
i
=
E
i
w
h
i
E
1
,
(11)
y
c
i
=
h
i
2
+
∑
j
=
1
i
-
1
h
i
,
(12)
t
=
∑
i
=
1
n
h
i
.Substituting (9)–(11) into (7), the moment of inertia of cross section of boron-doping silicon beam can be expressed as
(13)
I
e
=
∑
i
=
1
n
(
E
i
w
h
i
3
12
E
1
+
E
i
w
h
i
E
1
(
y
c
-
(
h
i
2
+
∑
j
=
1
i
-
1
h
i
)
)
2
)
.According to [14], the out-of-plane rigidity of folded beam can be written as
(14)
k
=
π
4
E
e
I
e
8
L
3
.The proposed FEL method avoids the time-consumed modeling and complicated meshing process and provides an analytical solution to the out-of-plane rigidity of boron-doping silicon beams based on the material characteristics of diffusion mechanism. This solution procedure can be programmed by any computer language and the model parameters can be adjusted easily through parametric design for many kinds of beams. Figure4 gives the calculation procedure of the FEL method. Because each layer after laminating is assumed to be with a uniform material property, the thickness of layers must be precisely controlled to acquire a more accurate model. This requirement is realized by defining the difference of Young’s modulus between two adjacent layers less than a set error, Δ
E
=
E
(
h
i
)
-
E
(
h
i
-
1
)
<
e.Figure 4
Procedure of rigidity calculation by FEL method.The results of FEL method are compared with those gained in commercial software COMSOL Multiphysics to verify its availability and effectiveness. The boron profile is determined by two diffusion depth points,z
s and z
b, which are gained by the self-stop etching process and their values are measured as 42 μm and 83 μm, respectively [3]. The parameters of simulated beams are shown in Table 1. Young’s modulus of monocrystalline silicon is 169 Gpa, while that of the heavy boron-doping silicon is 200 Gpa. The dimensions of beam are equal to those of the designed accelerometers. The rigidity constants of uniform beams with 169 Gpa and 200 GPa are 1.372 kN/m and 1.624 kN/m, respectively. The constant becomes 1.481 kN/m when the modulus variation induced by boron implantation is considered. These constants acquired by software are shown in Table 2, in which the relative errors are less than 7%. Figure 5 shows the curve of relationship of rigidity constant and layer number n from 100 to 1000, which indicates that a larger number, n, results in a more precise constant comparing with that of COMSOL.Table 1
Parameters of beams.
Property
Beam dimensions
E
0, GPa
E
m, GPa
L, µm
w, µm
h, µm
169
200
420
8
42Table 2
Rigidity constants of beams.
Beams
Young’s modulus
COMSOL
FEL
Relative error
E, GPa
kN/m
kN/m
1
169
1.292
1.372
6.19%
2
200
1.529
1.624
6.21%
3
Gauss
1.417
1.478
4.28%Figure 5
Constant of rigidity varies withn.The FEL method shows an effective solution to the constants of rigidity of microbeams with a nonuniform material distribution. Table3 lists the results from COMSOL and FEL method when the thickness of beams varies with process condition. The relative error decreasing with increasing thickness results from the ignored effects of hinges which, is not considered in simplification of microbeams. When the external force is applied on the folded beam, the single beam connecting to the anchors undergoes both bending and torsional deflection. The later will generate an additional displacement of the loaded point, of which the value is inversely proportional to the torsional stiffness. As the thickness of beam is increasing, the torsional stiffness increases linearly with thickness, while the bending stiffness increases cubically with thickness. Consequently, a thicker beam leads a higher torsional stiffness but a much higher bending stiffness, which means that the contribution of torsional displacement is declined and a more precise rigidity constant of folded beam is acquired with a smaller relative error.Table 3
Rigidity constants of boron-doping beams.
Thickness (µm)
30
35
40
45
50
COMSOL (kN/m)
0.528
0.831
1.228
1.732
2.350
FEL (kN/m)
0.560
0.875
1.285
1.802
2.435
Relative error
6.06%
5.29%
4.64%
4.04%
3.62%


## 4. Conclusion

To overcome the disadvantages of traditional methods, a finite element lamination model is developed to calculate the constant of rigidity of boron-doping silicon beam which shows a nonuniform material property along its thickness. The method divides the continuous beam into a series of layers which are glued rigidly together, and each layer has the same properties which are determined by the characteristics of boron-doping silicon material. The equivalent beams created by proposed method can be recognized as a multilayer structure. Its equivalent Young’s modulus and moment of inertia can be easily acquired by the theory of material mechanics. The constant of rigidity, therefore, can be derived by the small deformation theory. The results show a good agreement with those of commercial software COMSOL Multiphysics. The availability and effectiveness of this method can find its applications in MEMS fields, where the material nonuniformity and dimensional uncertainty are dominating, for a shortened modeling time and a higher calculating precise.The presented method now is just suitable and powerful enough for the beams with a uniform structural dimension and a variable material property in any manner. Besides, the proposed method is not superior to FEM simulation in the time consumed. But it can be programmed in any computer language after reasonable simplification of boron-doped beam which frequently appears in MEMS design. This convenience can save the engineers much time in learning complex software and sharply reduce the requirements of professional knowledge.



---

*Source: 894286-2014-02-27.xml*