High-Performance Numerical Modeling of Nanofabrics with Distinct Element Method

TL;DR

This paper introduces a high-performance, flexible numerical modeling framework for nanofabrics and composites, combining advanced fiber representation with efficient parallel computation, validated through a mechanical test on carbon nanotube textiles.

Contribution

The paper presents a novel, scalable modeling system for micro- and nanofiber composites using rigid body chains and enhanced vector formalism, integrated into the waLBerla framework.

Findings

High parallel performance demonstrated

Accurate modeling of nanofiber mechanics achieved

Effective simulation of carbon nanotube textiles

Abstract

A universal framework for modeling composites and fabrics of micro- and nanofibers, such as carbon nanotubes, carbon fibers and amyloid fibrils, is presented. Within this framework, fibers are represented with chains of rigid bodies, linked with elastic bonds. Elasticity of the bonds utilizes recently developed enhanced vector model formalism. The type of interactions between fibers is determined by their nature and physical length scale of the simulation. The dynamics of fibers is computed using the modification of rigid particle dynamics module of the waLBerla multiphysics framework. Our modeling system demonstrates exceptionally high parallel performance combined with the physical accuracy of the modeling. The efficiency of our approach is demonstrated with illustrative mechanical test on a hypothetical carbon nanotube textile.

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 the fiber interaction potential 5 .

$\epsilon$,$\left(eV\right)$	$A$	$B$	$\alpha$	$\beta$	$k,eV/{\text{Å}}^2$
$0.07124$	$0.0223$	$1.31$	$9.5$	$4.0$	$200$

Equations

$$r=\sqrt{2.5}{r}_f.$$

$$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}_f,$$

$$J=\pi h{r}_f(r_f^2+h^2/4),$$

$$J_p=2J.$$

$$U(r_{ij})=\begin{cases}\begin{array}[]{c}4\varepsilon\left(\frac{A}{(r_{ij}/r_{f}-2)^{\alpha}}-\frac{B}{(r_{ij}/r_{f}-2)^{\beta}}\right)\\ k(r_{ij}-r_{f})^{2}\end{array}&\begin{array}[]{c}r_{ij}>r_{f}\\ r_{ij}\leq r_{f}\end{array}\end{cases}$$

$$l_0=\sqrt{\frac{EJ}{\eta }}$$

Full text

1 Introduction

Fibrillar materials, based on biological fibrils, carbon fibers, nanofibers and carbon nanotubes (CNTs) [1, 2, 3, 4], and, in particular, textiles and fabrics made of individual fibrils or woven fibers, are of extreme interest for a number of military, aerospace, electronic and biomedical applications. However, intricate hierarchical structures of such materials, their discontinuous behavior with non-trivial inter-fiber interactions, as well as the prohibitively large sizes of representative volume elements in many cases prevent straightforward theoretical prediction of the mechanical, electrical and thermal properties of such materials. Understanding of the mesoscale behavior of these materials can be improved via numerical simulations. Atomic-level modeling techniques [5, 6, 7, 8], namely – tight-binding and molecular dynamics methods, were proved to be efficient numerical tools for modeling individual fibrils and their interactions at the nanoscale, however, the scalability of such techniques is insufficient for modeling large numbers of long fibers, necessary for studying the mechanics of the sufficiently large specimens of fibrillar materials. In order to address this problem, a number of mesoscale models were suggested. One of them, bead-spring model, employs the idea of coarse-grained molecular dynamics [9, 10, 11], initially proposed for modeling mesoscale mechanics of proteins. In this modeling concept, a chain of point masses represents a fibril interacting via classic potentials, representing either intra-fibril elasticity, or contact interactions between the neighboring fibrils. This model, despite its evident advantages, has certain limitations in a context of modeling fibrillar materials and fabrics. The most important of them is absent torsional stiffness of fibrils leading to unrealistic behaviors of fibrillar assemblies under certain loadings. In order to solve this issue, a different discretization concept [12, 13, 14, 15, 16, 17, 18, 19] was suggested using a representation of a thin fiber as a chain of rigid bodies, rather than point masses. Such a model allows not only bending of individual fibers, but their torsion as well. This simulation technique, known as mesoscopic distinct element method (MDEM), established itself in the field of modeling CNTs systems as one of the most efficient mesoscopic modeling tools, both computationally efficient and physically just. The technique can be successfully used for modeling a wide class of fibers and fibrillar materials, on the scales that admit athermal description of the fiber mechanics. Until now, the remaining obstacle on the route toward applications of MDEM to large-scale modeling of fibrillar assemblies was the absence of its scalable, parallel realization. Such realization was recently suggested in [19]. It is based on rigid particle dynamics module of the waLBerla multiphysics framework [20]. In the current work we illustrate the novel modeling approach in application to modeling nanofabrics – hypothetical textiles made of single wall carbon nanotubes (CNTs), ultimately strong nanofibers.

2 Method

Our model is based on the mesoscopic distinct element method, that computes the damped dynamics of a collection of interacting classical particles with certain mass and tensor of inertia. We utilize the spherical particles with the radius $r$, uniformly distributed mass $m$ and the scalar moment of inertia $I=\frac{2}{5}mr^{2}$. The state variables for each particle include translational positions and velocities, as well as rotational positions (in a shape of quaternions, as described in [20]) 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 explicit velocity Verlet time integration scheme. Two kinds of damping may be introduced in the system. The viscous damping forces, proportional to relative segment velocities, act in parallel with pair contact forces. In addition, PFC-style local damping forces [21] act at each body. In the simulations showcased below the viscous damping was absent and the constant of local damping was set to $0.4$.

Within our approach, undeformed fibers are partitioned into identical segments of finite lengths $T=2r_{f}$ and represented with chains of spherical rigid bodies (Fig. 1). Each spherical particle represents the inertial properties of a fiber segment - parameters $m$ and $I$ are equal with the mass and moment of inertia of a cylindrical segment taken with respect to the fiber axis. It follows that the spherical particle has a radius

[TABLE]

Elasticity of fibers in our model is represented with the formalism of enhanced vector model (EVM) [22, 23]. The EVM is based on a binding potential, describing the behavior of an elastic bond linking two rigid bodies. The formulation provides straightforward generalization on the case of large strains and accounts for a bending-twisting coupling. Consider two equal-sized spherical particles $i$ and $j$ with equilibrium separation $T$ and equilibrium orientation described in terms of orthogonal vectors $n_{ik}$, as depicted in Figure 1(B) (note that for an undeformed bond $n_{i1}=-n_{j1}$, $n_{i2}=n_{j2}$, $n_{i3}=n_{j3}$). Then the EVM bond potential is given as follows:

[TABLE]

Here $\mathbf{r}_{ij}$ is the radius vector connecting centers of bonded particles. The first term of the potential 2 accounts for the elastic strain energy stored due to axial tension/compression of a bond, second term is associated with shear of a bond, third term gives bond’s bending energy and the last term describes the energy associated with torsion of a bond. 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 papers [17, 22, 23] for more details):

[TABLE]

Here $E$ and $G$ are the bond material Young’s and shear moduli. Area $S$, moment of inertia $J$ and polar moment of inertia $J_{p}$ of a cylinder shell beam with radius $r_{f}$ and thickness $h$ are given by:

[TABLE]

As an example, consider here the parameterization of our discretization scheme for single-wall CNTs. Table 1 provides segment parameters for $(10,10)$ CNTs with diameter $2r_{f}=13.56$ $\mathring{A}$ and length $T=2r_{f}$. Each segment contains approximately $220$ carbon atoms. Microscopically computed Young’s $E=1,029$ GPa and shear modulus $G=459$ GPa [7] are used.

The interactions between fiber segments (e.g. Hertzian elastic repulsion, van der Waals (vdW) adhesion, wet surface tension, hydrogen bonds etc.) generally depend on the specifics of a particular problem. In our example, we utilize the combination of linear elastic repulsion between fiber segments at small distances, combined with vdW adhesion at large distances. The total potential of pair interaction is given by:

[TABLE]

For the distances that are less than two fiber radii the potential 5 describes linear elastic repulsion. The stiffness $k$ is fitted to ensure stable integration at a given timestep and the absence of fibers interpenetration. For the distances larger than two fiber radii the potential (5) describes the coarse-grained potential for vdW adhesion. The calibration of a coarse-grained isotropic vdW potential for $(10,10)$ CNTs is given in [13].

In order to capture geometric anisotropy of the cylindrical segments of fibers, we utilize simple numerical integration of the spherically-symmetric potential (5) over the length direction of segments. Three equispaced integration points along each segment’s axis are employed. Table 2 provides the parameterization of the potential.

Parallel damped dynamics simulations are based on the rigid particle dynamics module of waLBerla multiphysics framework, which is available under GPL license at (www.walberla.net). The parallelization is based on standard Message Passing Interface (MPI) [24] for distributed memory architectures. A complete description of the parallel algorithms and their realization [18] is beyond the scope of this paper. The simplified pipeline of our modeling framework is given in Fig. 2. The simulation starts with the generation of initial geometry of fibers, and imposition of boundary conditions. In the next step, the simulation domain is divided in a balanced manner into rectangular subdomains. These subdomains are distributed among the available MPI processes in such a way that every process is responsible for one or more subdomains. At the next stage, time integration cycles are performed on all MPI processes. The integration cycle consists of computation of pair interaction potentials, as well as the corresponding forces and moments at each contact. These forces and moments are then used to compute accelerations and angular accelerations that are then used in computing updated positions and velocities according to velocity Verlet time integration scheme. Particle migrations across subdomain borders are accounted via MPI communications. Then the list of contacts is updated. The contact detection scheme used in our simulations is based on hierarchical hash grids [25] and adapted for potential-based interactions. For correct contact detection of particles near the borders of a subdomain so-called ghost particles are introduced. These ghost particles mirror particles which touch the subdomain but are located at a different one. This way they are available for contact detection and force calculation. The configuration and traced quantities ( e.g. potential energy) are gathered at and saved periodically during the relaxation. The general scalability of waLBerla framework is proven up to almost half a million cores [20]. However, in our case some serial operations (e.g. gathering of the total potential and kinetic energy of the system) limit the efficiency of the parallelization.

3 Numerical results

As an example of application of our system, we consider here the relaxation and mechanical test on a hypothetical CNT textile material. Modern technologies do not yet allow to produce such textiles, however, our modeling framework allows to evaluate its properties in a mesoscopic simulation.

The simulations were performed at a computational cluster “Zhores” [26] using $20$ nodes. Every node uses two Intel Xeon 6136 Gold CPUs ($24$ cores, 3 GHz each). The high performance cluster network has the Fat Tree topology and is build from six Mellanox SB7890 (unmanaged) and two SB7800 (managed) switches that provide 100 Gbit/s connections between the nodes.

Consider a fragment of a CNT fabric, consisting of $400$ CNTs, $1.35$ $\mu m$ long each ($2.4\times 10^{6}$ model degrees of freedom), woven together into a square piece of a textile. Figure 2(A) shows an initial configuration, the magnified structure of a textile is shown on an inset. Such structure is stable when periodic boundary conditions are applied along $x$ and $y$ directions (infinite size approximation). However, it is interesting to answer the question about the stability of a finite-sized piece of such material.

The specimen shown in Figure 2(A) is allowed to relax in a damped dynamics simulation to a meta-stable state. At the initial stage of the simulation, edge CNTs start to separate from the fabric, since low vdW adhesion energy can not confine elastic strain energy, which is released during separation of side CNTs. Detached side CNTs form bundles comprising about $10-20$ tubes each. However, at the next stage of simulation the fabric disintegration process stops, while the relaxation slows down. CNTs do not separate from the fabric, since further separation is prevented by vdW adhesion. Similarly to the cases of other self-assembled CNT structures [14, 15, 16], CNT fabric achieves a meta-stable state, which is characterized by the balance between vdW adhesion energy and elastic strain energy. Such a balance is achieved for structure features of a certain size, characterized by the mesoscopic length scale

[TABLE]

where $EJ$ is the bending stiffness of a CNT, and $\eta$ is the vdW adhesion energy per unit length. This length scale arises explicitly in the analysis of elementary self-folded configurations - rings [13], rackets [16], multiple-winding rings [15]. For an individual $(10,10)$ CNTs, given the bending stiffness of $22350$ $eV\textup{\AA}$ and the adhesion energy of $0.22$ $eV/\textup{\AA}$, this length scale is equal to $0.032$ $\mu m$. For bundles, comprising multiple CNTs, both bending stiffness and characteristic length scale are somewhat larger.

Figure 2(C) gives the dependence of the elastic strain energy stored in separate CNTs during the relaxation. As we can see, the dependence is nearly exponential and strongly indicates the convergence to a stable state. Thus, we have a numerical evidence of the stability of hypothetical CNT fabrics.

In order to demonstrate in silico the exceptional mechanical properties of this material, we perform a numerical simulation of a large strain displacement controlled mechanical test on a specimen of CNT fabric material. Figures 4(A-C) illustrates the experiment setup. Self-assembled equilibrated CNT film specimen is subdivided into three regions - two grips, marked with green and red colors, and a gage region, marked with blue. Starting from the initial moment of the simulation, grips start to move in opposite directions, stretching the gage region. Grip velocity is kept constant, providing strain rate of $2\times 10^{8}$ $s^{-1}$, with exception for short constant acceleration period in the beginning of the simulation, necessary to avoid inertial peak at the beginning.

For this test, we introduced a simple breakage model for a CNT in assembly. An individual CNT breaks if it is stretched up to certain critical level. This critical strain has a normal distribution with the mean value $\epsilon_{c}$ and dispersion $\Delta\epsilon$; random distribution is introduced to qualitatively evaluate effects of finite temperature and CNT defects. In our test, $\epsilon_{c}=0.3$, $\Delta\epsilon=0.05$.

Figure 4(D) displays stress-strain curves during the simulation. Stress is defined in the assumption that the fabric’s thickness is equal to $17.1\textup{\AA}$ - two radii of a CNT plus equilibrium vdW separation. The initial step of the stress-strain curve is associated with dissipative response, conditioned by presence of local damping. The subsequent hardening region(strains of $2-20\%$) is associated with straightening of individual CNTs in fabric, oriented along the loading direction (Fig. 4(B)). At strains of $20-25\%$ we can see the elastic response conditioned by stretching of individual CNTs. Region of strains of $25-50\%$ features complex failure with correlated breakage of separate CNTs (Fig. 4(C)), first at the edges and then in the central part of the specimen. One can see that the hypothetical CNT fabric demonstrates the Yield strength of 15 GPa, which is three orders of magnitude higher that the strength of polymer films with similar flexibility and thickness [27]. Such exceptional properties are undoubtedly of a great interest for practical applications.

4 Conclusion

In this work we have presented the new general framework for modeling fibrous composites and textiles. Separate fibers are modeled as chains of interacting rigid bodies. The elasticity of individual fiber is represented with EVM formalism. The suggested framework is illustrated in the application to modeling of hypothetical CNT nanofabrics. In the benchmark example, the framework was capable to simulate relaxation and mechanical test on a CNT fabric specimen ( $2.4\times 10^{6}$ model degrees of freedom, approximately $10^{12}$ contact resolution computations) in approximately $20$ hours on $20$ nodes, with nearly linear scaling with the number of cores used. In the benchmark example considered, HPC capabilities of our framework allowed to discover stabilization of large specimens of hypothetical CNT fabric by vdW adhesion forces, and to perform the mechanical test indicating superior properties of hypothetical CNT textile. The proposed framework is capable to model efficiently any systems of interacting elastic fibers at any length and time scales admitting athermal description. Any types of contact interactions and nonlinearities in fiber’s constitutive behaviours can be straightforwardly incorporated into suggested modeling concept. Therefore, our framework can be straightforwardly applied to a wide class of problems, including composites, ropes and textiles for aerospace and military applications.

Acknowledgements

Authors acknowledge the financial support from the Russian Foundation for Basic Research (RFBR) under grants 16-31-60100 and 18-29-19198. Assistance of waLBerla package developers (S.Eibl and U.Rude) and Skoltech HPC team (R. Arslanov, I.Zacharov) is deeply appreciated. High-performance computations presented in the paper were carried out on Skoltech HPC cluster Zhores.