## Abstract

In this work, we present a computational study of the small strain mechanics of freestanding ultrathin carbon nanotube (CNT) films under in-plane loading. The numerical modeling of the mechanics of representatively large specimens with realistic micro- and nanostructure is presented. Our simulations utilize the scalable implementation of the mesoscopic distinct element method of the waLBerla multi-physics framework. Within our modeling approach, CNTs are represented as chains of interacting rigid segments. Neighboring segments in the chain are connected with elastic bonds, resolving tension, bending, shear, and torsional deformations. These bonds represent a covalent bonding within the CNT surface and utilize enhanced vector model (EVM) formalism. Segments of the neighboring CNTs interact with realistic coarse-grained anisotropic van der Waals potential, enabling a relative slip of CNTs in contact. The advanced simulation technique allowed us to gain useful insights on the behavior of CNT materials. It was established that the energy dissipation during CNT sliding leads to extended load transfer that conditions size-independent, material-like mechanical response of the weakly bonded assemblies of CNTs.

## 1 Introduction

Ultrathin freestanding carbon nanotube (CNT) films [1] possess a number of unique physical and mechanical properties, making them attractive for a number of interesting applications, in particular, in flexible electronics: transparent displays [2], pellicles for ultraviolet and soft X-ray radiation [3], wearable medical devices [4], etc. The major advantage of these materials, as compared with possible competitors, is their ability to withstand significant plastic deformations without irreversible structural changes. This property is conditioned by the discontinuous nature of the CNT films—they consist of individual CNTs that are chemically reactive and bound only by van der Waals (vdW) forces. Due to the molecular structure of CNTs, the interfaces between neighboring parallel CNTs resolve relative sliding of one CNT along the other, with very small static friction. Such sliding provides a mechanism for the plastic flow of the CNT material. This mechanism is non-destructive in a sense that small strain plastic deformations do not lead to significant changes in the microstructure of a CNT material, as well as in its mechanical, thermal, optical, and electric properties.

An interesting question to ask is whether the deformation localizes or remain distributed in freestanding CNT systems. The experimental observations indicate that the deformation remains distributed—one can apply significant (few percent) strains to the centimeter-sized ultrathin freestanding CNT films [5], implying that the overall specimen elongation is many orders of magnitude larger than the length of individual CNTs. Such a behavior implies significant load transfer between CNTs; however, this load transfer is hard to explain, given ultralow sliding friction in CNT systems (see, e.g., Ref. [6]).

The modern capabilities for experimental research on the properties of CNT materials are limited both in terms of technology and cost-efficiency. Thus, it is highly desirable to have a numerical model that could reliably reproduce the mechanics of thin CNT films. Although molecular-level modeling techniques are capable to provide useful insights on the properties of individual nanotubes [710], they are inefficient for modeling assemblies comprising large numbers of CNTs—the most extreme molecular dynamics calculations so far included only hundreds of CNTs [11].

In order to bypass the computational inefficiency of full molecular modeling techniques, a number of coarse-grained mesoscale models were developed [1218]. The general idea of these models is to reduce the number of model degrees of freedom in order to reach the computationally tractable model size. However, such coarse-graining often leads to significant artifacts in the model’s behavior. For example, pair vdW potential used in Refs. [12,13] leads to the absence of the CNT relative sliding, which makes the model inapplicable for modeling irreversible deformations of CNT assemblies. Similarly, giving up on torsional degrees of freedom in models [15,16] significantly alters the mechanical properties of CNT assemblies under certain types of loading. Moreover, neither model so far provides a realistic description of energy dissipation during relative sliding of CNTs, which makes it impossible to describe rate-dependent deformation properties of CNT materials.

In our works [1921], we have developed the mesoscopic modeling technique based on the distinct element method [22] (DEM). Distinct element formulation is well suited for the analysis of deformations in discontinuous systems, taking into account complex nonlinear interaction laws between structure parts. Within our approach, each CNT is represented as a chain of rigid bodies possessing tensile, bending, shear, and torsional stiffnesses. These stiffnesses are calibrated based on molecular dynamics simulations. CNTs interact via vdW potential tailored to provide an adequate description of relative CNT sliding. In our thermal model, energy dissipation was represented with viscous damping, appropriately calibrated in a top-down fashion [20]. The model showed to be efficient for the description of CNT self-folded configurations [23,24], weakly bonded and cross-linked CNT ropes and bundles [20,21], and qualitative description of large-strain deformations of CNT films [25]. Until recently, the model was lacking a scalable parallel implementation, which severely limited its capabilities in terms of achievable length and time scales of the simulation. We addressed this reservation with the development of a scalable parallel implementation of our model based on the rigid particle dynamics module of the waLBerla package [26].

In our present work, utilizing our new parallel computation capabilities, we investigate the mechanics of representatively large assemblies of thin films. We present the study of the representative volume element (RVE) size, dependence of the mechanics of CNT films on dissipation rate, film thickness, and length of individual CNTs. It appears that the presence of dissipation during relative CNT sliding provides the mechanism for a load transfer, leading to size-independent, material-like response of CNT materials.

## 2 Method

Our mesoscopic model mostly corresponds to the ones employed in our earlier works [19,20,26]. For clarity, we present a brief description of our modeling technique here. The model is based on the mesoscopic distinct element method that computes the damped dynamics of a collection of interacting rigid bodies of given mass, shape, and tensor of inertia—distinct elements. The state variables for each element include translational positions and velocities as well as rotational positions and angular velocities. The bodies change their velocities and angular velocities due to contact forces and moments arising in pair interactions as well as external forces and moments acting at each body. The system is evolved in time with an explicit velocity Verlet time integration scheme. The translational degrees of freedom is given as
$x(t+Δt)=x(t)+v(t)Δt+12a(t)Δt2v(t+Δt)=v(t)+a(t)+a(t+Δt)2Δt$
(1)
with x being the position, v the velocity, a the acceleration, t the current time, and Δt the time increment per time-step. Evolution of rotational degrees of freedom is analogous, and the rotations are stored as quaternions (see Ref. [27] for a detailed description).

Within our approach, undeformed CNTs are homogenized into cylindrical shells with finite thickness and then partitioned into identical segments of finite lengths T (Fig. 1). Each distinct element represents the inertial properties of a CNT segment. Table 1 provides segment parameters for (10, 10) CNTs with diameter 2rCNT = 13.560 Å and length T = 2rCNT. Each segment contains approximately 220 carbon atoms. The properties of such segments were parameterized in our previous works [19,20].

Table 1

Parameterization of the spherical particles and EVM bonds for a (10, 10) CNTs. m, r, I are the mass, radius, and moment of inertia of each spherical particle. B1, B2, B3, and B4 are EVM stiffnesses.

m (amu)r (Å)I (amu × Å2)B1 (eV/Å2)B2 (eV)B3 (eV)B4 (eV)
264910.721.218 × 10567.5919780−40321471
m (amu)r (Å)I (amu × Å2)B1 (eV/Å2)B2 (eV)B3 (eV)B4 (eV)
264910.721.218 × 10567.5919780−40321471
Fig. 1
Fig. 1
Close modal
The inertial properties of a segment are represented with mass m and the spherical moment of inertia $I=25mr2$. We utilize spherical tensors of inertia for simplicity, which is sufficiently a good approximation for the case of cylindrical segments with an aspect ratio of 1. Parameters m and I are matched with the mass and moment of inertia of a cylindrical segment taken with respect to the CNT axis by a proper choice of the radius of a distinct element:
$r=2.5rCNT$
(2)
The elements representing CNT segments are equispaced at a distance T.

Such a coarse-graining leads to the elimination of atomistic degrees of freedom. Therefore, the dissipative microscopic processes associated with CNT sliding are not explicitly captured and should be included in a phenomenological manner. In this respect, two independent channels of energy dissipation are introduced in our model, local damping and viscous damping, controlled by the parameters α and β, respectively.

PFC-style local damping [28] acts at each body. It is introduced to damp stiff interactions and stabilize time integration in dense particle assemblies. The components of the damping force $Fiα$ (moment $Miα$) are proportional to the corresponding components of the unbalanced force Fi (moment Mi) according to
$Fiα=−α|Fi|sign(vi),Miα=−α|Mi|sign(ωi)$
(3)
Here, vi and ωi are components of the translational and rotational velocity of an element and sign(x) is the sign function. In our simulations, the local damping is kept relatively low (α = 0.1) and its influence on CNT dynamics is insignificant.
Recent theoretical findings [29] predict a linear dependence between force and velocity in nanoscale thermal-induced friction. In our coarse-grained model, we introduce the viscous damping that provides a simplified model for nanotribological energy losses during CNT sliding. The viscous damping forces, proportional to relative segment velocities, act in parallel with vdW contact forces. These forces are controlled by parameter β. Normal $Fnβ$ and tangential $Fsβ$ viscous forces are proportional to normal vn and tangential vs relative velocities of elements in the vdW contact
$Fnβ=cnvn,Fsβ=csvs$
(4)
Viscosity coefficients cn and cs are related to β as follows:
$cn=2βmkn,cs=2βmks$
(5)
where kn and ks are linearized stiffnesses of the contact model, taken kn = 1 eV/Å2 and ks = 1 eV/Å2.

Damping factor β is kept as a free parameter in our simulations. The initial estimate of β is performed in a top-down way, as discussed in Ref. [20]. Unless otherwise noted β = 0.1 is assumed in this work.

The elasticity of CNTs in our model is represented with the formalism of the enhanced vector model (EVM) [30]. Unlike incremental parallel bonds, widely used in DEM community [31], EVM bonds explicitly conserve energy in symplectic time integration. EVM is based on a binding potential, describing the behavior of an elastic bond linking two arbitrary-shaped rigid bodies. The formulation provides straightforward generalization for the case of large strains and accounts for a bending-twisting coupling.

Consider two equal-sized particles i and j with equilibrium separation T and orientation described in terms of orthogonal vectors nik, as depicted in Fig. 2(a). Then, the EVM bond potential is given as follows:
$U(rij,nik,njk)=B12(rij−T)2+B22((nj1−ni1)rij/rij+2)+B3ni1nj1−B42(ni2nj2−ni3nj3)$
(6)
Here, rij is the radius vector connecting centers of bonded particles. Parameters B1B4 are directly related to longitudinal, shear, bending, and torsional rigidities of a bond, according to the Euler–Bernoulli beam theory (see Ref. [30] for more details):
$B1=ESTB2=12EJTB3=−2EJT−GJp2TB4=GJpT$
(7)
where S, I, and J are area, moment of inertia, and polar moment of inertia of a cylinder shell beam with radius rCNT and thickness h:
$S=2πhrfJ=πhrf(rf2+h2/4)Jp=2J$
(8)
Fig. 2
Fig. 2
Close modal

Table 1 summarizes the elastic properties of CNT segments.

The potential (6) with calibrations (7) reproduces the elastic behavior of a CNT in four admissible deformation modes (Fig. 2(b)).

Van der Waals interactions between CNT segments are represented with the coarse-grained anisotropic vdW potential as described in Refs. [19,20]. The potential is defined in terms of three variables R, θ, and ϕ, related to mutual position and orientation of two cylindrical segments (Fig. 2(c)):
$U(rij,θ,ϕ)=fc(rij)Vk(rij,θ)Φ(rij,ϕ)Vk(rij,θ)=ϵ′(A′(Dk(rij,θ))α′−B′(Dk(rij,θ))β′)Dk(rij,θ)=rijrCNTΘk(θ)−2Φ(rij,ϕ)=1+Wϕ(rij)(1−cos(2ϕ))Θk(θ)=1+∑i=1kCi((−1)i−1+cos(2iθ))Wϕ(rij)=Cϕ(rij/rCNT)δfc(rij)=∑i=03Qi(rij/8rCNT)i$
(9)

The potential describes interactions between CNT segments, taking into account their relative orientation, and providing symmetry for the interaction of two parallel infinite straight CNTs with respect to translation of one CNT along its axis. This property is crucially important for the correct representation of shear interactions in CNT bundles and load transfer in CNT materials. It is important to note that such vdW potential does not create any barriers to CNT sliding, which corresponds to the incommensurate mode of CNT slip. Such a model does not create any static friction between CNTs and therefore will underestimate the strength of CNT materials in comparison with the full atomistic model, where commensurate contacts between neighboring CNTs with nonzero static friction are possible.

The parameters of the anisotropic vdW potential for (10, 10) CNTs are provided in Table 2.

Table 2

Parameterization of vdW potential for (10, 10) CNTs

ε′, meVαβABδCϕ
149.39.540.02231.31−7.590
kC1C2C3C4C5
50.358180.03263−0.00138−0.000170.00024
Q1Q2Q3Q4
−80.0288.0−336.0128
ε′, meVαβABδCϕ
149.39.540.02231.31−7.590
kC1C2C3C4C5
50.358180.03263−0.00138−0.000170.00024
Q1Q2Q3Q4
−80.0288.0−336.0128

Large-scale simulations of CNT assemblies are possible using the waLBerla framework. The development of the framework is focused on providing a highly parallel and highly optimized basis for multi-physics applications [32]. It offers a rigid particle dynamics module that is capable of simulating up to 2.8 × 1010 particles with up to 1.8 × 106 parallel processes [27]. A detailed description of the algorithms and their implementation used in the framework can be found in Refs. [27,33,34]. The framework is freely available under GPLv3 license.2 Here, we only outline the basic features used for our simulations.

The simplified course of our simulation is presented in Fig. 3. The simulation begins with the generation of the initial geometry of CNTs and imposition of boundary conditions. Next, the simulation domain is partitioned using the distributed forest of octrees approach implemented in the waLBerla framework. Since no refinement is used, all the subdomains have exactly the same rectangular shape and are arranged on a regular three-dimensional grid. Subsequently, the subdomains are distributed among the available processes in a balanced manner. Every process is now responsible for one or more subdomains. In the following, we will denote the numbers of subdomains in the x, y, and z directions with a set of three integers (M, N, K). During the simulation, each process stores the information of all rigid bodies that are within its associated subdomains. At the next stage, time integration cycles are performed on all MPI processes. The integration cycle consists of four parts. First, all particle contacts are detected. To reduce the number of expensive contact checks, the advanced algorithms based on hierarchical hash grids and bounding volumes are used. In our implementation, the contact detection scheme used for the rigid body dynamics is adapted for potential-based interactions. This is done by decoupling the particle’s inertial radius, used to compute its moment of inertia, and its interaction radius, used in contact detection schemes. The latter in our case is equal to the potential interaction the cutoff radius. After all contacts are determined, the contact models described above are used to calculate the forces and moments acting at each contact. In the next step, the forces and moments are used to calculate accelerations and angular accelerations for each particle. Subsequently, the particles’ positions and rotations as well as its velocities and angular velocities are updated by a velocity Verlet time integration scheme. Since the whole simulation is conducted in parallel, information that is relevant for more than one processes has to be updated. Therefore, the final step of the integration cycle is the synchronization. It not only takes care of synchronizing contact lists but also migrates particles between processes if the subdomain they belong to has changed. Throughout the whole simulation, the information about the current particle configuration is gathered over regular time intervals and saved to log files for later analysis.

Fig. 3
Fig. 3
Close modal

## 3 Numerical Results

### 3.1 Simulations Setup.

Here, we present the results of our numerical study on the mechanics of thin CNT films. Every test consists of two stages: (i) self-assembly of the CNT film specimen and (ii) mechanical test on the specimen. Figure 4 illustrates the geometry and boundary conditions used in our simulations. At self-assembly stage, initially separate, straight CNTs are deposited into a cuboid of the sizes l × l × h with triple periodic boundary conditions (Fig. 4(a)) and then evolved to an equilibrated state. The computational domain is decomposed onto MPI blocks in (M, N, 2) manner. After equilibration, the obtained configuration is gathered to a master process, saved, and then used in mechanical test simulations.

Fig. 4
Fig. 4
Close modal

The mechanical test is performed in the following few steps.

First, we resize the computational domain to the size sufficient to fully contain the stretched specimen. For example, in uniaxial tension test, we resize x domain borders to the size (− (1 + ɛ)lm; (1 + ɛ)l + m), where ɛ is the tensile strain and m is the small margin.

Second, we replace the periodic boundary conditions with destructive boundary conditions in directions of stretching (x-direction in uniaxial test and x and y directions in the biaxial test).

Then, we specify the new static subdivision into MPI blocks. For uniaxial tests, we subdivide the domain into long blocks (1, N, 2) oriented along the direction of stretching. For biaxial tests, we save the decomposition (M, N, 2).

Next, we load the equilibrated configuration into a new MPI block configuration and scatter it over available MPI processes.

After that, we specify grip regions and prescribe the velocities of grips. In a uniaxial tension test, the grips are thin cuboid slices of the specimen (Fig. 4(b)), which are moving in opposite directions. They move with the constant velocity, except small constant acceleration period at the beginning of the simulation, introduced to avoid the initial inertial peak in the response:
$v(x,t)={2xlv0⋅ex,t>t02xlv0t0t⋅ex,t≤t0$
(10)
In the case of biaxial tension test, the grips region is a margin surrounding the gage region, consisting of four grips (Fig. 4(c)). The prescribed velocity field of the grip elements in this case is given by
$v(x,y,t)={2xlv0⋅ex+2ylv0⋅ey,t>t02xlv0t0t⋅ex+2ylv0t0t⋅ey,t≤t0$
(11)

In the expressions above, ex and ey are unit vectors directed along x and y coordinate axes. In the simulations described below, the time t0 is set to be 500 timesteps (10 ps).

In the uniaxial test, the stress is measured as an arithmetic average of x projections of forces acting on two grips, related to a cross-sectional area hl. For the case of biaxial tests, four independently measured force projections are averaged to compute stress.

### 3.2 Self-Assembly of a Carbon Nanotube Film Specimen.

Figures 5 and 6 summarize the initial relaxation of a CNT film specimen. Figure 5 demonstrates the part of $2×2μm$ specimen at the beginning of the relaxation process (a) and after 106 cycles (20 ns) of relaxation (b) and (c). One can see that the initially straight CNTs are bent by vdW adhesive forces and most of the CNTs are joined into small bundles.

Fig. 5
Fig. 5
Close modal
Fig. 6
Fig. 6
Close modal

Figure 6 demonstrates the evolution of vdW adhesion energy and elastic strain energy during the relaxation. One can see that the energy balance depends significantly on the CNT length, whereas the size of the specimen almost does not affect it. The difference between random realizations appears to be negligible for large enough specimens.

The specimens generated for the mechanical tests described below had periodic out-of-plane boundary conditions, in-plane sizes from 0.5 × 0.5 μm to 3 × 3 μm, and CNT lengths of 100, 200, 400, and 800 CNT segments. The same volumetric density (0.044 g/cm3) was prescribed. Unless otherwise noted, the specimens have a thickness of 20 nm and the surface density of 0.088 μg/cm2.

### 3.3 Mechanical Tests on the Carbon Nanotube Film Specimen.

Here, we summarize the results of small strain mechanical tests, performed as described above. The mechanics of a CNT film is characterized by the elastic deformation of individual CNTs and a relative slip of the CNTs in the specimen. The resistance to slip is conditioned by the change in surface energy and dissipation due to damping, introduced in our model to represent the effects of energy transfer to implicit (microscopic) degrees of freedom. Both local and viscous damping cause energy dissipation that is directly proportional to the strain rate [19,20]. In our simulations, local damping is kept low and has a negligible effect on the mechanics of CNT assembly. Overall mechanical response of the film is characterized by the initial viscoelastic deformation at small strains, conditioned by deformation of individual CNTs, followed by visco-plastic region, characterized by relative sliding of CNTs in the film. Due to the viscous part of the response, the observed stress-strain curves start from the nonzero stress.

Figure 7 summarizes the mechanics of CNT films observed in uniaxial tests. The instantaneous Young’s modulus of CNT films observed in our simulations varies in a large diapason of 5–30 GPa and is conditioned by many factors. Figure 7(a) demonstrates the dependence of the response on strain rate in an undamped simulation that originates mostly in finite time necessary for rearrangements in the CNT assembly. At slow enough strain rates and in the absence of viscous damping, the response of the assembly is nearly zero. Increasing CNT length (Fig. 7(b)) improves the load transfer in CNT assembly and significantly slows down the stress relaxation. Figure 7(c) demonstrates the size dependence of the response in the presence of viscous damping β = 0.1. As will be illustrated below, this dependence is strongly affected by the viscous damping coefficient β. The dependence of the Young’s modulus on thickness is rather negligible due to out-of-plane periodic boundary conditions (Fig. 7(d)).

Fig. 7
Fig. 7
Close modal

The analysis of the velocity fields in CNT films during the mechanical test indicates that, given sufficient load transfer between CNTs, a CNT assembly deforms nearly uniformly, however, with some heterogeneity associated with the random microstructure of the film. Relative slip of CNTs can be observed at large magnifications (Fig. 8).

Fig. 8
Fig. 8
Close modal

It appears that the dependence of the response on the size of the specimen is always associated with localization of strain close to grip regions. At large strain rates and low intertube sliding dissipation, the specimen deformation appears to be nonuniform (Fig. 9(a) gives the velocity field in the case of nonuniform, size-dependent deformation). Conversely, at slow enough strain rates and large enough dissipation, the specimen deformation is uniform (Fig. 9(b)).

Fig. 9
Fig. 9
Close modal

Figures 10(a)10(c) demonstrate the convergence to size-independent response at slow strain rates. It gives the comparison of stress-strain curves, obtained for two specimens of sizes 1.5 × 1.5 μm and 3.0 × 3.0 μm. At snapshots (a)–(c), the simulation is carried out at 100%, 50%, and 25% strain rate, respectively. The viscous dissipation coefficient β is increased proportionally: it is 0.1, 0.2, and 0.4, respectively. One can see that at slow enough strain rates, the response becomes nearly size-independent which is associated with a negligible role of inertia-dependent CNT rearrangement processes. It is therefore clear that strain localization and size dependence of stress-strain curves are associated with inertial effects that are negligible in slow and highly damped simulations.

Fig. 10
Fig. 10
Close modal

We can therefore summarize the conditions for size-independent, material-like mechanical response of an assembly of CNTs and absence of deformation localization. It appears that such response occurs at strain rates that are low enough to exclude inertial effects ($ε˙<2×106$), dissipation coefficient sufficient to provide necessary load transfer between CNTs (β > 0.2) long enough CNTs (lCNT > 1 μm) and RVE size larger than the length of a CNT.

It is interesting to compare the behavior of CNT films in uniaxial and biaxial tests. Figure 10(d) gives an example of such a comparison. It appears that the slope of stress-strain curve is slightly higher in a biaxial mechanical test. Assume the validity of linear isotropic elastic relations for these initial slopes:
$Yuniax=E1−ν2Ybiax=E1−ν$
(12)

Figure 10(d) indicates the values of Yuniax = 17.5 GPa and Ybiax = 18.75 GPa. Therefore, it follows that the film instantaneously behaves as linear in-plane isotropic elastic material with the Young’s modulus of 17.4 GPa and Poisson’s ratio of 0.071. However, in simulations with long CNTs and low damping, where CNT sliding dominates the deformation, the observed Poisson’s ratio drops down to −0.7. The transitions between negative and positive Poisson’s ratios in thin CNT films depending on the specimen parameters have been observed experimentally [35].

## 4 Related Work and Discussion

The mechanics of CNT films was studied with the toolkit of mesoscale modeling by a large community of authors. The pioneering works in the field are due to Buehler and co-workers [12,13]. The use of bead-spring model allowed them to model large assemblies of CNTs and study stress-strain curves of a CNT material. Similar approaches were later used in a number of publications of other authors [36,37]. However, due to the artifact of the vdW potential corrugation, that is described in detail in Ref. [19], such models are not capable to adequately represent shear interaction between CNTs, which led these authors to a conclusion that the only available mechanisms of deformation involve elastic deformations of individual CNTs, along with “zipping and unzipping” of individual CNTs and bundles, resulting in histeretic, temperature-independent homogenized viscoelastic response [37]. The other group of authors [1518] combined symplectic finite temperature integration with “slippery” vdW potential which resulted in nearly zero load transfer between separate CNTs and absent quasi-static mechanical response of a CNT material. The model assumptions were largely based on two observations: (i) the CNT sliding observed in MD simulations does not feature significant viscous friction and (ii) in the absence of static and dynamic friction between CNTs, the model predicts high mobility of CNTs, leading to the formation of large bundles, similar to ones observed in real CNT films and papers. The first argument is challenged by the recent work [38], which indicates that the realistic MD simulation does indicate significant static and dynamic friction in graphitic systems. Regarding the second argument, it should be observed that the bundling of CNTs in real CNT film production occurs mostly on a gas/liquid stage and is driven by inelastic scattering and charge exchange processes, rather than by short-ranged vdW interactions. It should be also observed that the remaining rearrangements in CNT films, driven by vdW adhesion, usually take hours and days, whereas the bundles formation in mesoscale simulations [1518] take nanoseconds. Finally, it is clear that the real CNT films behave as a material, with the well-defined stress-strain curves, while frictionless simulations give zero quasi-static mechanical response.

In our work [20], we have tried to model dissipation-based load transfer between CNTs, which led us to a conclusion, that, in spite of the absence of noticeable static shear interactions between CNTs, such systems demonstrate material-like response. In this work, we conclude that the energy dissipation during CNT sliding does have the same effect in two-dimensional CNT systems.

It is worth noting that the material limit response has been observed recently for cross-linked systems of CNT films in FEM-based models [39,40]. However, mechanisms of size-dependence of the response are rather different in case of stiff elastic cross-linked assemblies and compliant discontinuous assemblies linked by vdW interactions. In the first case, size dependence comes from varying fraction of short (edge) CNTs in a specimen; in the second case, size dependence is dominated by the deformation localization, as discussed earlier.

Molecular dynamics simulations indicate that the incommensurate CNT tribology is well described by negligibly small static friction barrier and viscous energy dissipation during sliding. It is also important to note that the modern theoretical works [29] on nanofriction confirm the linear dependence between the friction force and the relative sliding velocity in the “thermal” regime of friction between nanosurfaces. However, precise calibration of such dissipation is problematic, since the observed viscosity coefficient depends on multiple factors and requires accurate statistical averaging. In this work, we leave intertube viscosity as a free parameter—the calibration is performed in a top-down way, based on the macroscopic stress-strain curves.

Accurate bottom-up calibration of energy dissipation in CNT sliding would require the design of the mesoscale model for such dissipation calibrated with accurate ab initio and molecular dynamics simulations. Such a model should also mimic static friction at commensurate CNT surface contacts.

## 5 Concluding Remarks and Future Work

In this work, we presented the computational study on the properties of thin CNT films under small strain deformations. Our study reveals the presence of a representative volume element in such structures and gives an idea of its size. Similarly to the earlier studied CNT rope structures [20], we observed an important role of energy dissipation that is responsible for even load distribution along the specimen and its uniform deformation, leading to a material-like mechanical response. The stress-strain curves observed in our simulations are qualitatively similar to the ones observed in real mechanical tests performed on CNT films [3].

A number of questions remained beyond the scope of the present study. We provided the description of the mechanics for a narrow class of such films produced by deposition of perfectly separate, length, and type purified CNTs on a flat substrate, which are assembled and ordered only by short-range vdW interactions. It is clear that realistic CNT assemblies might have a lot more complicated morphology, conditioned not only by vdW adhesion on the deposition stage but also by kinetic interactions and aggregation on the aerosol/suspension stage.

The top-down calibration of energy dissipation in the model allowed us to qualitatively reproduce the behavior of the realistic CNT film under mechanical loading. However, the most important drawback of our model in its current shape is the absence of an accurate model for the energy dissipation and its bottom-up calibration, which would allow us to reproduce the mechanics for CNT materials from the first principles. The development of such a model, taking into account static friction for different CNT chiralities and commensurate/incommensurate contacts, is the matter of future work.

## Acknowledgment

IO acknowledges Russian Science Foundation (Grant No. 17-73-10442; Funder ID: 10.13039/501100006769) for the support of the development of the parallel system of modeling carbon nanotubes based on the waLBerla framework. The financial support from the Russian Foundation for Basic Research (Grant No. 18-29-19198; Funder ID: 10.13039/501100002261) is deeply appreciated. TD acknowledges the support from NASA’s Space Technology Research (Grant No. NNX16AE03G; Funder ID: 10.13039/100000104) and US Scholar Fulbright program. The simulations presented in the paper were performed on computational clusters Pardus and Zhores (Skolkovo Institute of Science and Technology).

## References

1.
Wu
,
Z.
,
Chen
,
Z.
,
Du
,
X.
,
Logan
,
J. M.
,
Sippel
,
J.
,
Nikolou
,
M.
,
Kamaras
,
K.
,
Reynolds
,
J. R.
,
Tanner
,
D. B.
,
Hebard
,
A. F.
, and
Rinzler
,
A. G.
,
2004
, “
Transparent, Conductive Carbon Nanotube Films
,”
Science
,
305
(
5688
), pp.
1273
1276
. 10.1126/science.1101243
2.
Park
,
S.
,
Vosguerichian
,
M.
, and
Bao
,
Z.
,
2013
, “
A Review of Fabrication and Applications of Carbon Nanotube Film-Based Flexible Electronics
,”
Nanoscale
,
5
(
5
), pp.
1727
1752
. 10.1039/c3nr33560g
3.
Timmermans
,
M. Y.arina
,
Mariano
,
M.
,
Pollentier
,
I.
,
Richard
,
O.
,
Huyghebaert
,
C.
, and
Gallagher
,
E. E.
,
2018
, “
Free-standing carbon nanotube films for extreme ultraviolet pellicle application
,”
J. Micro-Nanolith MEM.
,
17
(
4
), pp.
1
8
. DOI:10.1117/1.JMM.17.4.043504
4.
Wang
,
L.
, and
Loh
,
K. J.
,
2017
, “
Wearable Carbon Nanotube-Based Fabric Sensors for Monitoring Human Physiological Performance
,”
Smart Mater. Struct.
,
26
(
5
), p.
055018
. 10.1088/1361-665X/aa6849
5.
Nasibulin
,
A. G.
,
,
A.
,
Mustonen
,
K.
,
Anisimov
,
A. S.
,
Ruiz
,
V.
,
Kivistö
,
S.
,
Rackauskas
,
S.
,
Timmermans
,
M. Y.
,
Pudas
,
M.
,
Aitchison
,
B.
,
Kauppinen
,
M.
,
Brown
,
D. P.
,
Okhotnikov
,
O. G.
, and
Kauppinen
,
E. I.
,
2011
, “
Multifunctional Free-Standing Single-Walled Carbon Nanotube Films
,”
ACS Nano
,
5
(
4
), pp.
3214
3221
. 10.1021/nn200338r
6.
Kis
,
A.
,
Jensen
,
K.
,
Aloni
,
S.
,
Mickelson
,
W.
, and
Zettl
,
A.
,
2006
, “
Interlayer Forces and Ultralow Sliding Friction in Multiwalled Carbon Nanotubes
,”
Phys. Rev. Lett.
,
97
(
2
), p.
025501
. 10.1103/PhysRevLett.97.025501
7.
Yakobson
,
B. I.
,
Brabec
,
C. J.
, and
Bernholc
,
J.
,
1996
, “
Nanomechanics of Carbon Tubes: Instabilities Beyond Linear Response
,”
Phys. Rev. Lett.
,
76
(
14
), pp.
2511
2514
. 10.1103/PhysRevLett.76.2511
8.
Dumitrica
,
T.
,
Hua
,
M.
, and
Yakobson
,
B. I.
,
2006
, “
Symmetry-, Time-, and Temperature-Dependent Strength of Carbon Nanotubes
,”
Proc. Natl. Acad. Sci. U. S. A.
,
103
(
16
), pp.
6105
6109
. 10.1073/pnas.0600945103
9.
Zhang
,
D.-B.
, and
Dumitrică
,
T.
,
2008
, “
Elasticity of Ideal Single-Walled Carbon Nanotubes Via Symmetry-Adapted Tight-Binding Objective Modeling
,”
Appl. Phys. Lett.
,
93
(
3
), p.
031919
. 10.1063/1.2965465
10.
Nikiforov
,
I.
,
Zhang
,
D.-B.
,
James
,
R. D.
, and
Dumitrică
,
T.
,
2010
, “
Wavelike Rippling in Multiwalled Carbon Nanotubes Under Pure Bending
,”
Appl. Phys. Lett.
,
96
(
12
), p.
123107
. 10.1063/1.3368703
11.
Cornwell
,
C. F.
, and
Welch
,
C. R.
,
2011
, “
Very-High-Strength (60-GPa) Carbon Nanotube Fiber Design Based on Molecular Dynamics Simulations
,”
J. Chem. Phys.
,
134
(
20
), p.
204708
. 10.1063/1.3594197
12.
Buehler
,
M. J.
,
2006
, “
Mesoscale Modeling of Mechanics of Carbon Nanotubes: Self-Assembly, Self-Folding, and Fracture
,”
J. Mater. Res.
,
21
(
11
), pp.
2855
2869
. 10.1557/jmr.2006.0347
13.
Cranford
,
S. W.
, and
Buehler
,
M. J.
,
2010
, “
In Silico Assembly and Nanomechanical Characterization of Carbon Nanotube Buckypaper
,”
Nanotechnology
,
21
(
26
), p.
265706
. 10.1088/0957-4484/21/26/265706
14.
Mirzaeifar
,
R.
,
Qin
,
Z.
, and
Buehler
,
M. J.
,
2015
, “
Mesoscale Mechanics of Twisting Carbon Nanotube Yarns
,”
Nanoscale
,
7
(
12
), pp.
5435
5445
. 10.1039/C4NR06669C
15.
Volkov
,
A. N.
, and
Zhigilei
,
L. V.
,
2010
, “
Structural Stability of Carbon Nanotube Films: The Role of Bending Buckling
,”
ACS Nano
,
4
(
10
), pp.
6187
6195
. 10.1021/nn1015902
16.
Volkov
,
A. N.
, and
Zhigilei
,
L. V.
,
2010
, “
Mesoscopic Interaction Potential for Carbon Nanotubes of Arbitrary Length and Orientation
,”
J. Phys. Chem. C
,
114
(
12
), pp.
5513
5531
. 10.1021/jp906142h
17.
Volkov
,
A. N.
, and
Zhigilei
,
L. V.
,
2010
, “
Scaling Laws and Mesoscopic Modeling of Thermal Conductivity in Carbon Nanotube Materials
,”
Phys. Rev. Lett.
,
104
(
21
), p.
215902
. 10.1103/PhysRevLett.104.215902
18.
Wittmaack
,
B. K.
,
Volkov
,
A. N.
, and
Zhigilei
,
L. V.
,
2018
, “
Mesoscopic Modeling of the Uniaxial Compression and Recovery of Vertically Aligned Carbon Nanotube Forests
,”
Compos. Sci. Technol.
,
166
(
1
), pp.
66
85
. 10.1016/j.compscitech.2018.03.014
19.
Ostanin
,
I.
,
Ballarini
,
R.
,
Potyondy
,
D.
, and
Dumitrica
,
T.
,
2013
, “
A Distinct Element Method for Large Scale Simulations of Carbon Nanotube Assemblies
,”
J. Mech. Phys. Solids
,
61
(
3
), pp.
762
782
. 10.1016/j.jmps.2012.10.016
20.
Ostanin
,
I.
,
Ballarini
,
R.
, and
Dumitrica
,
T.
,
2014
, “
Distinct Element Method Modeling of Carbon Nanotube Bundles With Intertube Sliding and Dissipation
,”
ASME J. Appl. Mech.
,
81
(
6
), p.
061004
. 10.1115/1.4026484
21.
Ostanin
,
I.
,
Ballarini
,
R.
, and
Dumitrica
,
T.
,
2015
, “
Distinct Element Method for Multiscale Modeling of Cross-Linked Carbon Nanotube Bundles: From Soft to Strong Nanomaterials
,”
J. Mater. Res.
,
30
(
1
), pp.
19
25
. 10.1557/jmr.2014.279
22.
Cundall
,
P. A.
, and
Strack
,
O. D.
,
1979
, “
A Discrete Numerical Model for Granular Assemblies
,”
Geotechnique
,
29
(
1
), pp.
47
65
. 10.1680/geot.1979.29.1.47
23.
Wang
,
Y.
,
Gaidau
,
C.
,
Ostanin
,
I.
, and
Dumitrica
,
T.
,
2013
, “
Ring Windings From Single-wall Carbon Nanotubes: A Distinct Element Method Study
,”
Appl. Phys. Lett.
,
103
(
18
), p.
183902
. 10.1063/1.4827337
24.
Wang
,
Y.
,
Semler
,
M. R.
,
Ostanin
,
I.
,
Hobbie
,
E. K.
, and
Dumitrica
,
T.
,
2014
, “
Rings and Rackets From Single-Wall Carbon Nanotubes: Manifestations of Mesoscale Mechanics
,”
Soft Matter
,
10
(
43
), pp.
8635
8640
. 10.1039/C4SM00865K
25.
Wang
,
Y.
,
Xu
,
H.
,
Drozdov
,
G.
, and
Dumitrică
,
T.
,
2018
, “
Mesoscopic Friction and Network Morphology Control the Mechanics and Processing of Carbon Nanotube Yarns
,”
Carbon
,
139
(
1
), pp.
94
104
. 10.1016/j.carbon.2018.06.043
26.
Ostanin
,
I.
,
Zhilyaev
,
P.
,
Petrov
,
V.
,
Dumitrica
,
T.
,
Eibl
,
S.
,
Ruede
,
U.
, and
Kuzkin
,
V.
,
2018
, “
Toward Large Scale Modeling of Carbon Nanotube Systems With the Mesoscopic Distinct Element Method
,”
Lett. Mater.
,
8
(
3
), pp.
240
245
. 10.22226/2410-3535
27.
Preclik
,
T.
, and
Rüde
,
U.
,
2015
, “
Ultrascale Simulations of Non-Smooth Granular Dynamics
,”
Comput. Part. Mech.
,
2
(
2
), pp.
173
196
. 10.1007/s40571-015-0047-6
28.
Itasca Consulting group Inc.
,
2008
,
PFC3D (Particle Flow Code in Three Dimensions). Version 4.0.
,
Itasca Consulting Group Inc.
,
Minneapolis, MN
.
29.
Brilliantov
,
N. V.
,
Budkov
,
Y. A.
, and
Seidel
,
C.
,
2017
, “
Theoretical and Numerical Analysis of Nano-Actuators Based on Grafted Polyelectrolytes in an Electric Field
,”
,
199
(
1
), pp.
487
510
. 10.1039/C6FD00240D
30.
Kuzkin
,
V. A.
, and
Asonov
,
I. E.
,
2012
, “
Vector-Based Model of Elastic Bonds for Simulation of Granular Solids
,”
Phys. Rev. E
,
86
(
5
), p.
051301
. 10.1103/PhysRevE.86.051301
31.
Potyondy
,
D. O.
, and
Cundall
,
P. A.
,
2004
, “
A Bonded-Particle Model for Rock
,”
Int. J. Rock Mech. Min. Sci.
,
41
(
8
), pp.
1329
1364
. 10.1016/j.ijrmms.2004.09.011
32.
Godenschwager
,
C.
,
Schornbaum
,
F.
,
Bauer
,
M.
,
Köstler
,
H.
, and
Rüde
,
U.
,
2013
, “
A Framework for Hybrid Parallel Flow Simulations with a Trillion Cells in Complex Geometries
,”
Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, ACM
,
Denver, CO
,
Nov. 17–22
,
ACM Press
, p.
35
.
33.
Schornbaum
,
F.
, and
Rüde
,
U.
,
2016
, “
Massively Parallel Algorithms for the Lattice Boltzmann Method on Non-Uniform Grids
,”
SIAM J. Sci. Comput.
,
38
(
2
), pp.
C96
C126
. 10.1137/15M1035240
34.
Eibl
,
S.
, and
Rüde
,
U.
,
2018
, “
A Local Parallel Communication Algorithm for Polydisperse Rigid Body Dynamics
,”
Parallel Comput.
,
80
(
1
), pp.
36
48
. 10.1016/j.parco.2018.10.002
35.
Ma
,
Y. J.
,
Yao
,
X. F.
,
Zheng
,
Q. S.
,
Yin
,
Y. J.
,
Jiang
,
D. J.
,
Xu
,
G. H.
,
Wei
,
F.
, and
Zhang
,
Q.
,
2010
, “
Carbon Nanotube Films Change Poisson’s Ratios From Negative to Positive
,”
Appl. Phys. Lett.
,
97
(
6
), p.
061909
. 10.1063/1.3479393
36.
Xie
,
B.
,
Liu
,
Y.
,
Ding
,
Y.
,
Zheng
,
Q.
, and
Xu
,
Z.
,
2011
, “
Mechanics of Carbon Nanotube Networks: Microstructural Evolution and Optimal Design
,”
Soft Matter
,
7
(
21
), pp.
10039
10047
. 10.1039/c1sm06034a
37.
Li
,
Y.
, and
Kröger
,
M.
,
2012
, “
Viscoelasticity of Carbon Nanotube Buckypaper: Zipping–Unzipping Mechanism and Entanglement Effects
,”
Soft Matter
,
8
(
30
), pp.
7822
7830
. 10.1039/c2sm25561h
38.
Sinclair
,
R. C.
,
Suter
,
J. L.
, and
Coveney
,
P. V.
,
2018
, “
Graphene–Graphene Interactions: Friction, Superlubricity, and Exfoliation
,”
,
30
(
13
), p.
1705791
39.
Chen
,
Y.
,
Pan
,
F.
,
Guo
,
Z.
,
Liu
,
B.
, and
Zhang
,
J.
,
2015
, “
Stiffness Threshold of Randomly Distributed Carbon Nanotube Networks
,”
J. Mech. Phys. Solids
,
84
(
1
), pp.
395
423
. 10.1016/j.jmps.2015.07.016
40.
Pan
,
F.
,
Chen
,
Y.
,
Liu
,
Y.
, and
Guo
,
Z.
,
2017
, “
Out-of-Plane Bending of Carbon Nanotube Films
,”
Int. J. Solids Struct.
,
106–107
(
1
), pp.
183
199
. 10.1016/j.ijsolstr.2016.11.020