# Elastic Properties of Boron-Nitride Nanotubes through an Atomic Simulation Method

Authors: Jixiao Tao; Guangmin Xu; Yuzhou Sun
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

The elastic properties of the boron-nitride nanotubes are studied based on an atomic simulation method that is called atomic-scale finite element method. The Tersoff-Brenner potential is used to describe the interaction between boron and nitrogen atoms, and the computational method is established in an atomic-scale scheme similar to the classical finite element method. Young’s modulus is evaluated for the boron-nitride nanotubes, and their buckling behavior is analyzed. It is shown that the diameter has an obvious influence on Young’s modulus of BNNTs, and the buckling is little related to the length of the nanotubes.


---

## Body

## 1. Introduction

Boron-nitride nanotubes (BNNTs) are a very promising one-dimensional material and have a structural analogy to carbon nanotubes [1–3]. Some research reports have shown that BNNTs possess unique structural, mechanical, thermal, electrical, and chemical properties. For example, Young’s modulus of BNNTs is on the order of 1 TPa [4, 5]. The thermal conductivity along the nanotube is also very high. BNNTs have always large band gaps regardless of the chirality and diameter and are therefore semiconductors which are contrary to carbon nanotubes. The atomic simulation methods are important to the development of nanotechnology and to the study of nanomaterials and nanosystems. Molecular dynamics [6, 7] is a commonly used atomic-scale method, in which the conjugate gradient method is used for the energy minimization that consumes a large amount of computational resources, and is only available for the very small size. The atomic-scale finite element method (AFEM) is proposed by Liu et al. [8, 9], and it can achieve a high computational efficiency with the same accuracy as molecular dynamics. In the present study, AFEM is used to obtain the stiffness matrix and nonequilibrium force vector of the system, and the equilibrium state is determined with the nonlinear iteration.In the field of theoretical research, only a few researches have been reported about the tensile and compressive properties of BNNTs. Oh used the continuum lattice approach to estimate elastic properties of BNNTs, in which the Tersoff-Brenner potential was used to describe the interaction between boron and nitrogen atoms, but some potential parameters have been modified to fit the cohesive energy and the bond length of boron nitride [10]. Song et al. employed an atomic-based continuum theory to study Young’s modulus, stress-strain curve, and nonlinear bifurcation in BNNTs [11]. They pointed out that the mechanical behavior of BNNTs is virtually independent of the diameter and length of BNNTs but has a strong dependence on chirality. Wei et al. used classical molecular dynamics simulations to investigate compressive and tensile behaviors of the carbon nanotubes and boron-nitride nanotubes [12]. From the computational analyses, they found that the chirality is the main factor affecting the behavior of the nanotubes, and the nanotubes in different materials but with the same chirality have similar deformation patterns. Liao et al. investigated the deformation behaviors of an (8, 8) boron-nitride nanotube under axial tensile strains, in which Tersoff potential was employed with the appropriate potential parameters [13]. According to their results, the BNNT starts to fail at the failure strain of 26.7% and the local elongation dominates the tensile failure of the BNNT. In this paper, the widely used Tersoff-Brenner potential is employed in the atomistic simulation, and a set of potential parameters modified by Oh are used to investigate the elastic properties of boron-nitride nanotubes [10]. AFEM is used to obtain the equilibrium states.


## 2. Atomic-Scale Modeling Method

### 2.1. Boron-Nitride Nanotubes

An undeformed BNNT can be visualized as a hollow cylinder that is formed by rolling up a BN sheet into a cylindrical shape. It can be uniquely characterized by a chiral vectorΓ in terms of a set of two integers (
n
,
m
) corresponding to BN sheet unit vectors a
1 and a
2 (Figure 1):(1)
Γ
=
n
a
1
+
m
a
2
.This tube is denoted as an (
n
,
m
) tube with its diameter given by (2)
D
=
3
π
a
B
–
N
n
2
+
m
n
+
m
2
,where a
B
–
N is the bond length in the BN sheet. The tubes with n
=
m are commonly referred to as armchair tubes and those with m
=
0 as zigzag tubes. Other tubes are called chiral tubes in general with the chiral angle θ which is defined as that between the vector Γ and the zigzag direction a
1:(3)
θ
=
t
a
n
-
1
3
m
2
n
+
m
,where θ ranges from 0° for zigzag (m
=
0) to 30° for armchair (m
=
n) tubes (n
≥
m is used for convention). Shown in Figure 2 are Zigzag and armchair BNNTs.Figure 1
The definitions for the basis unit vectors, chiral vector, and chiral angle.Figure 2
The atomic structure of the zigzag and armchair BNNTs.
(a)
(b)


### 2.2. The Atomic-Scale Finite Element

The present research adopts AFEM that was proposed by Liu et al. [8, 9]. The basic idea is to divide nanotubes into finite number of elements, and each element is characterized by a set of discrete atoms. The positions of all atoms are determined by minimizing the energy in the system.For a system ofN atoms, the energy stored in the atomic bond can be denoted by the function of each atom coordinate:(4)
U
t
=
U
t
x
1
,
x
2
,
…
,
x
n
=
∑
i
<
j
N
V
B
x
j
-
x
i
.Tersoff-Brenner potential [14, 15] is a multibody potential and can better describe the interaction between C, B, H, and N atoms. In the present study, Tersoff-Brenner potential V
B
r
i
j is used to describe the interaction between the boron and nitrogen atoms:(5)
V
B
r
i
j
=
V
R
r
i
j
-
B
i
j
V
A
r
i
j
,where V
R and V
A are the repulsive pair potential and attractive pair potential; r
i
j is the distance from atom i to atom j;  B
i
j is the bond order function. The sets of potential parameters modified by Oh [10] are used which are listed in Table 1.Table 1
Potential parameters.
Parameter
D
(
e
) (eV)
S
β
R
(
e
)
R
(
1
) (Å)
R
(
2
) (Å)
δ
a
0
c
0
d
0
6.36
1.0769
22.0
0.133
1.9
2.1
0.382
2.0813
330
3.5The total energy is thus evaluated as(6)
E
t
x
=
U
t
x
-
∑
i
=
1
N
F
i
·
x
i
,where x
=
(
x
1
,
x
2
,
…
,
x
n
)
T, F
i is the external force exerted on atom i. The state of minimal energy corresponds to(7)
∂
E
t
∂
x
=
0
.Giving Taylor expansion of E
t
x and substituting it into (7) yield the following equation:(8)
K
Δ
u
=
P
,where Δ
u is displacement increment and K and P are, respectively, the stiffness matrix and nonequilibrium force vector given by(9)
K
=
∂
2
E
t
∂
x
∂
x
=
∂
2
U
t
∂
x
∂
x
,
P
=
-
∂
E
t
∂
x
=
F
-
∂
U
t
∂
x
.Newton iteration method can be used to solve the present problem. It is much faster than the widely used conjugate gradient method because the first and second order derivatives were used. Materials may display softening behavior when they were under axial compression. For problems involving material softening,K is nonpositive definite and K may be replaced by K
∗
=
K
+
I
α to ensure the convergence, where I is the identity matrix and α is a positive number slightly larger than the absolute value of the minimum negative eigenvalue of the stiffness matrix [16, 17].


### 2.3. The Simulation Process

The above method has been written as a Fortran code for BNNTs, in which the following steps are used to compute the elastic properties of BNNTs and determine their buckling deformation.Step 1.
First, construct the initial configuration of BNNTs with uniform bond length using a separate program; store the coordinates and each piece of bond information in an array.Step 22.
Using the coordinate and bond number of arrays, find the first and second neighbor atoms information and store them in an array.Step 3.
Call a separate program to make the system back to the equilibrium coordinates for given initial coordinates and boundary conditions of BNNTs.Step 4.
Apply displacement field to equilibrium coordinates. This process adopts the constant displacement values at each load step.Step 5.
Store the potential of BNNTs against applied displacement filed. Using the polynomial curve fitting, fit data for equation of potential in terms of displacement filed (strain).Step 6.
Calculate Young’s modulus using equations in Section3.1.Boundary conditions in Step3 are to restrain one side of BNNTs and make the other side free until the system returns to the equilibrium configuration. During the process of axial tensile or compression, one end of the BNNT is completely fixed, and the tensile or compression is achieved by incrementally imposing an axial movement at another end. The length of the tube is changed by 0.01 nm per loading step until material appears buckling.


## 2.1. Boron-Nitride Nanotubes

An undeformed BNNT can be visualized as a hollow cylinder that is formed by rolling up a BN sheet into a cylindrical shape. It can be uniquely characterized by a chiral vectorΓ in terms of a set of two integers (
n
,
m
) corresponding to BN sheet unit vectors a
1 and a
2 (Figure 1):(1)
Γ
=
n
a
1
+
m
a
2
.This tube is denoted as an (
n
,
m
) tube with its diameter given by (2)
D
=
3
π
a
B
–
N
n
2
+
m
n
+
m
2
,where a
B
–
N is the bond length in the BN sheet. The tubes with n
=
m are commonly referred to as armchair tubes and those with m
=
0 as zigzag tubes. Other tubes are called chiral tubes in general with the chiral angle θ which is defined as that between the vector Γ and the zigzag direction a
1:(3)
θ
=
t
a
n
-
1
3
m
2
n
+
m
,where θ ranges from 0° for zigzag (m
=
0) to 30° for armchair (m
=
n) tubes (n
≥
m is used for convention). Shown in Figure 2 are Zigzag and armchair BNNTs.Figure 1
The definitions for the basis unit vectors, chiral vector, and chiral angle.Figure 2
The atomic structure of the zigzag and armchair BNNTs.
(a)
(b)


## 2.2. The Atomic-Scale Finite Element

The present research adopts AFEM that was proposed by Liu et al. [8, 9]. The basic idea is to divide nanotubes into finite number of elements, and each element is characterized by a set of discrete atoms. The positions of all atoms are determined by minimizing the energy in the system.For a system ofN atoms, the energy stored in the atomic bond can be denoted by the function of each atom coordinate:(4)
U
t
=
U
t
x
1
,
x
2
,
…
,
x
n
=
∑
i
<
j
N
V
B
x
j
-
x
i
.Tersoff-Brenner potential [14, 15] is a multibody potential and can better describe the interaction between C, B, H, and N atoms. In the present study, Tersoff-Brenner potential V
B
r
i
j is used to describe the interaction between the boron and nitrogen atoms:(5)
V
B
r
i
j
=
V
R
r
i
j
-
B
i
j
V
A
r
i
j
,where V
R and V
A are the repulsive pair potential and attractive pair potential; r
i
j is the distance from atom i to atom j;  B
i
j is the bond order function. The sets of potential parameters modified by Oh [10] are used which are listed in Table 1.Table 1
Potential parameters.
Parameter
D
(
e
) (eV)
S
β
R
(
e
)
R
(
1
) (Å)
R
(
2
) (Å)
δ
a
0
c
0
d
0
6.36
1.0769
22.0
0.133
1.9
2.1
0.382
2.0813
330
3.5The total energy is thus evaluated as(6)
E
t
x
=
U
t
x
-
∑
i
=
1
N
F
i
·
x
i
,where x
=
(
x
1
,
x
2
,
…
,
x
n
)
T, F
i is the external force exerted on atom i. The state of minimal energy corresponds to(7)
∂
E
t
∂
x
=
0
.Giving Taylor expansion of E
t
x and substituting it into (7) yield the following equation:(8)
K
Δ
u
=
P
,where Δ
u is displacement increment and K and P are, respectively, the stiffness matrix and nonequilibrium force vector given by(9)
K
=
∂
2
E
t
∂
x
∂
x
=
∂
2
U
t
∂
x
∂
x
,
P
=
-
∂
E
t
∂
x
=
F
-
∂
U
t
∂
x
.Newton iteration method can be used to solve the present problem. It is much faster than the widely used conjugate gradient method because the first and second order derivatives were used. Materials may display softening behavior when they were under axial compression. For problems involving material softening,K is nonpositive definite and K may be replaced by K
∗
=
K
+
I
α to ensure the convergence, where I is the identity matrix and α is a positive number slightly larger than the absolute value of the minimum negative eigenvalue of the stiffness matrix [16, 17].


## 2.3. The Simulation Process

The above method has been written as a Fortran code for BNNTs, in which the following steps are used to compute the elastic properties of BNNTs and determine their buckling deformation.Step 1.
First, construct the initial configuration of BNNTs with uniform bond length using a separate program; store the coordinates and each piece of bond information in an array.Step 22.
Using the coordinate and bond number of arrays, find the first and second neighbor atoms information and store them in an array.Step 3.
Call a separate program to make the system back to the equilibrium coordinates for given initial coordinates and boundary conditions of BNNTs.Step 4.
Apply displacement field to equilibrium coordinates. This process adopts the constant displacement values at each load step.Step 5.
Store the potential of BNNTs against applied displacement filed. Using the polynomial curve fitting, fit data for equation of potential in terms of displacement filed (strain).Step 6.
Calculate Young’s modulus using equations in Section3.1.Boundary conditions in Step3 are to restrain one side of BNNTs and make the other side free until the system returns to the equilibrium configuration. During the process of axial tensile or compression, one end of the BNNT is completely fixed, and the tensile or compression is achieved by incrementally imposing an axial movement at another end. The length of the tube is changed by 0.01 nm per loading step until material appears buckling.


## 3. Results and Discussions

### 3.1. Young’s Modulus

In continuum mechanics, the constitutive response between the load and deformation is established prior to solving a specific problem. If the material is homogeneous and isotropic, the material can be represented by two independent constants, namely, Young’s modulusY and Poisson’s ratio ν. For a material undergoing a uniaxial deformation, Y is defined as(10)
Y
=
2
V
δ
-
V
0
V
o
l
L
2
δ
2
=
k
L
2
V
o
l
,
k
=
2
V
δ
-
V
0
δ
2
,where V
o
l is the volume, L is the initial nanotube length, and V
0 and V
δ are the equilibrium (minimised) potentials corresponding to the initial and deformed equilibrium configurations, respectively. δ is the length change in BNNTs, and k is the stiffness of the BNNT, as δ
→
0  k
→
k
0.The thickness is often taken as 0.34 nm [11–13, 18] in the evaluation of V
o
l. V
δ can be expressed as(11)
V
δ
=
a
1
+
a
2
δ
+
a
3
δ
2
+
a
4
δ
3
.The stiffness constant k
0 is then obtained as(12)
k
0
=
lim
δ
→
0
⁡
∂
2
V
δ
∂
δ
2
=
2
a
3
.Deformation behavior of an (8, 0) BNNT under axial tensile strains is first investigated in this paper. It has 23 hexagonal cells along the axis, and its initial length and radius are 5.037 nm and 0.326 nm, respectively. Figure3 shows the stress-strain curve. Plot of equilibrium potential energy versus the length change is displayed in Figure 4.Figure 3
Stress-strain curve of an (8, 0) BNNT.Figure 4
Potential versus the length change of an (8, 0) BNNT for axial tensile/compressive deformation.The potential equation (11) is obtained from Figure 4 using polynomial curve fitting. Replace k
0 with 2
a
3
=
745.8. Substituting (12) into (10), the obtained Y is 863.85 GPa.Using the above steps, armchair(
n
,
n
) and zigzag (
n
,
0
) BNNTs are also simulated. Figure 5 shows the variation of the axial Young’s moduli with the tube radius. It can be seen that Young’s moduli of both armchair and zigzag BNNTs increase with increasing tube diameter, and this trend becomes much obvious for small radii zigzag tubes. With an increasing tube radius, the axial Young’s moduli tend to the same constant. Employing an atomistic-based continuum theory, Song et al. [11] obtained the normalized Young’s modulus of BNNT for the armchair (
n
,
n
) and zigzag (
n
,
0
) BNNTs with the varying tube diameter, and they compared their results with tight binding [4] and ab initio calculation [5]. The present results agree well with those from [4, 5, 10, 11].Figure 5
Young’s moduli versus tube radius for(
n
,
0
) and (
n
,
n
) BNNTs.


### 3.2. Axial Buckling of BNNTs

AFEM is also applied to present a complete numerical simulation of buckling behavior. At each loading step, the stable state is solved with Newton’s method, and then the further compressive displacement is used. The penalty function method [19–21] is used to enforce the essential boundary condition. When BNNTs appear buckling, K is nonpositive definite. The iterative convergence is achieved to replace K with K
∗
=
K
+
I
α, where I is the identity matrix and α is a positive number larger than the minimum negative eigenvalue of the stiffness matrix.The buckling deformation of an (8, 0) BNNT is shown in Figure6 and it is very similar to the single-wall carbon nanotube studied by some researchers. The molecular dynamic method is the most popular atomic simulation method, in which the conjugate gradient method is used to achieve the energy minimization. The present method has the same accuracy with the molecular dynamic method because they both are atomic-scale methods. The conjugate gradient method is an order-N
2 method, and its computational cost is very huge. In the present AFEM, Newton iteration method is applied to obtain the equilibrium state, in which the first and second order derivatives are used and 3–5 iterative steps can achieve a good convergence. The computation is far faster than molecular dynamic method. Liew et al. showed that MD simulation of the buckling behavior of a (10, 10) SWCNT containing 2,000 atoms required 36 hours in a single CPU of SGI origin 2000, whereas the computation for a four-walled MWCNT containing 15,097 atoms took four months [2, 6]. The computational time in AFEM scales linearly with number of atoms and the numbers of iteration steps is approximately independent of the atom number, which implies that AFEM is an order-N method and is very effective for the nanostructure with a larger number of atoms.Figure 6
Buckling deformation.Our results are in good agreement with Wei et al. who have used classical molecular dynamics simulations to investigate compressive behaviors of the boron-nitride nanotubes [12]. In order to investigate the relationship between the bulking and length of single-walled boron-nitride nanotubes, some (
n
,
0
) BNNTs were simulated. Figure 7 shows the buckling strain versus length for several zigzag BNNTs. It is observed that the curve is nearly flat so that the nanotube length has little influence on buckling.Figure 7
The buckling strain versus the length for several zigzag boron-nitride nanotubes.


## 3.1. Young’s Modulus

In continuum mechanics, the constitutive response between the load and deformation is established prior to solving a specific problem. If the material is homogeneous and isotropic, the material can be represented by two independent constants, namely, Young’s modulusY and Poisson’s ratio ν. For a material undergoing a uniaxial deformation, Y is defined as(10)
Y
=
2
V
δ
-
V
0
V
o
l
L
2
δ
2
=
k
L
2
V
o
l
,
k
=
2
V
δ
-
V
0
δ
2
,where V
o
l is the volume, L is the initial nanotube length, and V
0 and V
δ are the equilibrium (minimised) potentials corresponding to the initial and deformed equilibrium configurations, respectively. δ is the length change in BNNTs, and k is the stiffness of the BNNT, as δ
→
0  k
→
k
0.The thickness is often taken as 0.34 nm [11–13, 18] in the evaluation of V
o
l. V
δ can be expressed as(11)
V
δ
=
a
1
+
a
2
δ
+
a
3
δ
2
+
a
4
δ
3
.The stiffness constant k
0 is then obtained as(12)
k
0
=
lim
δ
→
0
⁡
∂
2
V
δ
∂
δ
2
=
2
a
3
.Deformation behavior of an (8, 0) BNNT under axial tensile strains is first investigated in this paper. It has 23 hexagonal cells along the axis, and its initial length and radius are 5.037 nm and 0.326 nm, respectively. Figure3 shows the stress-strain curve. Plot of equilibrium potential energy versus the length change is displayed in Figure 4.Figure 3
Stress-strain curve of an (8, 0) BNNT.Figure 4
Potential versus the length change of an (8, 0) BNNT for axial tensile/compressive deformation.The potential equation (11) is obtained from Figure 4 using polynomial curve fitting. Replace k
0 with 2
a
3
=
745.8. Substituting (12) into (10), the obtained Y is 863.85 GPa.Using the above steps, armchair(
n
,
n
) and zigzag (
n
,
0
) BNNTs are also simulated. Figure 5 shows the variation of the axial Young’s moduli with the tube radius. It can be seen that Young’s moduli of both armchair and zigzag BNNTs increase with increasing tube diameter, and this trend becomes much obvious for small radii zigzag tubes. With an increasing tube radius, the axial Young’s moduli tend to the same constant. Employing an atomistic-based continuum theory, Song et al. [11] obtained the normalized Young’s modulus of BNNT for the armchair (
n
,
n
) and zigzag (
n
,
0
) BNNTs with the varying tube diameter, and they compared their results with tight binding [4] and ab initio calculation [5]. The present results agree well with those from [4, 5, 10, 11].Figure 5
Young’s moduli versus tube radius for(
n
,
0
) and (
n
,
n
) BNNTs.


## 3.2. Axial Buckling of BNNTs

AFEM is also applied to present a complete numerical simulation of buckling behavior. At each loading step, the stable state is solved with Newton’s method, and then the further compressive displacement is used. The penalty function method [19–21] is used to enforce the essential boundary condition. When BNNTs appear buckling, K is nonpositive definite. The iterative convergence is achieved to replace K with K
∗
=
K
+
I
α, where I is the identity matrix and α is a positive number larger than the minimum negative eigenvalue of the stiffness matrix.The buckling deformation of an (8, 0) BNNT is shown in Figure6 and it is very similar to the single-wall carbon nanotube studied by some researchers. The molecular dynamic method is the most popular atomic simulation method, in which the conjugate gradient method is used to achieve the energy minimization. The present method has the same accuracy with the molecular dynamic method because they both are atomic-scale methods. The conjugate gradient method is an order-N
2 method, and its computational cost is very huge. In the present AFEM, Newton iteration method is applied to obtain the equilibrium state, in which the first and second order derivatives are used and 3–5 iterative steps can achieve a good convergence. The computation is far faster than molecular dynamic method. Liew et al. showed that MD simulation of the buckling behavior of a (10, 10) SWCNT containing 2,000 atoms required 36 hours in a single CPU of SGI origin 2000, whereas the computation for a four-walled MWCNT containing 15,097 atoms took four months [2, 6]. The computational time in AFEM scales linearly with number of atoms and the numbers of iteration steps is approximately independent of the atom number, which implies that AFEM is an order-N method and is very effective for the nanostructure with a larger number of atoms.Figure 6
Buckling deformation.Our results are in good agreement with Wei et al. who have used classical molecular dynamics simulations to investigate compressive behaviors of the boron-nitride nanotubes [12]. In order to investigate the relationship between the bulking and length of single-walled boron-nitride nanotubes, some (
n
,
0
) BNNTs were simulated. Figure 7 shows the buckling strain versus length for several zigzag BNNTs. It is observed that the curve is nearly flat so that the nanotube length has little influence on buckling.Figure 7
The buckling strain versus the length for several zigzag boron-nitride nanotubes.


## 4. Conclusions

This paper has used AFEM to study the elastic properties of boron-nitride nanotubes based on interatomic potentials for boron and nitrogen atoms. It is shown that Young’s moduli of both armchair and zigzag BNNTs closely related to tube diameter, especially for small radii zigzag tubes. When radius exceeds 1 nm, Young’s moduli tend to the same constant. For zigzag BNNTs, the buckling strain is virtually independent of the nanotube length, and its average strain is 12%. This is consistent with the conclusions that the mechanical behavior of BNNTs is independent of the diameter and length of BNNTs by some researchers. AFEM is an efficient and accurate computation method and it is also readily applicable for solving many physics related optimization problems.



---

*Source: 240547-2015-10-11.xml*