Size-Independent Mechanical Response of Ultrathin Carbon Nanotube Films in Mesoscopic Distinct Element Method Simulations

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 [7–10], 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 [12–18]. 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 [19–21], 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.

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 Å and length T = 2r_CNT. Each segment contains approximately 220 carbon atoms. The properties of such segments were parameterized in our previous works [19,20].

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.

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.

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)).

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.

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 × 10¹⁰ particles with up to 1.8 × 10⁶ 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.² 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.

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.

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 + ɛ)l − m; (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.

In the expressions above, e_x and e_y are unit vectors directed along x and y coordinate axes. In the simulations described below, the time t₀ 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.

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 10⁶ 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.

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/cm³) was prescribed. Unless otherwise noted, the specimens have a thickness of 20 nm and the surface density of 0.088 μg/cm².

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)).

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).

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)).

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.

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 (l_CNT > 1 μm) and RVE size larger than the length of a CNT.

Figure 10(d) indicates the values of Y^uniax = 17.5 GPa and Y^biax = 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].

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 [15–18] 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 [15–18] 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.

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.

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).