Computational Study of Ultrathin CNT Films with the Scalable Mesoscopic Distinct Element Method

TL;DR

This study uses a scalable mesoscopic simulation method to analyze the small strain mechanics of ultrathin carbon nanotube films, revealing how energy dissipation during CNT sliding influences their mechanical behavior.

Contribution

It introduces a detailed mesoscopic modeling approach for large CNT assemblies, incorporating realistic inter-tube interactions and deformation mechanisms.

Findings

Energy dissipation during CNT sliding extends load transfer.

CNT films exhibit material-like mechanical response due to weak bonding.

The simulation provides insights into the micro-mechanical behavior of CNT assemblies.

Abstract

In this work we present a computational study of the small strain mechanics of freestanding ultrathin 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 CNT surface and utilize Enhanced Vector Model (EVM) formalism. Segments of the neighboring CNTs interact with realistic coarse-grained anisotropic vdW potential, enabling relative slip of CNTs in contact. The advanced simulation technique allowed us to gain useful insights on the behavior of CNT materials. In particular, it was established that the energy dissipation during CNT sliding leads to extended load transfer that conditions material-like mechanical response of the weakly bonded assemblies of CNTs.

Figures

Tables

Table 1. Table 1: Parameterization of the spherical particles and EVM bonds for a ( 10 , 10 ) 10 10 (10,10) CNTs. m 𝑚 m , r 𝑟 r , I 𝐼 I are the mass, radius, moment of inertia of each spherical particle. B 1 subscript 𝐵 1 B_{1} , B 2 subscript 𝐵 2 B_{2} , B 3 subscript 𝐵 3 B_{3} , B 4 subscript 𝐵 4 B_{4} are EVM stiffnesses.

$m$	$r$	$I$	$B_1$	$B_2$	$B_3$	$B_4$
$\left(amu\right)$	$\left(\text{Å}\right)$	$\left(amu\times {\text{Å}}^2\right)$	$\left(eV/{\text{Å}}^2\right)$	$\left(eV\right)$	$\left(eV\right)$	$\left(eV\right)$
$2,649$	$10.72$	$1.218\times {10}^5$	$67.59$	$19780$	$-4032$	$1471$

Table 2. Table 2: Parameterization of vdW potential for ( 10 , 10 ) 10 10 (10,10) CNTs.

${ϵ}^{\prime }$, meV	${\alpha }^{\prime }$	${\beta }^{\prime }$	$A^{\prime }$	$B^{\prime }$	$\delta$	$C_{\varphi }$
$149.3$	$9.5$	$4$	$0.0223$	$1.31$	$-7.5$	$90$
$k$	$C_1$	$C_2$	$C_3$	$C_4$	$C_5$
$5$	$0.35818$	$0.03263$	$-0.00138$	$-0.00017$	$0.00024$
$Q_1$	$Q_2$	$Q_3$	$Q_4$
$-80.0$	$288.0$	$-336.0$	$128$

Equations

$$r=\sqrt{2.5}{r}_{CNT}.$$

$$F_i^{\alpha }=-\alpha \left|F_i\right|sign(v_i),M_i^{\alpha }=-\alpha \left|M_i\right|sign({\omega }_i)$$

$$F_n^{\beta }=c_n{v}_n,F_s^{\beta }=c_s{v}_s,$$

$$c_n=2\beta \sqrt{m{k}_n},c_s=2\beta \sqrt{m{k}_s},$$

$$U(r_{ij},n_{ik},n_{jk})=\frac{B_1}{2}(r_{ij}-T{)}^2+\frac{B_2}{2}({\mathbf{n}}_{j1}-{\mathbf{n}}_{i1}){\mathbf{r}}_{ij}/r_{ij}+B_3{\mathbf{n}}_{i1}{\mathbf{n}}_{j1}-\frac{B_4}{2}({\mathbf{n}}_{i2}{\mathbf{n}}_{j2}-{\mathbf{n}}_{i3}{\mathbf{n}}_{j3})$$

$$B_1=\frac{ES}{T},$$

$$B_2=\frac{12EJ}{T},$$

$$B_3=-\frac{2EJ}{T}-\frac{G{J}_p}{2T},$$

$$B_4=\frac{G{J}_p}{T}.$$

$$S=2\pi h{r}_{CNT},$$

$$J=\pi h{r}_{CNT}(r_{CNT}^2+h^2/4),$$

$$J_p=2J.$$

$$U(r_{ij},\theta ,\varphi )$$

$$V^k(r_{ij},\theta )$$

$$D^k(r_{ij},\theta )$$

$$\Phi (r_{ij},\varphi )$$

$${\Theta }^k(\theta )$$

$$W_{\varphi }(r_{ij})$$

$$f_c(r_{ij})$$

$$v(t)=\{\begin{array}{ll}\frac{2x}{l}{v}_0& t>t_0\\ \frac{2x}{l}\frac{v_0}{t_0}t& t\le t_0\end{array}$$

$$v(t,x,y)=\{\begin{array}{ll}\frac{2x}{l}\frac{2y}{l}{v}_0& t>t_0\\ \frac{2x}{l}\frac{2y}{l}\frac{v_0}{t_0}t& t\le t_0\end{array}$$

$$Y^{uniax}=\frac{E}{1-{\nu }^2},Y^{biax}=\frac{E}{1-\nu }$$

Full text

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], mirrors for soft X-ray radiation [3], wearable medical devices [4] etc. The major advantage of these materials, as compared to possible competitors, is their ability to withstand significant plastic deformations without irreversible structural changes. This property is conditioned by discontinuous nature of CNT films - they consist of individual CNTs that are chemically reactive and bound only by van der Waals forces. Due to 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 mechanism for plastic flow of CNT material. This mechanism is non-destructive in a sense that small strain plastic deformations do not lead to significant changes in microstructure of a CNT material, as well as in its mechanical, thermal, optical and electric properties. Another interesting feature is that such plastic deformations do not immediately localize in a small damage zone, but for a considerable strain diapason remain distributed over the length of a tested specimen.

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 [5, 6, 7, 8] they’re inefficient for modeling assemblies comprising large numbers of CNTs - the most extreme molecular dynamics calculations so far included only hundreds of CNTs [9].

In order to bypass the computational inefficiency of full molecular modeling techniques, a number of coarse-grained mesoscale models were developed [10, 11, 12, 13, 14, 15, 16]. The general idea of these models is to reduce the number of model degrees of freedom (DOF) in order to reach the computationally tractable model size. However, such coarse-graining often leads to significant artifacts in model’s behavior. For example, pair vdW potential used in [10, 11] leads to absence of 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 [13, 14] 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 [17, 18, 19] we have developed the mesoscopic modeling technique based on distinct element method [20]. 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 adequate description of relative CNT sliding. In our athermal model, energy dissipation was represented with viscous damping, appropriately calibrated in a top-down fashion [18]. The model showed to be efficient for the description of CNT self-folded configurations [21, 22], weakly bonded and cross-linked CNT ropes and bundles [18, 19] and qualitative description of large-strain deformations of CNT films [23]. 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 [24].

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. We restrict our consideration with the case of small strain deformation of CNT materials, with strains less than few percent. We could see that even at such strains CNT films exhibit complex nonlinear behavior, conditioned by CNT sliding within the tested specimen. Our observations are in good agreement with recent experimental studies on the mechanics of CNT films.

2 Method

Our mesoscopic model mostly corresponds to the ones employed in our earlier works [17, 18, 24]. For clarity, we present a brief description of our modeling technique here. The model is based on the Mesoscopic Distinct Element Method (MDEM), 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.

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 $2r_{CNT}=13.560$ $\mathring{A}$ and length $T=2r_{CNT}$. Each segment contains approximately $220$ carbon atoms.

The inertial properties of a segment are represented with mass $m$ and the moment of inertia $I=\frac{2}{5}mr^{2}$ (a spherical tensor of inertia of a distinct element is assumed here for simplicity, however, a number of non-spherical shapes are immediately available in the framework). 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 proper choice of the radius of a distinct element:

[TABLE]

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 $\alpha$ and $\beta$ respectively.

PFC-style local damping [25] acts at each body. It is introduced to damp stiff interactions and stabilize time integration in dense particle assemblies. The components of damping force $F^{\alpha}_{i}$ (moment $M^{\alpha}_{i}$ ) are proportional to the corresponding components of unbalanced force $F_{i}$ (moment $M_{i}$) according to:

[TABLE]

Here $v_{i}$ and $\omega_{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 ($\alpha=0.1$) and its influence on CNT dynamics is insignificant.

Recent theoretical findings [26] 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 $\beta$. Normal $F^{\beta}_{n}$ and tangential $F^{\beta}_{s}$ viscous forces are proportional to normal $v_{n}$ and tangential $v_{s}$ relative velocities of elements in vdW contact:

[TABLE]

Viscosity coefficients $c_{n}$ and $c_{s}$ are related to $\beta$ as follows:

[TABLE]

where $k_{n}$ and $k_{s}$ are linearized stiffnesses of the contact model, taken $k_{n}=100eV/nm^{2}$, $k_{s}=100eV/nm^{2}$.

Damping factor $\beta$ is kept as a free parameter in our simulations. Initial estimate of $\beta$ is performed in a top-down way, as discussed in [18]. Unless otherwise noted $\beta=0.1$ is assumed in this work.

Elasticity of CNTs in our model is represented with the formalism of enhanced vector model (EVM) [27]. Unlike incremental parallel bonds, widely used in DEM community [28], 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 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 $n_{ik}$, as depicted in fig. 2(A). Then the EVM bond potential is given as follows:

[TABLE]

Here $\mathbf{r}_{ij}$ is the radius vector connecting centers of bonded particles. Parameters $B_{1}...B_{4}$ are directly related to longitudinal, shear, bending, and torsional rigidities of a bond, according to Euler-Bernoulli beam theory (see [27] for more details):

[TABLE]

where $S,I,J$ are area, moment of inertia and polar moment of inertia of a cylinder shell beam with radius $r_{CNT}$ and thickness $h$:

[TABLE]

Table 1 summarizes the elastic properties of CNT segments.

The potential 5 with calibrations 2 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 [17, 18]. The potential is defined in terms of three variables $R$, $\theta$ and $\phi$, related to mutual position and orientation of two cylindrical segments (fig. 2(C)):

[TABLE]

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 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 incommensurate mode of CNT slip. Such model does not create any static friction between CNTs and therefore will underestimate the strength of CNT materials in comparison with 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.

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 [29]. It offers a rigid particle dynamics module which is capable of simulating up to $2.8\text{\times}{10}^{10}$ particles with up to $1.8\text{\times}{10}^{6}$ parallel processes [30]. A detailed description of the algorithms and their implementation used in the framework can be found in [30, 31, 32]. The framework is freely available under GPLv3 license at www.walberla.net. 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 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 $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 advanced algorithms based on hierarchical hash grids and bounding volumes are used. The contact detection scheme used for 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 potential interaction cutoff radius. After all contacts are determined the contacts are used to calculate the forces and moments acting at each contact. For the calculation of the forces and moments the interaction model described above is used. In the next step the forces and moments are used to calculate accelerations and angular accelerations and the particle position and rotation gets updated using 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 all information but also migrates particles between processes if the subdomain they belong to has changed. Throughout the whole simulation information about the current particle configuration is gathered on a regular interval and saved to log files for later analysis.

3 Numerical results

3.1 Simulations setup

Below 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. Fig. 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\times l\times 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.

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+\varepsilon)l-m;(1+\varepsilon)l+m)$, where $\varepsilon$ is tensile strain and $m$ is 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)), that 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:

[TABLE]

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:

[TABLE]

In the simulations described below the time $t_{0}$ is set to be 500 timesteps (10 ps).

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

3.2 Self assembly of a CNT film specimen

Fig. 5 and fig. 6 summarize the initial relaxation of a CNT film specimen. Fig. 5 demonstrates the part of $2\times 2$ $\mu m$ specimen at the beginning of the relaxation process (A) and after $10^{6}$ cycles ($20$ ns) of relaxation (B,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. 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 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\times 0.5$ $\mu m$ to $3\times 3$ $\mu m$ and CNT lengths of $100$, $200$, $400$ and $800$ CNT segments. The same volumetric density ($0.044g/cm^{3}$ was prescribed. Unless otherwise noted, the specimens have thickness of $20nm$ and surface density of $0.088\mu g/cm^{2}$.

3.3 Mechanical tests on CNT 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 relative slip of the CNTs in the specimen. The resistance to slip is conditioned by 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 strain rate [17, 18]. 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 viscous part of the response, the observed stress-strain curves start from nonzero stress.

Fig. 7 summarizes the mechanics of CNT films observed in uniaxial tests. The Young’s modulus of CNT films observed in our simulations varies in a large diapason of $5-30$ GPa, and conditioned by many factors, including strain rate (fig. 7(A)), CNT length (fig. 7(B)) and specimen size (fig. 7(C)). The dependence of the Young’s modulus on thickness is rather negligible due to out-of-plane periodic boundary conditions (fig. 7(D)).

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 random microstructure of the film. Relative slip of CNTs can be observed at large magnifications (fig. 8).

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. 10(A-C) demonstrates 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\times 1.5\mu m$ and $3.0\times 3.0\mu m$. At snapshots (A-C) the simulation is carried out at $100\%$,$50\%$ and $25\%$ strain rate, respectively. The viscous dissipation coefficient $\beta$ 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 negligible role of inertia-dependent CNT rearrangement processes. It is therefore clear that strain localization and size-dependence of stress-strain curves is associated with inertial effects, that are negligible in slow and highly damped simulations.

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 validity of linear isotropic elastic relations for these initial slopes:

[TABLE]

Fig. 10(D) indicates the values of $Y^{uniax}=17.5GPa$ and $Y^{biax}=18.75GPa$ Therefore, it follows that the film 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 has been observed experimentally [33].

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 M. Buehler and co-workers [10, 11]. 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 [34, 35]. However, due to the artifact of vdW potential corrugation, that is described in detail in [17], 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 [35]. The other group of authors [13, 14, 15, 16] 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. This model assumptions were largely based on two observations: i) The CNT sliding observed in MD simulations does not feature significant viscous friction, and ii) given the absent 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 recent work [36], that indicates that 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 then 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 [13, 14, 15, 16] 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 [18] 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.

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 [26] on nanofriction confirm linear dependence between friction force and relative sliding velocity in ”thermal” regime of friction between nanosurfaces. However, precise calibration of such dissipation is problematic, since the observed viscosity coefficient depends on multiple factors and require accurate statistic 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 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 on its size. Similarly to earlier studied CNT rope structures [18], we observe an important role of energy dissipation, that is responsible for even load distribution along the specimen and its uniform deformation, leading to 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, assembled and ordered only by short-range vdW interactions. In is clear, that realistic CNT assemblies have a lot more complicated morphology, conditioned not only by vdW adhesion on the deposition stage, but also 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 realistic CNT film under mechanical loading. However, the most important deficiency of our model in its current shape is the absence of accurate model for the energy dissipation and its bottom-up calibration, that would allow us to reproduce the mechanics for CNT materials from the first principles. The development of such model, taking into account static friction for different CNT chiralities and commensurate/incommensurate contacts is the matter of future work.

Acknowledgements

I.O. acknowledges Russian Science Foundation (grant 17-73-10442) for the support of the development of the parallel system of modeling carbon nanotubes based on waLBerla framework. The financial support from the Russian Foundation for Basic Research under grant 18-29-19198 is deeply appreciated. The simulations presented in the paper were performed on computational clusters Pardus and Zhores (Skolkovo Institute of Science and Technology).