After rolling, the atomic bonds form a tetrahedron with lengths b1 to b6 (secants of the cylinder). As the ... Ch = Ch. = n a1 + m a2. = √n2a2. 1 + m2a2. 2 + nm (a2. 1 + a2. 2 − a2. 3)§,. (2.1) ..... model a simplified energy expression. Employed ...... nxt += r[h] ∗ d[h];. 253. i f ((k == 8) | | (nxt
Christian Walchshofer
Atomistic and Continuum Simulations for Single Walled Carbon Nanotubes
Atomistic and Continuum Simulations for Single Walled Carbon Nanotubes Christian Walchshofer
Supervisor: O.Univ.-Prof. Dipl.-Ing. Dr.techn. Christian C. Celigoj, Institute for Strength of Materials, Graz University of Technology
Submitted at Graz University of Technology for the degree Diplom Ingenieur (Master of Science)
Graz, June 2007
iii
Abstract
Nanotubes are an essential component of modern nanosystems. Profound knowledge of mechanical properties and deformation behaviour is crucial for further progress and developments. The present work deals with the modelling as well as with the techniques necessary to simulate carbon nanotubes. Special attention is paid to the transformation of the (real existing) multibody system into a continuum system, that can be solved by using the finite element method. The simulation method which is presented is much more efficient to compute in comparison to the conventional (atomistic) multibody system, which allows for simulations of nanosystems in the whole.
Kurzfassung
Nanor¨ ohren sind ein elementarer Bestandteil moderner Nanotechnologien. Tiefgehendes Verst¨ andnis f¨ ur die mechanischen Eigenschaften dieser Strukturen ist von wesentlicher Bedeutung f¨ ur weitere Entwicklungen. Die vorliegende Arbeit befasst sich sowohl mit der Modellierung, als auch mit der notwendigen Simulationsmethodik von Kohlenstoffnanor¨ohren. Beson¨ ders bedeutsam ist bei letzterem die Uberleitung des (reell vorhandenen) Mehrk¨ orpersystems in ein Kontinuumssystem, das mit Hilfe von finiten Elementen simuliert werden kann. Die dargestellte Berechnungsmethodik ist wesentlich effizienter als herk¨ommliche Mehrk¨orpersysteme, was schließlich Simulationen von ganzheitlichen Nanosystemen erm¨oglicht.
5.1. Parameter set for numerical example . . . . . . . . . . . . . . . . . . . .
22
1. Introduction 1.1. Carbon nanotubes Carbon nanotubes (CNTs) are tubes consisting of pure graphene, which exhibit diameters from 0,4 to 100 nm (10−9 m). Those structures have been of great interest for scientists and companies since their discovery in 1991∗, as they show remarkable mechanical, chemical and electrical properties. Young’s modulus for CNTs has been determined by experiments as well as calculations to be around 106 N/mm2 . This feature, combined with a tensile strength of 45 000 N/mm2 and a mass density of 1, 3 kg/dm3 makes CNTs one of the most promising materials for the future in mechanical design. Literature gives coefficients for heat conductance in the order of 6 000 W/mK, which is double as much as that of diamond. This behaviour is strongly interrelated with the fact that CNTs are either conductors or semiconductors, depending on their structure† . The aforementioned material parameters succeed in various applications: CNT fibre reinforced compound materials have been introduced e.g. in cycling competitions‡ . A first nanotube integrated memory circuit was presented in 2004, demonstrating the useability for IT components. Displays, electron microscopy, hydrogen repositories (tanks), Nanomachines (-motors), cooling elements and solar cells are just a few CNT applications that have already been implemented or are expected to be implemented in the next future.
∗
S. Iijima, Helical microtubules of graphite carbon, Nature 354, p. 56 - 58. Chirality, see Sec. 1.2 ‡ Floyd Landis’ winning bike at the Tour de France 2006 †
2
Figure 1.1.: SEM (Scanning electron microscopy) image of CNTs Image courtesy of Micro- and Nanotech. Laboratory, University of Illinois at Urbana-Champain http://www.micro.uiuc.edu/Equipment/hitachi4800.htm
1.2. Chirality Nanotubes subdivide into a variety of classifications. This thesis deals exemplarily with single walled carbon nanotubes (SWNTs), but there are also multi walled - as well as non-carbon tubes§ . SWNTs are further subdivided in terms of chirality (n,m)¶ . A so called “chiral” vector Ch = n a1 + m a2 (where n ≥ |m| ≥ 0) is determined by the plane sheet of graphene,
see Fig. 1.2. After rolling this sheet into a SWNT, the chiral vector encompasses the circumference of the tube in such a manner that the vector’s head matches its tail. The SWNT is defined if the integers n and m are known, which classify tubes as Zigzag (n,0), Armchair (n,n) and Chiral (n,m) types, see Fig. 1.3. The structure of SWNTs affects mechanical and electrical properties: Chapter 3 shows differing coupling of extension and twist. Jiang et al. [7] give distinct stiffnesses for the three SWNT types. It is also interesting that SWNTs behave like conductors if the fraction
n−m 3
is an integer, while they are semiconducting if it is a floating point
number.
§ ¶
E.g. silicon, boron nitride and titanium dioxide nanotubes have already been synthesized. Greek: Handedness; a system is called “chiral” if its mirrored image does not match with, nor superimpose on the original object.
3
Figure 1.2.: Chirality of nanotubes Image courtesy of Wikipedia Commons, http://commons.wikimedia.org/wiki/Image:CNTnames.png, 2005
Figure 1.3.: Zigzag, Armchair and Chiral SWNTs Image courtesy of Saito Laboratory, Department of Quantum Engineering, Nagoya University www.surf.nuqe.nagoya-u.ac.jp
4
1.3. Outline This thesis deals with the deformation behaviour of carbon nanotubes. Chapter 2 refers to the calculation of the equilibrium geometry for SWNTs. Mapping of a graphene sheet into a SWNT as well as the interatomic potential for carbon are discussed. The coupling of extension and twist is shown in Chapter 3. We apply a direct map to calculate an energy potential, which is minimized over the relevant degrees of freedom. Numerical results and their interpretations are giving insights on the mechanics of nanotubes. Chapter 4 sketches the fundamentals required for a continuum formulation for nanotubes. The standard treatment for space filling crystals is shown, its deficiencies are discussed and an adaption for bodies of reduced dimensionality (e.g. nanotubes) is given. The implementation of the constitutive model that has been derived is presented in Chapter 5. We develop an atomistic (multibody) description for chains in 2D at first, to be able to create reliable benchmark tests. Then, we develop a finite element model that incorporates the given theory. For this task, we will use Hermitean finite elements. Finally, results of both simulations are displayed, compared and discussed. Here, special attention is paid to the observed locking behaviour of the finite element model.
2. Equilibrium for unloaded, single walled carbon nanotubes The equilibrium configuration of arbitrary single walled carbon nanotubes is of major importance as a starting point for further calculations. This chapter describes how to calculate atomic equilibrium positions on SWNTs in an efficient way∗ . The underlying procedure is as follows: We consider a plane, relaxed sheet of graphene, which displays a hexagonal lattice structure. This lattice can be reduced into recurring triangular sublattices, that are defined by 5 independent lengths ai (with i = 1,2,...,5), see Fig. 2.1. The graphene sheet is then rolled up to obtain a tube, on which the end of the circumferential (chiral) vector Ch = n a1 + m a2 † is identical with its origin. After rolling, the atomic bonds form a tetrahedron with lengths b1 to b6 (secants of the cylinder). As the length of bonds and the angles in between change throughout this process, the bond potential changes as well, and the equilibrium positions have to be recalculated. For reasons of periodicity, only one representative triangle BCD has to be taken into account for further investigations.
2.1. Mapping of plane graphene into a SWNT The calculation of equilibrium positions is done by the help of a mapping. As the SWNT is defined by the lengths ai as well as its chirality (n,m), all geometry data can be expressed by using these input values. The lengths ai can be considered as the system’s degrees of freedom, and can be used for energy minimization. ∗ †
Deductions in this Chapter are closely following the work of Jiang et al. [7] n and m are integers that define the chirality of a SWNT
6
D a2
a3 ϕ1
D b3 B
b6
A
b1
ϕ2
θ
C a1
Ch
B b5
b4
A a5
a4
C b2
Figure 2.1.: Graphene sheet map and tetrahedron on a SWNT
We start mapping‡ by defining the length of the chiral vector Ch = k Ch k = k n a1 + m a2 k =
q
n2 a21 + m2 a22 + n m (a21 + a22 − a23 )§ ,
(2.1)
which defines the SWNT diameter dt =
Ch . π
(2.2)
The relevant angles on the map are ϕ1 = arccos
a21 + a23 − a22 , 2 a1 a3
(2.3)
ϕ2 = arccos
a21 + a24 − a25 , 2 a1 a4
(2.4)
θ = arccos
‡ §
na21 +
m 2
For the derived expressions see Fig. 2.1 Note that 2 a1 · a2 = a21 + a22 − a23 (law of cosines)
(a21 + a22 − a23 ) . a1 Ch
(2.5)
7
Now cylindrical coordinates for the atoms A,B,C and D in tubular configuration can be obtained ZA = a4 sin (ϕ2 + θ) ,
ZB = 0 ,
ZC = a1 sin θ ,
ZD = a3 sin (ϕ1 + θ) , (2.6)
ΘA =
2a4 cos (ϕ2 + θ) , dt
ΘB = 0 ,
ΘC =
2a1 cos θ , dt
ΘD =
2a3 cos (ϕ1 + θ) , dt (2.7)
which are used for the evaluation of secants (bond lengths) bi ¶ bXY =
r
d2t [1 − cos (ΘY − ΘX )] + (ZY − ZX )2 . 2
(2.8)
Angles between bonds are also needed for the calculation of bond energies. They are γABC = arccos
b24 + b25 − b21 , 2 b4 b5
(2.9)
γACD = arccos
b25 + b26 − b22 , 2 b5 b6
(2.10)
γADB = arccos
b24 + b26 − b23 . 2 b4 b6
(2.11)
2.2. Brenner’s interatomic potential for carbon With the geometry data from section 2.1 we need to calculate bond energies. This is done using Brenner’s interatomic potential for carbon [3]. Brenner defines the interatomic potential for atoms i and j as a superposition of repulsive (VR ) and attractive (VA ) terms V (rij ) = VR (rij ) − Bij VA (rij ) , ¶
Notation is given for arbitrary atoms X and Y in this case.
(2.12)
8
which are defined as VR (r) =
D (e) −√2Sβ (r−R(e) ) e fc (r) , S −1
(2.13)
VA (r) =
D (e) S −√ S2 β (r−R(e) ) e fc (r) . S−1
(2.14)
D (e) , S, β and R(e) are parameters, for which two different sets are given. Set 1 fits for the bond lengths measured in experiments, while the second one matches better with observed stretching force constants. fc (r) is a smooth cut-off function, which truncates the interaction of distant atoms. It is given by
fc (r) =
1 1 2
1 + cos
0
π (r−R(1) ) R(2) −R(1)
r < R(1) R(1) < r < R(2)
(2.15)
r > R(2)
We set fc (r) = 1 because we do not model distant atoms but exclusively the 3 bonds AB,AC and AD in figure 2.1. To risk modelling these bonds with zero assigned potential would not make sense, but cause numerical problemsk as well as raise the numerical cost. In this context it has to be mentioned that carried out calculations show, that fc (r) has no significant influence on the outcome. The term Bij is a factor that weakens the attractive term in dependence of bond angles with neighbouring atoms k. It is given by
Bij =
Bij + Bji , 2
Bij = 1 +
X
k(6=i,j)
G (γijk ) fc (rik ) ,
(2.16)
where G (γijk ) is defined as G(γijk ) = a0
"
c20 c2 1 + 02 − 2 d0 d0 + (1 + cos γijk )2
#
.
Constants a0 , c0 and d0 are given as well in Brenner’s parameter sets. k
flat spots and volatile derivatives
(2.17)
9
2.3. Trigonometric relations D b3 b6 B1
b4
b2 A b5 b1
C
B B2 Figure 2.2.: Trigonometric relations on a SWNT
For the calculation of Brenner’s multibody coupling term Bij , we use the antisymmetry concerning the bonds, here shown for b4 . Fig. 2.2 shows the equivalency of angles γABC with γBAB1 and γABD with γBAB2 , which results in equal factors G(γijk ). As a consequence, the averaging of factor Bij in Eq. (2.16) can be avoided, only 2 neighbouring bonds have to be considered. There is also a physical interpretation for this circumstance. The factor Bij can be considered as an additional torsional stiffness between bonds, because shear between bonds alters the state of energy. Evaluating Bij solely at point A indicates, that all torsional stiffness is assigned to angles γABC , γACD and γABD . Hence we can consider the representative triangle as mechanically uncoupled from the SWNT during deformation (relaxation).
10
2.4. Minimization of energy To calculate equilibrium positions for the atoms, we minimize the energy W = V (rAB ) + V (rAC ) + V (rAD ) = W (a1 , ..., a5 )
(2.18)
over the degrees of freedom (a1 , ..., a5 ) by using conjugate gradients coupled with a simple, but precise line search. As starting values for the calculation of the degrees of freedom we assign the respective values taken from a relaxed, plane sheet of graphene (compare: Tab. 2.2, values for graphene) Dof a1 a2 a3 a4 a5
Table 2.1.: Starting values for energy minimization An abbreviated version of the C code that was employed in this thesis can be found in Appendix B.1.
2.5. Results The table below gives equilibrium positions for Brenner’s parameter set 2. Chirality Diameter (n,m)
dt [nm]
graphene
Bonds b4 [nm]
Bonding angles
b5 [nm]
b6 [nm]
a2 [nm]
a4 [nm]
a5 [nm]
∞
0.145070 0.145070
0.145070
120.000
120.000
120.000
-7.37562
0.251269
0.251269 0.251269 0.145070
0.145070
(30,0)
2.402061
0.145109 0.145109
0.145081
119.802
119.963
119.963
-7.36805
0.251543
0.251293 0.251293 0.145158
0.145158
(20,0)
1.603585
0.145160 0.145160
0.145096
(10,0)
0.807711
0.145443 0.145443
0.145176
119.553
119.917
119.917
-7.35855
0.251891
0.251330 0.251330 0.145272
0.145272
118.198
119.686
119.686
-7.30673
0.253750
0.251547 0.251547 0.145896
0.145896
(7,0)
0.571028
0.145857 0.145851
0.145277
116.281
119.401
119.406
-7.23351
0.256277
0.251910 0.251909 0.146792
0.146786
(5,0)
0.415220
0.146691 0.146691
0.145422
112.586
118.991
118.991
-7.09231
0.260890
0.252777 0.252777 0.148562
0.148562
(18,18)
2.495036
0.145088 0.145114
0.145088
119.866
120.016
119.866
-7.36862
0.251392
0.251392 0.251320 0.145093
0.145196
(12,12)
1.664610
0.145114 0.145172
0.145114
119.699
120.036
119.699
-7.35985
0.251552
0.251551 0.251391 0.145125
0.145357
(6,6)
0.835639
0.145251 0.145490
0.145250
118.800
120.165
118.800
-7.31237
0.252409
0.252408 0.251791 0.145297
0.146235
(5,5)
0.697974
0.249804 0.249804
0.252044
118.274
120.255
118.274
-7.28444
0.252909
0.252909 0.252044 0.145396
0.146759
(4,4)
0.560704
0.145478 0.146037
0.145478
117.311
120.448
117.311
-7.23292
0.253824
0.253825 0.252541 0.145578
0.147740
(25,9)
2.442463
0.145098 0.145114
0.145082
119.822
120.002
119.913
-7.36830
0.251493
0.251317 0.251309 0.145120
0.145188
(16,6)
1.578610
0.145140 0.145178
0.145103
119.574
120.008
119.789
-7.35804
0.251806
0.251395 0.251373 0.145191
0.145358
(9,6)
1.049730
0.145205 0.145330
0.145164
119.133
120.082
119.365
-7.33569
0.252224
0.251769 0.251572 0.145264
0.145786
(9,3)
0.871116
0.145318 0.145436
0.145179
118.558
120.014
119.355
-7.31701
0.253106
0.251669 0.251618 0.145504
0.146020
(6,2)
0.585661
0.145644 0.145918
0.145309
116.716
120.071
118.608
-7.24274
0.255360
0.252238 0.252116 0.146068
0.147247
(4,2)
0.433551
0.146034 0.146745
0.145554
114.186
120.589
116.918
-7.12834
0.257809
0.253846 0.253294 0.146571
0.149500
γABD
[◦ ]
Degrees of freedom a1 [nm]
γACD
[◦ ]
Potential [eV]
γABC
[◦ ]
a3 [nm]
Table 2.2.: Equilibrium positions on a nanotube, parameter set 2
11
3. Coupling of extension and twist As we look at the structure of nanotubes it is apparent, that not only their material behaviour is non-linear, but that they are also anisotropic. Considering this fact it is likely, that they will twist under imposed extension and vice versa, i.e. they will extend under imposed twist. The methodology for the calculation of deformation is very similar to Chapter 2. First we need cylindrical coordinates for the atoms B, C, D ∗ and A for the case of an unloaded SWNT† , as calculated in the previous chapter. We employ a direct map according to [5] (see Fig. 3.1), using the following 5 parameters: Axial strain ǫ, twist κ, radial shift u (assumed to be equal for all atoms), angular shift ηΘ for atom A, and axial shift ηZ , also concerning atom A. A shift vector η for atom A is necessary, because the lattice structure of nanotubes is not centrosymmetric, which causes unbalanced forces [7]. Now we expose the SWNT to a certain axial strain ǫ. The undeformed configuration is mapped into the deformed one, and the bond lengths (and angles) can be recalculated. Once again, we formulate the energy for the representative triangle BCD (Fig. 2.1), and minimize it over the remaining degrees of freedom (κ, u, ηΘ and ηZ ). The process is analogous for the imposed twist case, where κ is prescribed, and ǫ is the sought degree of freedom.
∗ †
The representative triangle In their state of equilibrium
13
3.1. The direct map for coupling of extension and twist The formulation of the direct map r = Φ (R), which maps an arbitrary atom I into its deformed configuration i, is straightforward and is visualized by Fig. 3.1. R 7→ r = R + u ,
Θ 7→ θ = Θ + κZ ,
Z 7→ z = (1 + ǫ) Z
(3.1)
For atom A, we have to incorporate the aforementioned shift vector [ηΘ , ηZ ] R 7→ r = R + u ,
Θ 7→ θ = Θ + κZ + ηΘ ,
Z 7→ z = (1 + ǫ) Z + ηZ
(3.2)
With the obtained (mapped) coordinates, we can formulate the potential energy for the representative triangle as shown in Section 2.1 (Eqs. (2.8) ÷ (2.17)).
eY , ey R(I)
κZ (I)
Θ(I)
B
r(i)
(undeformed) i (deformed)
I u
eX , ex Figure 3.1.: Direct map for extension / twist coupling
r(i)
= Φ [R(I)] = [R(I) + u] [cos (Θ(I) + κZ(I))] ex + [R(I) + u] [sin (Θ(I) + κZ(I))] ey + (1 + ǫ) Z(I) ez
r(i, j) = r(j) − r(i)
r 2 (i, j) =|r(i, j)|2 =
(3.3) (3.4)
2 [R(I) + u]2 [1 − cos (Θ(J) − Θ(I)) cos (κZ(J) − κZ(I))]
+ 2 [R(I) + u]2 [
sin (Θ(J) − Θ(I)) sin (κZ(J) − κZ(I))]
+ (1 + ǫ)2 (Z(J) − Z(I))2
(3.5)
14
3.2. Numerical treatment The implementation written in C language can be found in Appendix B.2. A gradient method‡ coupled with a modified (bracketing method-)§ line search algorithm is employed for energy minimization. Convergence for this solver is slow, because the gradient algorithm is diffusive on the one hand, and the bracketing method takes a lot of function evaluations on the other hand. Anyhow, there are crucial reasons for the application of this solver, namely stability and precision. The simulation of large deformations often comes close to instabilities, where structured solvers loose their orientation. Sophisticated solvers¶ are also in danger of loosing their path if second derivatives fluctuate intensively. Analyses show, that both inconveniences are relevant to the present energy formulation, which averts the use of quicker solvers.
3.3. Results and discussion Fig. 3.2(a) shows calculated deformations for the representative triangle of an Armchair SWNT under imposed twist. The necessity for an additional shift vector [ηΘ , ηZ ] for atom A is obvious: Whilst atoms B, C and D are performing homogeneous deformationsk within the circumferential plane, atom A is not, but relaxes to a point of lower potential. The direct map without shift vector would not be able to display such behaviour. Figures 3.3 to 3.6 compare the macroscopic behaviour of parameter sets 1 and 2, Zigzag, Armchair and Chiral SWNTs∗∗ under imposed extension and twist. The initial configuration of Armchair (7,7) and Zigzag (13,0) tubes is symmetric, while this is not the case for the Chiral (9,6) tube. This results in a non-symmetric, macroscopic deformation behaviour for Chiral SWNTs, as can be seen in all given graphs. ‡
A “canned algorithm” can be found in [10]. Bracketing method for minima, proposed by [8]. ¶ E.g. Newton Raphson line searches k ǫ, κ and u are homogeneous. ∗∗ To compare influences of chirality, SWNTs must exhibit comparable shear in circumferential plane (for same values of κ), thus SWNTs of approximately the same diameter are displayed. §
15
Differences in Brenner’s parameter sets 1 and 2 can be seen regarding the (13,0) Zigzag SWNT. For set 1, the Zigzag tube elongates when exposed to twist (Fig. 3.3(a)), while set 2 shows the contrary behaviour (Fig. 3.5(a)). Different parametrizations for factor Bij within the interatomic potential formulation are responsible for this behaviour. Parametric studies and comparison with experimental data would be necessary to prove applicability of Brenner’s potential formulation to nanotubes. There is another interesting fact that can be observed when twisting Zigzag SWNTs, regardless of the employed parameter set: Fig. 3.2(b) shows, that if the angle γACD is reduced, the system responds with a relatively large elongation of bond AD. From a mechanical point of view, this system can be considered very stiff concerning the bending of bonds one to another, but relatively soft regarding extensional stiffness.
0.25
Κ=-0.70 Κ=-0.35 Κ= 0.00 Κ= 0.35 Κ= 0.70
Z @nmD
0.2 0.15 0.1
A
C
D
0.2 Κ=0 Κ=0.235 Κ=0.47 Κ=0.7
0.15 Z @nmD
D
0.1 A 0.05
0.05
C
B
B
0 -0.05
0
0
0.05 0.1 0.15 0.2 0.25 Θ @nmD
(a)
0
0.05
0.1 0.15 Θ @nmD
0.2
0.25
(b)
Figure 3.2.: Deformation visualization for imposed twist case, Armchair (7,7) (a) and Zigzag (13,0) SWNT (b), Brenner parameter set 1
16
0.005 H9,6L H7,7L H13,0L
0.015
0.01
0.005
Radial shift u @nmD
Axial strain Ε
0.02
0 -0.6
-0.4
-0.2
0
0.2
0.4
0 -0.005 -0.01 -0.015 -0.02
H9,6L H7,7L H13,0L
-0.025 -0.03 -0.6
0.6
-0.4
Twist Κ @radnmD
-0.2
0
0.2
0.4
0.6
Twist Κ @radnmD
(a)
(b)
Figure 3.3.: Axial strain ǫ and radial shift u for imposed twist κ, parameter set 1 0.02
H9,6L H7,7L H13,0L
Twist Κ @radnmD
-0.002 -0.004 -0.006 -0.008 -0.01 -0.012
Radial shift u @nmD
0 H9,6L H7,7L H13,0L
0.01 0 -0.01 -0.02 -0.03
-0.014 0
0.1
Axial strain Ε
(a)
0.2
0.3
0
0.1
0.2
0.3
Axial strain Ε
(b)
Figure 3.4.: Twist κ and radial shift u for imposed axial strain ǫ, parameter set 1
17
0.02
Axial strain Ε
0.01 0.005 0
-0.6
-0.4
-0.2
0
0.2
0.4
0
Radial shift u @nmD
H9,6L H7,7L H13,0L
0.015
-0.01
-0.02 H9,6L H7,7L H13,0L
-0.03
-0.6
0.6
-0.4
Twist Κ @radnmD
-0.2
0
0.2
0.4
0.6
Twist Κ @radnmD
(a)
(b)
Figure 3.5.: Axial strain ǫ and radial shift u for imposed twist κ, parameter set 2
0.02
H9,6L H7,7L H13,0L
-0.005 -0.01 -0.015 -0.02 -0.025 0
0.1
Axial strain Ε
(a)
0.2
0.3
Radial shift u @nmD
Twist Κ @radnmD
0
H9,6L H7,7L H13,0L
0.01 0 -0.01 -0.02 -0.03 -0.04 0
0.1
0.2
0.3
Axial strain Ε
(b)
Figure 3.6.: Twist κ and radial shift u for imposed axial strain ǫ, parameter set 2
4. Continuum formulation for nanotubes The following chapter formulates a continuum description for SWNTs, following closely the derivations given by Arroyo [1].
4.1. Standard Cauchy-Born rule for space filling crystals The methodology for the establishment of continuum models for space filling crystals is intuitive and obvious. We consider a single elementary cell, which can be thought of as a finite, but small volume. At first, the lattice points (respectively vectors) are mapped into the deformed configuration. An expression for the continuum strain energy density can now be obtained by calculating the bond potentials for the single cell, and dividing them by the cell’s volume. Review this operation once again, now in reverse order, from the continuum- to the atomistic description: Take an arbitrary point within the undeformed continuum body, and pretend it is an atom. We model a unit cell around this atom, map the lattice vectors into the deformed configuration, and calculate their bond potentials. Their sum of energy is now spread over the affiliated volume to obtain the strain energy density, which is minimized to identify the systems state of equilibrium. The reversed point of view that is given makes clear, that we do not necessarily model a congruent system to the original one, as the point we have chosen is not generally coinciding with a real atom.
19
One thing that we see, is that the material behaviour of examined crystallites is depending on the geometry of the unit cell, as well as on the formulation of bond potentials. In addition to the material parameters, we need a mapping from the initial configuration into the deformed one. In most cases, this mapping is realized in the form of the Cauchy-Born rule a = F A,
(4.1)
where a is the mapped lattice vector, A is its initial counterpart and F is the deformation gradient (Fig. 4.1). Φ F B
γ
b
Γ
a
A
Figure 4.1.: Cauchy-Born rule The standard Cauchy-Born rule comes with certain difficulties. The deformation gradient maps dX to dx : dx = F dX, where F=
∂ [Φ (X)] ∂x = ∂X ∂X
(4.2)
In the present context, we use it for the mapping of (finite) lattice vectors. It has to be mentioned that F consists of the first derivatives of the deformation map Φ with respect to X, which points out that we are considering a linear mapping. It can be shown that F is exact for the mapping of finite vectors, if the deformation is homogeneous. This condition is fitting to “the central hypothesis behind molecular theories at finite strains”, namely that “at the scale of atomic spacing, the deformation of the crystal is homogeneous”∗ . Finally, for the calculation of bond potentials we need to extract bond lengths (a,b) and angles (γ) after they have been mapped. This can be done using the right Cauchy-Green deformation tensor C = FT F and the relations k a k= ∗
Arroyo [1]
√
AT C A ,
cos γ =
AT C B . kakkbk
(4.3)
20
4.2. The exponential Cauchy-Born rule In the case of surfaces in 3D and curves in 2D the standard Cauchy-Born rule becomes inoperable. For surfaces of reduced dimensionality it is not clear how to define a homogeneous deformation, which makes it impossible to map finite vectors exactly into a deformed configuration when using exclusively the deformation gradient F. The deformation gradient itself causes further troubles, because it is only capable of mapping tangent spaces of the undeformed into tangent spaces of the deformed surface, see Fig. 4.2. This problem can be visualized by a small thought example: Imagine a plane sheet of graphene, as shown in chapter 2. This sheet is rolled into a tubular configuration without stretch. As we have already seen, the bond lengths and angles change throughout this process, as the bonds are chords between atoms on the surface. Using the deformation gradient to map the plane sheet into the deformed configuration, we face difficulties: Two bonds emanating from the same atom lie in the same tangent plane, before as well as after the mapping. This fact shows that bond angles are not changed. Further, the deformation gradient is blind to the fact, that mapped vectors are chords on a surface. As the sheet is rolled without stretch, bond length A remains unchanged too. Noticing that no change in bond potentials is caused by this mapping† , we see that the system shows zero bending stiffness in this case. Hence, use of the standard Cauchy-Born rule fails for arbitrary deformation of surfaces in 3D (respectively curves in 2D).
F?
A
=⇒ a
Figure 4.2.: Difficulties with standard Cauchy-Born rule applied to surfaces
†
Neither bond lengths nor angles in between are changed by this mapping.
21
Arroyo [1] presents the exponential Cauchy-Born rule F to overcome these limitations. Fig. 4.3 sketches the proposed procedure: First, we are mapping the undeformed bond A into the tangent plane to obtain W. This is done by the inversion exp−1 X (A) of the so called exponential map. W can now be mapped legitimately into w by using the deformation gradient F. Finally, we obtain a by applying the exponential map expΦ(X) to vector w. In a short form, we can write FX (A) = expΦ(X) F (X) exp−1 X (A) .
a = FX (A) ,
(4.4)
expΦ(X) W w
X A
Y
F
=⇒
y x a
exp−1 X Figure 4.3.: Exponential Cauchy-Born rule Formulations of the exponential maps expΦ(X) and exp−1 X are not trivial for surfaces in 3D‡ . For curves in 2D, the expressions can be determined in the following way: The curve is locally replaced by its osculating circle. We can now obtain an expression for the ⌢
exponential map expΦ(X) by rolling vector w onto the circle (arc xy) to determine the bond vector a. An expression for the inverse map exp−1 X can be derived accordingly. A complete formulation for the map expΦ(X) can be found in the following chapter.
‡
For further discussion see [1].
5. Numerical example - Chain in 2D To test the validity of Arroyo’s model, we simulate a 1 dimensional chain deforming in 2 dimensions. Sec. 5.2 describes a “quasi multi-body simulation” of the atomistic system, which supplies accurate results needed to benchmark the continuum simulation derived in Sec. 5.3. For reasons of straightforwardness, we do not use Brenner’s interatomic potential, but model a simplified energy expression. Employed parameters therefore are represented by terms known from macroscopic continuum mechanics (E,I,Area ), which is convenient, though mathematically equivalent∗ . The system’s energy potential is defined as Πint
# n " # " nB A X X E Area (ai − Ai )2 E I (γi − Γi )2 = + , 2 A 2 A i i i=1 | {z } i=1 | {z } stretching potential
(5.1)
bending potential
with ai denoting the length of bond i, and γi denoting the bond angle at atom i. The parameter set for Eq. (5.1) that was employed in this work is given in the table below. Parameter nB nA Area E I Ai Γi
Description Number of bonds Number of atoms Cross section area Young’s modulus Moment of inertia Equilibrium bond length Equilibrium bond angle
Value 122 2.06 e5 123 98.135349 0
Unit [1] [1] [mm2 ] [N/mm2 ] [mm4 ] [mm] [rad]
Table 5.1.: Parameter set for numerical example ∗
Note: The present formulation is not equivalent from a physical point of view; e.g. atomic bonds do not display any thickness.
23
5.1. Incremental - iterative method for large deformations As we are considering large deformations, we need to apply an iterative method for the calculation of balanced conditions. The relation for the minimization of energy† t+∆t
δΠ = t+∆t δΠint + t+∆t δΠext = 0
(5.2)
ends up in a set of nonlinear equations, which has to be linearized and solved iteratively‡ . We use a first-order Taylor series approximation for t+∆t
t+∆t
δΠint
δΠint = t δΠint + ∆t δΠint ,
(5.3)
to obtain the linear relation ∆t δΠint = −t δΠint − t+∆t δΠext .
(5.4)
This can be written in the form of: δxT K(i) ∆x(i) = δxT
(i) fext − fint ,
(5.5)
where x(i) is the vector of nodal positions (degrees of freedom) at iteration i, δx denotes (i)
its variation and ∆x(i) its increment. K(i) is the global stiffness matrix and fext − fint is the residual force vector. This system is solved for ∆x(i) , and new nodal positions can be found by t+∆t (i)
x
= t+∆t x(i−1) + ∆x(i) .
(5.6)
As we have solved a linearized system, Eq. (5.2) is not yet equally satisfied, and t+∆t x(i) is not sufficiently close to the final position
t+∆t (n)§
x
. We have to recompute K(i) and
(i)
fint for the new configuration, and solve (5.5) and (5.6) again. This is done repeatedly (n times), until an abortion criterion is applied. In the present work we use a criterion †
The upper left index t + ∆t represents the system’s state after equilibrium has been achieved; t indicates the system’s actual configuration. ‡ Superscript (i) will indicate any iteration § Superscript (n) indicates the final iteration.
24 (i) T
with respect to the “unbalanced energy” ∆x ∆x(i)
(i) fext − fint ≤ ǫ
T
(i) fext − fint ǫ = 10 e−20
(5.7)
f t + ∆t t+∆t
fext
fint
unbalanced energy for 3rd iteration t
fext
t ∆x(1) t
x
t+∆t
t+∆t (1)
x
x
x
Figure 5.1.: Incremental iterative method Fig. 5.1 illustrates the described method for a 1 dimensional case. It is obvious, that the procedure is mathematically equivalent to the Newton-Raphson method. Just like Newton-Raphson, the incremental iterative method converges quadratically in a neighbourhood of the solution, which can be seen in Fig. 5.2.
1010 (i)
fext − fint
105 100 10- 5
∆x(i)
T
quadratic fit
10-10 10-15 10-20 1
2
3
4
5
6
7
Iterations
Figure 5.2.: Convergence behaviour of the incremental iterative method; exemplary case taken from a computation of the multibody system described in the next chapter
25
5.2. Atomistic simulation i+1 ai i
αi+1 ai+1
i−1
αi
Figure 5.3.: Multibody simulation element In this section we will derive the formulations necessary for the simulation of the multibody system. We need to determine the components of Eq. (5.4) ∆t δΠint = −t δΠint − t+∆t δΠext to obtain the desired set of linear equations. For this purpose, we employ a version of the potential Πint (see Eq. (5.1)), which is modified to fit the multibody element depicted in Fig. 5.3 Πint =
X i
(
" # ) E Area 1 (ai − Ai )2 (ai+1 − Ai )2 E I (γi − Γi )2 + , + 2 2 Ai Ai 2 Ai
(5.8)
where γi = αi+1 − αi . The terms of Eq. (5.4) are therefore evaluated to: δΠint ∆δΠint
δΠext
X E Area 1 EI = [(ai − Ai ) δai + (ai+1 − Ai ) δai+1 ] + (γi − Γi ) δγi A 2 A i i i X E Area 1 = [∆ai δai + (ai − Ai ) ∆δai + ∆ai+1 δai+1 A 2 i i EI (∆γi δγi + (γi − Γi ) ∆δγi ) + (ai+1 − Ai ) ∆δai+1 ] + Ai X = − δxT fext . i
, (5.9) (5.10) , (5.11)
26
To expand Eqs. (5.9) ÷ (5.11) further, we need to find expressions for the bond lengths ai , the bond angles γi, and their respective variations in terms of the nodal point coordinates. We define
xi−1
yi−1 xi , yi xi+1 yi+1 | {z } xi
which allows us to express the bond length a as the Euclidean norm of vector a ai = k ai k =
q
aTi
ai ,
ai+1 = k ai+1 k =
q
aTi+1 ai+1 .
(5.13)
The following derivations will be given exclusively for bond ai . The expressions for bond ai+1 will end up analogously, with BI replaced by BII . The first variation δai is given as
δaT ai = δaTi ni , δai = p i T ai ai
with
δai = BI δxi
and
δxi =
δxi−1
δyi−1 δxi . δyi δxi+1 δyi+1
(5.14)
Note that ni is the unit vector pointing into the direction of bond ai . The increment of the first variation ∆δai follows as I − ni nTi ∆δai = δaTi ∆ni = δaTi p ∆ai , aTi ai
where I is the second order identity matrix, and δ is commutative with ∆¶ . ¶
Schwarz’s theorem
(5.15)
27
Expressions concerning the bond angle γi are somewhat more laborious to obtain. γi is defined as γi = αi+1 − αi ,
(5.16)
where αj (j = i, i + 1) is calculated using the unit vector nj of the respective bond. Note that this calculation has to be unique and invertible for an interval of 2π, which demands for a case distinction when α is computed. As for the variations δγi = δαi+1 − δαi ,
(5.17)
∆δγi = ∆δαi+1 − ∆δαi , we can define αi by αi = arctan
(5.18)
yi − yi−1 ni,y = arctan , xi − xi−1 ni,x
(5.19)
with ni,x and ni,y being the x and y components of vector ni , respectively. The variation δαi becomes δαi = 1+
where M is the rotation matrix for an angle of π/2 [rad]. From a differential geometry point of view, we may also replace δni by δai /k ai k, which can be computed more efficiently (see Fig. 5.4).
The increment of δα is more complex to derive and is given by ∆δαi = −δaTi
ni nTi MT + M ni nTi p ∆ai . aTi ai
(5.21)
With the obtained relations, we can set up the element stiffness matrix Ke as well as the element vector fint,e by backward substitution of Eqs. (5.12) ÷ (5.21) into Eqs. (5.9) ÷ (5.11). The full expressions for the resulting tensors can be found on the following page.
∆t δΠint = −t δΠint − t+∆t δΠext X E Area 1 EI [∆γi δγi + (γi − Γi ) ∆δγi ] [∆ai δai + (ai − Ai ) ∆δai + ∆ai+1 δai+1 + (ai+1 − Ai ) ∆δai+1 ] + Ai 2 Ai i q E A T T T xTi BTI BI xi − Ai X BI xi xTi BTI rea 1 BI BI xi xi BI BI T T BI + q BI I − T T δxi = Ai 2 xTi BTI BI xi xi BI BI xi xT BT B x i
∆δΠint =
∆δΠint
i
I
I
(5.22)
i
q T T T T T xTi BTII BII xi − Ai B x x B BII BII xi xi BII BII II i BTII I − T T i II BII + + q xTi BTII BII xi xi BII BII xi T T xi BII BII xi T T T T xi BII M BII BTI MBI xi xTi BTI MT BI EI BII MBII xi − (αi+1 − αi − Γi ) − T T − T T + Ai xTi BTII BII xi xi BI BI xi xTi BTII BII xi xi BI BI xi !#) BTII BII xi xTi BTII MT + MBII xi xTi BTII BII BTI BI xi xTi BTI MT + MBI xi xTi BTI BI · ∆xi − xTi BTII BII xi xTi BTI BI xi X δxTi Ke ∆xi =
(5.23)
(5.24)
i
E Area 1 EI (γi − Γi ) δγi − δxTi fext [(ai − Ai ) δai + (ai+1 − Ai ) δai+1 ] + Ai 2 Ai i q q T BT B x − A T BT B x − A E A 1 x x X I i i II i i i i I II rea q −δxTi = BTI BI xi + q BTII BII xi Ai 2 T T T T x B B x x B B x i
X − − t δΠint + t+∆t δΠext =
i
EI (αi+1 − αi − Γi ) + Ai
=
X
δxTi fe
I
I
i
i
BTI MBI xi BTII MBII xi − xTi BTII BII xi xTi BTI BI xi
II
II
− fext
i
(5.25)
(5.26)
(5.27)
i
28
29 These expressions are now assembled to obtain a linear set of equations K(i) ∆x(i) = (i)
fext − fint , which has to be repeatedly solved according to Sec. 5.1. Note that K(i) and (i)
fint have to be re-evaluated for every iteration step.
The exemplary implementation in C language can be found in Appendix B.3. Results will be discussed later in Sec. 5.4, to be able to compare with the following finite element formulation.
ey
δni δαi
ni
αi
ez
ex
Figure 5.4.: δni = δαi ez × ni
30
5.3. Finite element formulation Y, y X, x
A
Figure 5.5.: Reference configuration for the exponential Cauchy-Born rule To establish a finite element (FE) formulation, we need to derive an expression for the strain energy densityk first. For this purpose, we have to determine the components of the exponential Cauchy-Born rule, according to Chapter 4. We start by assigning an initial (undeformed) equilibrium configuration to the chain. It is defined to be equal to the chosen 1-dimensional reference configuration, which is depicted in Fig. 5.5. The first procedure of the exponential Cauchy-Born rule is the (“inverse exponential”) mapping exp−1 X of bond vectors from the initial configuration into its tangent space. With the present initial state, this mapping results in an identity map which can be omitted. Next, the deformation gradient F has to be applied. Note that the deformation map Φ (X)∗∗ is a vector valued operator, that is exclusively depending on the reference coordinate X. Thus the deformation gradient reduces to be a vector F =
dΦx dX
dΦy dX
T
,
and the right Cauchy-Green deformation tensor C = FT F becomes a scalar C. Using F, we can map lattice vector A from the undeformed tangent space into the deformed one to obtain w, see Fig. 4.3. The resulting length of vector w can be easily calculated by
where
√
w = k FA k =
k
CA ,
(5.28)
C is identified to be the stretch √
∗∗
√
C=
s
dΦx dX
2
+
dΦy dX
2
.
in dependence of the (local) continuum measures strain and curvature for this special reference configuration
(5.29)
31
r γ
a
w
Figure 5.6.: Osculating circle To complete the “exponential Cauchy-Born” mapping, we need to find the exponential map expΦ(X) , which maps the tangent vector w into the deformed configuration a. For this mapping, we have to “roll” w onto the osculating circle (Fig. 5.6), to determine the resulting chord (bond) length a as well as the angle γ formed with neighbouring chords. The radius r of this circle can be derived as the reciprocal of the curvature κ: 1 r= κ
κ=−
1 3
C2
dΦy d2 Φx dΦx d2 Φy − dX dX2 dX dX2
(5.30)
The triangle formed by the center of the circle and the chord is isosceles. The angle found at the center of the circle can be identified to be the angle γ enclosed by neighbouring bonds, and it is given by γ=
√ w = κw = κ CA . r
(5.31)
The bond length a is represented by the unequal side of the triangle, and is a = k a k = 2r sin
2 γ γ = sin . 2 κ 2
(5.32)
With expressions for a and γ, we can now formulate a potential for the deformed bond in √ terms of the continuum deformation measures stretch ( C) and curvature (κ). Accordingly, the elastic energy density†† W can be found as the energy for one bond, divided by its undeformed length A (see Eq. (5.1)) " # 1 E Area (a − A)2 E I γ 2 . + W (C, κ) = A 2 A 2 A ††
per undeformed volume
(5.33)
32
The internal potential can now be found as the energy density W , integrated over the undeformed configuration Ω0 Πint =
Z
W dΩ0 .
(5.34)
Ω0
With the obtained expression, we are now able to discretize atomistic systems using Hermitean finite elements‡‡ . Similar to section 5.2, we need to define the components of Eq. (5.4) ∆t δΠint = −t δΠint − t+∆t δΠext . This time, the elements are evaluated to δΠint ∆δΠint
Z
E Area EI = (a − A) δa + 2 γ δγ dΩ0 , A2 A Ω0 Z EI E Area [∆a δa + (a − A) ∆δa] + 2 (∆γ δγ + γ ∆δγ) dΩ0 , = A2 A
(5.35) (5.36)
Ω0
δΠext = −
Z
δxT fext dΩ0 .
(5.37)
Ω0
Expressions for the variations of a and γ can be derived from Eq. (5.32) and (5.31): cos γ2 γ 2 δa = − δκ 2 sin + δγ κ 2 κ
(5.38) γ
cos 2 4 γ 2 γ sin ∆κ − ∆δκ 2 sin − δκ ∆γ 3 κ 2 κ 2 κ2 sin γ2 cos γ2 cos γ2 ∆κ − δγ ∆γ + ∆δγ − δγ κ2 2κ κ √ 1 1 δγ = δκ CA + δC κ C − 2 A 2 √ 1 1 ∆δγ = ∆δκ CA + δκ C − 2 A ∆C 2 3 1 1 −1 κ κ + δC C 2 A ∆κ − δC C − 2 A ∆C + ∆δC C − 2 A 2 4 2 ∆δa = δκ
‡‡
Bounded second derivatives are necessary, because of the curvature κ.
(5.39)
(5.40)
(5.41)
33
y
2 x2 , y2 , φ2 ,
1 x1 , y1 , φ1 ,
ds dξ 1
x
Y
ds dξ 2
0
≤
ξ
≤
2
1
0 ≤ s ≤ 1
1
L X 2
Figure 5.7.: Hermitean element To be able to evaluate the continuum deformation measures stretch
√
C and curvature
κ as well as their variations, we need the Hermitean interpolations for the position coordinates first. They are dx dx x (ξ) =N1 x1 + N2 + N3 x2 + N4 , dξ 1 dξ 2 dy dy y (ξ) =N1 y1 + N2 + N3 y2 + N4 , dξ 1 dξ 2
(5.42) (5.43)
with N1 =1 − 3ξ 2 + 2ξ 3 , N2 = ξ − 2ξ 2 + ξ 3 , N3 = 3ξ 2 − 2ξ 3 , N4 = −ξ 2 + ξ 3, (5.44) where ξ is the natural coordinate that 0 to 1. For a neat formulation, runs from we dy dx dx ds ds replace the positional derivatives dξ by ds dξ = cos φ1,2 dξ and dξ by 1,2 1,2 1,2 1,2 1,2 dy ds ds ds dξ = sin φ1,2 dξ , where φ1,2 are the angles of the tangents counted from 1,2
1,2
1,2
the abscissa x, and s is the arc length.
34
The formulations for the Cauchy-Green scalar C (Eq. (5.29)) and curvature κ (Eq. (5.30)) are " 2 2 2 # 2 dx dx dξ dy dξ 1 dy , (5.45) + = + = 2 + dξ dX dξ dX L dξ dξ 1 dx d2 y 1 dy d2 x dy d2 x dx d2 y =− 3 , (5.46) − − κ=− 3 C 2 dX dX2 dX dX2 C 2 L3 dξ dξ 2 dξ dξ 2
dx C= dX
2
dy dX
2
where L is the length of a finite element in its reference configuration∗ . The variations of C and κ are given by 2 dx dx dy dy δC = 2 δ , +δ L dξ dξ dξ dξ 2 dx dy dy dy dy dx dx dx ∆δC = 2 δ ∆ + , ∆δ +δ ∆ + ∆δ L dξ dξ dξ dξ dξ dξ dξ dξ 3κ d2 x dy dx d2 y d2 y dx dy d2 x + δC +δ 2 −δ −δ 2 , δ δκ = − 3 2 2 3 dξ dξ dξ dξ dξ dξ dξ dξ 2C C2L 1
∆δκ = −
1 3
C 2 L3
(5.48)
(5.49)
dy d2 x dy d2 x d2 x dy d2 x dy ∆δ + δ ∆ 2 + ∆δ 2 + δ 2 ∆ dξ 2 dξ dξ dξ dξ dξ dξ dξ
d2 y dx dx d2 y dx d2 y d2 y dx − 2 ∆δ − δ ∆ 2 − ∆δ 2 − δ 2 ∆ dξ dξ dξ dξ dξ dξ dξ dξ dy d2 x 3 1 d2 x dy dx d2 y d2 y dx ∆C + δ +δ 2 −δ −δ 2 dξ dξ 2 dξ dξ dξ dξ 2 dξ dξ 2 C 25 L3 − δC
(5.47)
(5.50)
3 κ 3κ 31 ∆C + δC ∆κ + ∆δC , 2 2C 2C 2C
As we can see, we will need to determine the first and second derivatives of the Hermitean interpolations for x and y (with respect to ξ), as well as their variations. Respective expressions can be found in Appendix A.
∗
Note: In general, L is not equal to the undeformed / reference bond length A.
35 (i)
Finally, we can build a linear system of equations K(i) ∆x(i) = fext − fint by assembling
Eq. (5.4). This is done in the following way:
First, we are substituting the shape functions† into Eqs. (5.45) ÷ (5.50) to obtain the
continuum measures C and κ as well as their variations δ and increments ∆ in terms of the element’s degrees of freedom (xi )‡ and the natural coordinate ξ. The resulting
expressions are again substituted into Eqs. (5.38) ÷ (5.41) (a, γ, respective variations and increments), which are further substituted into Eqs. (5.35) ÷ (5.37) (internal and
external potential expressions).
The necessary integration over the natural coordinate ξ is carried out using a Gauss quadrature. Note that shape functions are contained within trigonometric functions§ , which makes an exact quadrature impossible. In the present work we employ 8 Gauss points for standard cases, to obtain high precision results. Note that further remarks concerning this topic will be made in the following section. We have now obtained the desired set of equations for arbitrary Hermitean elements, (i)
which are assembled to obtain the global system K(i) ∆x(i) = fext − fint . The numerical treatment for solving this system is equivalent to the procedure for the atomistic system (i)
(see Sec. 5.2): The system has to be repeatedly solved for ∆x(i) , and K(i) and fint have to be re-evaluated for every iteration step i. The implementation in C language can be found in Appendix B.4. Because quadratic matrices for the FE system are not necessarily positive definite, and the chosen benchmark systems are relatively small, we employ a LU decomposition¶ to solve the resulting set of linearized equations.
†
in their respective tensor forms; see Appendix A iT h ds ds ‡ x y1 φ1 dξ x2 y2 φ2 dξ 1 1 2 § See e.g. Eq. (5.38) ¶ according to [8]
36
5.4. Results and discussion
4000
y @mmD
3000
2000
1000
0 -1000
0
1000
2000 3000 x @mmD
4000
5000
6000
Figure 5.8.: Rolling of a straight FE - chain into a circle; 64 finite elements The following section illustrates benchmark simulations that have been carried out to test the developed programs. For all given figures, a red chain indicates a multibody simulation, where diamonds represent the atoms. All given multibody systems consist of 64 bonds / atoms, and parameters are set to the values that are given in Tab. 5.1. Finite element simulations are usually depicted in blue, where boxes represent node points. Discretizations for the finite element simulations are identical with those for the multibody systems. Figure 5.8 shows a finite element simulation of the rolling of a straight chaink into a circle. Employed constraints for this calculation are fixed coordinates x, y and φ for the node located at (0,0). The applied load is a progressively increased moment at the end of the chain. The number of load increments is 24. This benchmark demonstrates, that the equivalent continuum potential (5.33) employing the exponential Cauchy-Born rule works. As we have discussed in Chapter 4, the standard Cauchy-Born rule would not assign any bending stiffness for the given initial configuration, which would render the model numerically instable. A multibody reference to test the precision of the result cannot be given for this case, because we cannot apply a moment onto the atomistic system. k
in its reference configuration
37
2000
2000
1500
1500
y @mmD
y @mmD
1000
1000
500
500 0
0 -500
-1000
-500
0 x @mmD
500
1000
1500
(a)
-1000
-500
0 x @mmD
500
1000
1500
(b)
Figure 5.9.: Comparison of FE and multibody simulation; 64 finite elements; (a) shows a deformation driven case, while (b) is force driven. Comparisons to test the accuracy of the finite element model are shown in Fig. 5.9. Fig. 5.9(a) gives a deformation driven case. Constraints for both the multibody and the finite element systems are fixed coordinates at the circle’s floor (0,0), while the roof (0,2000) is succeedingly shifted downwards. Fig. 5.9(b) depicts a force driven case, where the floor (0,0) is constrained (x = y = φ = 0), and the load is a uniformly distributed weight. The deformation driven case (Fig. 5.9(a)) shows a very good correlation for the multibody and the finite element systems. Variations of the coordinates are hard to notice. The force driven case (Fig. 5.9(b)) shows some minor discrepancies. The finite element discretization with 64 elements ends up slightly stiffer than the multibody system. A very interesting feature appears, when the same system is discretized using only 16 finite elements (Fig. 5.10): The finite element model (indicated by triangles and boxes) experiences considerable stiffening, especially when the chain is severely bent in local regions. This behaviour is very similar to locking phenomena observed with continuum beam models . In these cases, a common approach is the application of reduced integration. Figure 5.10 demonstrates the benefits of a cutback in Gauss points: A reduction from 14 (boxes) to 3 Gauss points (triangles) improves the positional error by approximately 30% compared to the initial error.
38
2000
1500
y @mmD
1000
500
0
-500
-1000
-500
0 x @mmD
500
1000
1500
Figure 5.10.: Locking phenomenon and benefit of reduced integration; blue: 14 Gauss points per element; black: 3 Gauss points per element.
2000
1500
y @mmD
1000
500
0
-500
-1000
-500
0 x @mmD
500
1000
1500
Figure 5.11.: Reduced locking for the case of reduced axial (bond) stiffness; blue: 14 Gauss points per element; black: 3 Gauss points per element.
39
Figure 5.12.: Comparison of atomistic (left) and continuum model (right). As it is evident that we are dealing with a certain kind of locking, it is not yet clear where to locate the root of our problems. Arroyo’s proposed model [1] distributes the energy of the atomistic (multibody) system into a continuum, which usually results in softer systems than the original (discrete) ones. If we illustrate the present model, we can find an explanation for the aforementioned behaviour: Think of the given benchmark case. 16 finite elements simulating 64 bonds makes one finite element for every 4 bonds. Using 14 Gauss points for one finite element, we model 14 bonds for this chain section to obtain the strain energy density for every point. Doubtless, we are weighting these energies accordingly, but note that we do not model a physically equivalent one, but only a similar system. The discussed relations are shown in figure 5.12. Considering this figure, we can see the following: As the bonds are overlapping, the system resembles a framework made of atomic bonds∗∗ . The structure of the framework displays a certain spatial thickness, especially for the case of severe curvature. This architecture resembles that of bones, where bending forces are conducted using the extensional stiffness of very small “trusses”. The following conclusion can be drawn: Our model transfers bending energy into axial strain energy†† . With the present material (potential-) parameters, the model is very stiff with regard to stretching, but weak with regard to bending. Axial (bond-) forces arise that are bending out the chain and stiffen the present model.
∗∗ ††
In our case, the bonds are not unconstrained, but coupled via the Hermitean shape functions. Shape functions are coupling bending and stretching potentials, which are physically uncoupled.
40
The presented theory is confirmed by figure 5.11. This benchmark is fully identical to the previous one (Fig. 5.10), but axial stiffness‡‡ is reduced by a factor of 100∗ . Now, bending forces cannot be absorbed by axial ones to the same extent, and the FE model gains accuracy. It is necessary to mention that for the simulation of structural behaviour at the atomic scale, Arroyo’s model [1] is quite accurate. Note that the values of parameters for the potential energy are arbitrary ones in the present thesis, with bonds relatively stiff concerning axial elongation. For realistic interatomic potential formulations† , the relation of bending stiffness‡ to axial stiffness is almost reversed, thus relevant locking effects do not occur.
‡‡
In this case represented by parameter Area The axial stiffness for the multibody simulation (depicted by diamonds) has accordingly been reduced. † E.g. Brenners interatomic potential for carbon, see Sec. 2.2 ‡ ... of bonds one to another
∗
References [1] M. Arroyo, Finite Crystal Elasticity for Curved Single Layer Lattices: Applications to Carbon Nanotubes, Ph.D. thesis, Northwestern Universtity, 2003. [2] J. Bonet and R.D. Wood, Nonlinear continuum mechanics for finite element analysis, Cambridge University Press, 1997. [3] D.W. Brenner, Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films, Physical Review B 42 (1990), 9458 – 9471. [4] C. C. Celigoj, Finite Element Method, Institute for Strength of Materials, Graz University of Technology, 1998. [5] K. Chandraseker and S. Mukherjee, Coupling of extension and twist in single-walled carbon nanotubes, ASME Journal of Applied Mechanics 73 (2006), 315 – 326. [6] K. Chandraseker, S. Mukherjee, and Y.X. Mukherjee, Modifications to the CauchyBorn rule: Applications in the deformation of single walled carbon nanotubes, International Journal of Solids and Structures 43 (2006), 7128 – 7144. [7] H. Jiang, P. Zhang, B. Liu, Y. Huang, P.H. Geubelle, H. Gao, and K.C. Hwang, The effect of nanotube radius on the constitutive model for carbon nanotubes, Computational Materials Science 28 (2003), 429 – 442. [8] W. Press, S. Teukolsky, W. Vetterling, and B. Flannery, Numerical Recipes in C, The Art of Scientific Computing, Cambridge University Press, 1990. [9] K. Schmaranz, Softwareentwicklung in C, Springer, 2001. [10] J. R. Shewchuck, An Introduction to the Conjugate Gradient Method Without the Agonizing Pain, www.cs.berkeley.edu/∼jrs/jrspapers.html, 1994.
A. Hermitean element: Interpolations, respective derivatives and their variations The fundamental shape functions for the Hermitean element are (Eqs. (5.42) ÷ (5.44)) ds ds + N3 x2 + N4 cos φ2 , x (ξ) =N1 x1 + N2 cos φ1 dξ 1 dξ 2 ds ds y (ξ) =N1 y1 + N2 sin φ1 + N3 y2 + N4 sin φ2 , dξ 1 dξ 2 with N1 =1 − 3ξ 2 + 2ξ 3 , N2 = ξ − 2ξ 2 + ξ 3 , N3 = 3ξ 2 − 2ξ 3, N4 = −ξ 2 + ξ 3 . The first derivatives with respect to ξ are ds ds dx + N3,ξ x2 + N4,ξ cos φ2 , =N1,ξ x1 + N2,ξ cos φ1 dξ dξ 1 dξ 2 dy ds ds + N3,ξ y2 + N4,ξ sin φ2 , =N1,ξ y1 + N2,ξ sin φ1 dξ dξ 1 dξ 2
Expressions for the second derivatives are just the same, but using the second derivatives of shape functions (N,ξξ ) instead of the first ones. The C code given in Appendix B.4 makes use of the following notation: d2 x δ 2 = δxT Rx , dξ d2 x ∆δ 2 = δxT Sx ∆x , dξ
d2 y δ 2 = δxT Ry , dξ d2 y ∆δ 2 = δxT Sy ∆x . dξ
(A.13) (A.14)
B. Source code listings B.1. Equilibrium positions on a SWNT 1 2 3
// nanotube . c // C a l c u l a t e s atom c o o r d i n a t e s , d i a m eter , bonding l e n g t h s and a n g l e s f o r // s i n g l e w a l l e d ca r b o n n a n o tu b es (SWNTs)
