High-Performance Numerical Modeling of Nanofabrics with Distinct Element Method
HHigh-Performance Numerical Modeling of Nanofabrics withDistinct Element Method
Igor Ostanin Skolkovo Institute of Science and Technology, Nobel St. 3, Moscow, Russia
Abstract
A universal framework for modeling composites and fabrics of micro- and nanofibers, such as carbonnanotubes, carbon fibers and amyloid fibrils, is presented. Within this framework, fibers are representedwith chains of rigid bodies, linked with elastic bonds. Elasticity of the bonds utilizes recently developedenhanced vector model formalism. The type of interactions between fibers is determined by their natureand physical length scale of the simulation. The dynamics of fibers is computed using the modification ofrigid particle dynamics module of the waLBerla multiphysics framework. Our modeling system demonstratesexceptionally high parallel performance combined with the physical accuracy of the modeling. The efficiencyof our approach is demonstrated with illustrative mechanical test on a hypothetical carbon nanotube textile.
Keywords:
Nanofibers, Carbon Nanotubes, Distinct Element Method
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 extremeinterest for a number of military, aerospace, electronic and biomedical applications. However, intricate hier-archical structures of such materials, their discontinuous behavior with non-trivial inter-fiber interactions, aswell as the prohibitively large sizes of representative volume elements in many cases prevent straightforwardtheoretical prediction of the mechanical, electrical and thermal properties of such materials. Understand-ing of the mesoscale behavior of these materials can be improved via numerical simulations. Atomic-levelmodeling techniques [5, 6, 7, 8], namely – tight-binding and molecular dynamics methods, were proved tobe 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 forstudying the mechanics of the sufficiently large specimens of fibrillar materials. In order to address thisproblem, a number of mesoscale models were suggested. One of them, bead-spring model, employs the ideaof coarse-grained molecular dynamics [9, 10, 11], initially proposed for modeling mesoscale mechanics ofproteins. In this modeling concept, a chain of point masses represents a fibril interacting via classic poten-tials, representing either intra-fibril elasticity, or contact interactions between the neighboring fibrils. Thismodel, despite its evident advantages, has certain limitations in a context of modeling fibrillar materials andfabrics. The most important of them is absent torsional stiffness of fibrils leading to unrealistic behaviorsof 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 aswell. This simulation technique, known as mesoscopic distinct element method (MDEM), established itselfin the field of modeling CNTs systems as one of the most efficient mesoscopic modeling tools, both com-putationally efficient and physically just. The technique can be successfully used for modeling a wide class Corresponding author, e-mail:[email protected]
Preprint submitted to Elsevier a r X i v : . [ phy s i c s . c o m p - ph ] F e b igure 1: a) Schematics of a fiber discretization b) representation of an elastic bond between two bodies of fibers and fibrillar materials, on the scales that admit athermal description of the fiber mechanics. Untilnow, the remaining obstacle on the route toward applications of MDEM to large-scale modeling of fibrillarassemblies 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 cur-rent work we illustrate the novel modeling approach in application to modeling nanofabrics – hypotheticaltextiles 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 ofa collection of interacting classical particles with certain mass and tensor of inertia. We utilize the sphericalparticles with the radius r , uniformly distributed mass m and the scalar moment of inertia I = mr . Thestate 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 velocitiesand angular velocities due to contact forces and moments arising in pair interactions, as well as externalforces and moments, acting at each body. The system is evolved in time with explicit velocity Verlet timeintegration 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 wasabsent and the constant of local damping was set to 0 . T = 2 r f and represented with chains of spherical rigid bodies (Fig. 1). Each spherical particle represents the inertialproperties of a fiber segment - parameters m and I are equal with the mass and moment of inertia of acylindrical segment taken with respect to the fiber axis. It follows that the spherical particle has a radius r = √ . r f . (1)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 rigidbodies. The formulation provides straightforward generalization on the case of large strains and accounts fora 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) (note2 able 1: Parameterization of the spherical particles and EVM bonds for a (10 ,
10) CNTs. m , r , I are the mass, radius, momentof inertia of each spherical particle. B , B , B , B are EVM stiffnesses. m r I B B B B ( amu ) (˚A) ( amu × ˚A ) ( eV / ˚A ) ( eV ) ( eV ) ( eV )2 ,
649 10 .
72 1 . × .
59 19780 − Table 2: Parameterization of the fiber interaction potential 5. ε ,( eV ) A B α β k, eV / ˚A . . .
31 9 . . n i = − n j , n i = n j , n i = n j ). Then the EVM bond potential is given asfollows: U ( r ij , n ik , n jk ) = B r ij − T ) + B n j − n i ) r ij /r ij + B n i n j − B n i n j − n i n j ) (2)Here r ij is the radius vector connecting centers of bonded particles. The first term of the potential 2accounts for the elastic strain energy stored due to axial tension/compression of a bond, second term isassociated with shear of a bond, third term gives bond’s bending energy and the last term describes theenergy associated with torsion of a bond. Parameters B ...B 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): B = EST ,B = 12 EJT ,B = − EJT − GJ p T ,B = GJ p T . (3)Here E and G are the bond material Young’s and shear moduli. Area S , moment of inertia J and polarmoment of inertia J p of a cylinder shell beam with radius r f and thickness h are given by: S = 2 πhr f ,J = πhr f ( r f + h / ,J p = 2 J. (4)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 2 r f = 13 .
56 ˚ A and length T = 2 r f .Each segment contains approximately 220 carbon atoms. Microscopically computed Young’s E = 1 ,
029 GPaand 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 ourexample, 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:3 igure 2: Simplified diagram of the modeling framework pipeline. U ( r ij ) = (cid:40) ε (cid:16) A ( r ij /r f − α − B ( r ij /r f − β (cid:17) k ( r ij − r f ) r ij > r f r ij ≤ r f (5)For the distances that are less than two fiber radii the potential 5 describes linear elastic repulsion. Thestiffness 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 vdWadhesion. 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 numericalintegration of the spherically-symmetric potential (5) over the length direction of segments. Three equispacedintegration points along each segment’s axis are employed. Table 2 provides the parameterization of thepotential.Parallel damped dynamics simulations are based on the rigid particle dynamics module of waLBerlamultiphysics framework, which is available under GPL license at ( ). The parallelizationis based on standard Message Passing Interface (MPI) [24] for distributed memory architectures. A completedescription of the parallel algorithms and their realization [18] is beyond the scope of this paper. Thesimplified pipeline of our modeling framework is given in Fig. 2. The simulation starts with the generationof initial geometry of fibers, and imposition of boundary conditions. In the next step, the simulation domainis divided in a balanced manner into rectangular subdomains. These subdomains are distributed among theavailable MPI processes in such a way that every process is responsible for one or more subdomains. At thenext stage, time integration cycles are performed on all MPI processes. The integration cycle consists ofcomputation 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 thenused 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 ofcontacts is updated. The contact detection scheme used in our simulations is based on hierarchical hashgrids [25] and adapted for potential-based interactions. For correct contact detection of particles near theborders of a subdomain so-called ghost particles are introduced. These ghost particles mirror particles whichtouch the subdomain but are located at a different one. This way they are available for contact detectionand force calculation. The configuration and traced quantities ( e.g. potential energy) are gathered at andsaved periodically during the relaxation. The general scalability of waLBerla framework is proven up toalmost half a million cores [20]. However, in our case some serial operations ( e.g. gathering of the totalpotential and kinetic energy of the system) limit the efficiency of the parallelization.4 . Numerical results
As an example of application of our system, we consider here the relaxation and mechanical test on ahypothetical 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 nodeuses two Intel Xeon 6136 Gold CPUs (24 cores, 3 GHz each). The high performance cluster network has theFat Tree topology and is build from six Mellanox SB7890 (unmanaged) and two SB7800 (managed) switchesthat provide 100 Gbit/s connections between the nodes.Consider a fragment of a CNT fabric, consisting of 400 CNTs, 1 . µm long each (2 . × model degreesof freedom), woven together into a square piece of a textile. Figure 2(A) shows an initial configuration, themagnified structure of a textile is shown on an inset. Such structure is stable when periodic boundaryconditions are applied along x and y directions (infinite size approximation). However, it is interesting toanswer 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-stablestate. At the initial stage of the simulation, edge CNTs start to separate from the fabric, since low vdWadhesion 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 ofsimulation the fabric disintegration process stops, while the relaxation slows down. CNTs do not separatefrom the fabric, since further separation is prevented by vdW adhesion. Similarly to the cases of otherself-assembled CNT structures [14, 15, 16], CNT fabric achieves a meta-stable state, which is characterizedby the balance between vdW adhesion energy and elastic strain energy. Such a balance is achieved forstructure features of a certain size, characterized by the mesoscopic length scale l = (cid:115) EJη (6)where EJ is the bending stiffness of a CNT, and η is the vdW adhesion energy per unit length. Thislength 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 ˚A andthe adhesion energy of 0 . eV / ˚A, this length scale is equal to 0 . µm . For bundles, comprising multipleCNTs, 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 re-laxation. As we can see, the dependence is nearly exponential and strongly indicates the convergence to astable 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 anumerical simulation of a large strain displacement controlled mechanical test on a specimen of CNT fabricmaterial. Figures 4(A-C) illustrates the experiment setup. Self-assembled equilibrated CNT film specimenis subdivided into three regions - two grips, marked with green and red colors, and a gage region, markedwith 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 × s − , with exceptionfor short constant acceleration period in the beginning of the simulation, necessary to avoid inertial peak atthe beginning.For this test, we introduced a simple breakage model for a CNT in assembly. An individual CNT breaksif it is stretched up to certain critical level. This critical strain has a normal distribution with the mean value (cid:15) c and dispersion ∆ (cid:15) ; random distribution is introduced to qualitatively evaluate effects of finite temperatureand CNT defects. In our test, (cid:15) c = 0 .
3, ∆ (cid:15) = 0 . . − −
25% we can see the elastic response5 igure 3: Spontaneous disintegration of CNT fabric with subsequent stabilization. (A) initial configuration (on the inset -magnified structure of the fabric). (B) Final relaxed configuration of a finite-sized piece of CNT fabric. (C) Elastic strainenergy stored in CNTs as a function of time (D) Magnified piece an equilibrated configuration of a CNT fabric, featuringbundling of CNTs within the fabric structure. igure 4: Mechanical test on a CNT fabric specimen (A,B,C) - structure snapshots for engineering strains of 0%, 20%, and40% respectively. (D) Stress-strain curve monitored during the test −
50% features complex failure withcorrelated breakage of separate CNTs (Fig. 4(C)), first at the edges and then in the central part of thespecimen. One can see that the hypothetical CNT fabric demonstrates the Yield strength of 15 GPa, whichis 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 rep-resented with EVM formalism. The suggested framework is illustrated in the application to modeling ofhypothetical CNT nanofabrics. In the benchmark example, the framework was capable to simulate relax-ation and mechanical test on a CNT fabric specimen ( 2 . × model degrees of freedom, approximately10 contact resolution computations) in approximately 20 hours on 20 nodes, with nearly linear scalingwith the number of cores used. In the benchmark example considered, HPC capabilities of our frameworkallowed to discover stabilization of large specimens of hypothetical CNT fabric by vdW adhesion forces, andto perform the mechanical test indicating superior properties of hypothetical CNT textile. The proposedframework is capable to model efficiently any systems of interacting elastic fibers at any length and timescales admitting athermal description. Any types of contact interactions and nonlinearities in fiber’s con-stitutive behaviours can be straightforwardly incorporated into suggested modeling concept. Therefore, ourframework can be straightforwardly applied to a wide class of problems, including composites, ropes andtextiles 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 computationspresented in the paper were carried out on Skoltech HPC cluster Zhores.
ReferencesReferences [1] S. Iijima. Helical microtubules of graphitic carbon. Nature 354, 56–58, 1991.[2] R. Baughman. Carbon nanotubes – the route towards application. Science 297, 787–792, 2002.[3] S. Eppell, B. Smith, H. Kahn, R. Ballarini. Nano measurements with micro-devices: mechanical properties of hydratedcollagen fibrils. Journ. Roy. Soc. Int. 3(6), 117–121, 2006.[4] S. Chand. Review carbon fibers for composites. Journ. Mater. Sci. 35(6), 1303-1313, 2000.[5] B. Yakobson, C. Brabec, J. Bernholc. Nanomechanics of Carbon Tubes: Instabilities beyond Linear Response. Phys. Rev.Lett. 76(14), 2511, 1996.[6] T. Dumitrica, M. Hua, B. Yakobson. Symmetry-, time-, and temperature-dependent strength of carbon nanotubes. Proc.Natl. Acad. Sci. U.S.A. 103(16), 6105, 2006.[7] D. Zhang, T. Dumitrica. Elasticity of ideal single-walled carbon nanotubes via symmetry-adapted tight-binding objectivemodeling. Appl. Phys. Lett., 93, 031919, 2008.[8] I. Nikiforov, D. Zhang, R. James, T. Dumitrica. Wavelike rippling in multiwalled carbon nanotubes under pure bending.Appl. Phys. Lett. 96, 123107, 2010.[9] M. Buehler. Mesoscale modeling of mechanics of carbon nanotubes: Self-assembly, self-folding and fracture. Journ. Mat.Res. 21(11), 2855, 2006.[10] S. Cranford, M. Buehler. In silico assembly and nanomechanical characterization of carbon nanotube buckypaper. Nan-otechnology 21, 265706, 2010.[11] R. Mirzaeifar, Z. Qin, M. Buehler. Mesoscale mechanics of twisting carbon nanotube yarns. Nanoscale 7(12), 5435, 2015.[12] T. Anderson, E. Akatyeva, I. Nikiforov, D. Potyondy, R. Ballarini, T. Dumitrica. Toward distinct element methodsimulations of carbon nanotube systems. Journ. Nanotech. Eng. Med. 1(4), 0410009, 2010.[13] I. Ostanin, R. Ballarini, D. Potyondy, T. Dumitrica. A distinct element method for large scale simulations of carbonnanotube assemblies. Mech. Phys. Sol. 61(3), 762, 2013.
14] I. Ostanin, R. Ballarini, T. Dumitrica. Distinct element method modeling of carbon nanotube bundles with intertubesliding and dissipation. Appl. Mech. 81(6), 061004, 2014.[15] Y. Wang, C. Gaidau, I. Ostanin, T. Dumitrica. Ring windings from single-wall carbon nanotubes: A distinct elementmethod study. Appl. Phys. Lett. 103 (18), 183902, 2013.[16] Y. Wang, M. Semler, I. Ostanin, E. Hobbie, T. Dumitrica. Rings and rackets from single-wall carbon nanotubes: mani-festations of mesoscale mechanics. Soft Matter 10 (43), 8635-8640, 2014.[17] I. Ostanin, R. Ballarini, T. Dumitrica. Rings and rackets from single-wall carbon nanotubes: manifestations of mesoscalemechanics. Journ. Mat. Res. 30(1), 19, 2015.[18] Y. Wang, I. Ostanin, C. Gaidau, T. Dumitrica. Twisting carbon nanotube ropes with the mesoscopic distinct elementmethod: Geometry, packing, and nanomechanics. Langmuir, 31(45), 12323, 2015.[19] I. Ostanin, P. Zhilyaev, V. Petrov, T. Dumitrica, S. Eibl, U. Ruede, V. Kuzkin. Toward large scale modeling of carbonnanotube systems with the mesoscopic distinct element method. Mater. 8(3), 240-245, 2018.[20] T. Preclik, U. Ruede. Ultrascale simulations of non-smooth granular dynamics. Comp. Part. Mech., 2, 173, 2015.[21] Itasca Consulting Group Inc., 2015. PFC3D (Particle Flow Code in Three Dimensions). Version 5.0. Itasca ConsultingGroup Inc., Minneapolis.[22] V. Kuzkin, I. Asonov. Vector-based model of elastic bonds for simulation of granular solids. Phys. Rev. E, 86(5), 051301,2012.[23] V. Kuzkin, A. Krivtsov. Enhanced vector-based model for elastic bonds in solids. Lett. Mat. 7(4), 455, 2017.[24] MPI Forum. MPI: A message-passing interface standard. Technical report, Knoxville, TN, USA, 1994.[25] C. Ericson, Real-time collision detection. CRC Press, 2004.[26] I. Zacharov, R. Arslanov, M. Gunin, D. Stefonshin, S. Pavlov, O. Panarin, A. Maliutin, S. Rykovanov. “Zhores” —new PFlops supercomputer for data-driven modeling, machine learning and artificial intelligence installed in SkolkovoInstitute of Science and Technology. Preprint arXiv:1902.07490, 2019.[27] C. K. Huang, W. M. Lou, C. J. Tsai, T.-C. Wu, H.-Y. Lin. Mechanical properties of polymer thin film measured by thebulge test. Thin Sol. Films 515, 7222, 200714] I. Ostanin, R. Ballarini, T. Dumitrica. Distinct element method modeling of carbon nanotube bundles with intertubesliding and dissipation. Appl. Mech. 81(6), 061004, 2014.[15] Y. Wang, C. Gaidau, I. Ostanin, T. Dumitrica. Ring windings from single-wall carbon nanotubes: A distinct elementmethod study. Appl. Phys. Lett. 103 (18), 183902, 2013.[16] Y. Wang, M. Semler, I. Ostanin, E. Hobbie, T. Dumitrica. Rings and rackets from single-wall carbon nanotubes: mani-festations of mesoscale mechanics. Soft Matter 10 (43), 8635-8640, 2014.[17] I. Ostanin, R. Ballarini, T. Dumitrica. Rings and rackets from single-wall carbon nanotubes: manifestations of mesoscalemechanics. Journ. Mat. Res. 30(1), 19, 2015.[18] Y. Wang, I. Ostanin, C. Gaidau, T. Dumitrica. Twisting carbon nanotube ropes with the mesoscopic distinct elementmethod: Geometry, packing, and nanomechanics. Langmuir, 31(45), 12323, 2015.[19] I. Ostanin, P. Zhilyaev, V. Petrov, T. Dumitrica, S. Eibl, U. Ruede, V. Kuzkin. Toward large scale modeling of carbonnanotube systems with the mesoscopic distinct element method. Mater. 8(3), 240-245, 2018.[20] T. Preclik, U. Ruede. Ultrascale simulations of non-smooth granular dynamics. Comp. Part. Mech., 2, 173, 2015.[21] Itasca Consulting Group Inc., 2015. PFC3D (Particle Flow Code in Three Dimensions). Version 5.0. Itasca ConsultingGroup Inc., Minneapolis.[22] V. Kuzkin, I. Asonov. Vector-based model of elastic bonds for simulation of granular solids. Phys. Rev. E, 86(5), 051301,2012.[23] V. Kuzkin, A. Krivtsov. Enhanced vector-based model for elastic bonds in solids. Lett. Mat. 7(4), 455, 2017.[24] MPI Forum. MPI: A message-passing interface standard. Technical report, Knoxville, TN, USA, 1994.[25] C. Ericson, Real-time collision detection. CRC Press, 2004.[26] I. Zacharov, R. Arslanov, M. Gunin, D. Stefonshin, S. Pavlov, O. Panarin, A. Maliutin, S. Rykovanov. “Zhores” —new PFlops supercomputer for data-driven modeling, machine learning and artificial intelligence installed in SkolkovoInstitute of Science and Technology. Preprint arXiv:1902.07490, 2019.[27] C. K. Huang, W. M. Lou, C. J. Tsai, T.-C. Wu, H.-Y. Lin. Mechanical properties of polymer thin film measured by thebulge test. Thin Sol. Films 515, 7222, 2007