Inter-molecular interactions in ultrahigh molecular weight polyethylene single crystals

Inter-molecular interactions in ultrahigh molecular weight polyethylene single crystals

Computational Materials Science xxx (xxxx) xxxx Contents lists available at ScienceDirect Computational Materials Science journal homepage: www.else...

5MB Sizes 0 Downloads 33 Views

Computational Materials Science xxx (xxxx) xxxx

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Inter-molecular interactions in ultrahigh molecular weight polyethylene single crystals ⁎



Sanjib C. Chowdhurya, , Subramani Sockalingamb,c, , John W. Gillespie Jr.a,d,e,f a

Center for Composite Materials, University of Delaware, DE, USA McNAIR Center for Aerospace Innovation and Research, University of South Carolina, SC, USA c Department of Mechanical Engineering, University of South Carolina, SC, USA d Department of Materials Science and Engineering, University of Delaware, DE, USA e Department of Mechanical Engineering, University of Delaware, DE, USA f Department of Civil and Environmental Engineering, University of Delaware, DE, USA b

A R T I C LE I N FO

A B S T R A C T

Keywords: Polymer fibers Shear lag Molecular dynamics simulation

This paper investigates the inter-molecular interactions during tensile loading in ultrahigh molecular weight polyethylene (UHWMPE) single crystals at the atomistic scale. Molecular dynamics (MD) simulations of velocity controlled chain pullout are employed to study inter-molecular load transfer mechanisms. The transfer of tensile load is governed by van der Waals forces that dominate the inter-molecular shear interactions. The tensile stress build up occurs over a length of approximately 40c, where c is the lattice constant along the chain axis. Atomistic MD models incorporate the influence of surrounding neighboring atoms. Therefore, a nonlocal shear lag continuum model is developed for the first time to bridge length scales by extending the classical shear lag model of stress transfer in composites. The nonlocal model predictions correlate better with the MD results compared to the classical shear lag formulation. Combining the MD and shear lag model results, a bilinear mode II cohesive traction-separation behavior is identified to describe the inter-molecular interactions of the continuum with interface stiffness (2.38 GPa/nm), peak traction (0.14 GPa) and mode II fracture toughness (17 mJ/m2).

1. Introduction Polymeric fibers such as Kevlar® KM2 and ultrahigh molecular weight polyethylene (UHMWPE) Dyneema® SK76 are used in soldier protection systems [1,2] due to their high specific tensile strength and tensile stiffness. During impact, fibers are subjected to complex dynamic multiaxial loading [3,4]. Axial tensile specific toughness and longitudinal wave speed are important fiber properties contributing to the ballistic performance. Dyneema® SK76 fibers possess a heterogeneous fibrillar structure as shown in Fig. 1. The fibrillar structure is comprised of a complex hierarchy of fibril diameters which range from nanometer scale, to micrometer scale, up to millimeter scale that form under high draw ratios during processing [5,6]. The microfibrils are suggested to be connected by weak tie-chain molecules and van der Waals (vdW) type forces [7,8]. The UHMWPE crystals are present within the fibrils with the long chain molecules aligned in the fiber direction. In order to design and synthesize next generation fibers, it is important to establish a fundamental understanding of load transfer mechanisms under dynamic tensile loading at multiple length scales.



Both inter-molecular and inter-fibrillar interaction [9] mechanisms are thought to play an important role in tensile load transfer within the fiber, which is similar to the collagen fibrils [10,11]. These load transfer mechanisms are somewhat analogous to the classical shear lag [12] stress transfer at higher length scales that occurs between fiber and matrix in fiber-reinforced composite materials. There are several molecular dynamics (MD) simulations in the literature on determining the anisotropic elastic and strength properties, wave propagation and ballistic performance of crystalline and semicrystalline polyethylene [14–17]. Lee and Rutledge [18] and Kim et el. [19] studied plastic deformation of semi-crystalline polyethylene using MD simulations. They reported plastic deformation behavior during fast deformation under uniaxial extension at large strains. They also noted stress–strain responses are similar and nonlinear in transverse shear deformation in both zx and zy (z is the length direction, x- and y- are lateral directions) indicating transverse isotropy. Due to the weakly connected non-bonded vdW inter-chain interaction between finite length chains, most of the effective tensile properties of the crystal and fibrils in reality are controlled by the inter-chain interactions (i.e.

Corresponding authors at: McNAIR Center for Aerospace Innovation and Research, University of South Carolina, SC, USA. (S. Sockalingam). E-mail addresses: [email protected] (S.C. Chowdhury), [email protected] (S. Sockalingam).

https://doi.org/10.1016/j.commatsci.2019.109360 Received 26 May 2019; Received in revised form 15 October 2019; Accepted 17 October 2019 0927-0256/ © 2019 Elsevier B.V. All rights reserved.

Please cite this article as: Sanjib C. Chowdhury, Subramani Sockalingam and John W. Gillespie Jr., Computational Materials Science, https://doi.org/10.1016/j.commatsci.2019.109360

Computational Materials Science xxx (xxxx) xxxx

S.C. Chowdhury, et al.

Fig. 1. Hierarchical structure of UHMWPE fiber (microfibril image from [13]).

heterogeneity [26]. This approach is ideally suited for the present problem of chain pull-out since the interaction zone is comprised of weakly connected non-bonded vdW inter-chain interactions. Open source code LAMMPS [27] is used for MD simulations and the interatomic interactions are modeled with the Adaptive Intermolecular Reactive Empirical Bond Order (AIREBO) force field [28].

shear) versus the tensile response of the chains. However, there are limited studies on understanding the load transfer mechanisms at the molecular-scale and at higher length scales. Heo et al. [20] studied inter-fibril tribological properties of crystalline polyethylene and they found that tribological properties are dependent on the sliding direction with respect to the fibril orientation. The sliding that occurs along the direction of the chains exhibit a lower friction than sliding in an orthogonal direction. Coussens et al. [9] carried out single chain pullout using all-atom MD simulations to understand the molecular mechanism that limits the strength of UHMWPE fibers. They found that the force required to pull a polyethyelene chain is in the order of 1.3 nN which corresponds to approximately 7.0 GPa stress. This stress value is in very good agreement with the ultimate strength of the UHMWPE fibers produced and tested in the laboratories [21]. Moreover, they observed that the mechanism of pull-out involves a characteristic de-bonding zone of about 40C-atoms ( 5.1nm ). Recently, using all-atom MD simulations, O’Connor et al. [22] studied the yielding mechanism of polyethylene considering chain ends in the simulation models. They identified that polyethylene yields by chain slip instead of bond breakage and the ineffective chain end length over which stress transfer occurs is 2.5 nm (approximately 10c, where c = 0.253 nm, is lattice constant representing the spacing of the corresponding adjacent C2H4 groups in the PE chain). They reported that yield stress increases with chain length and saturates at 6.3 GPa agreeing well with experiments [21]. To design high performance ballistic UHMWPE fibers, it is essential to understand the load transfer mechanism at different length scales and develop numerical models capturing these mechanisms to guide the manufacturing process and experimentalists. In this work we study the inter-molecular interactions in polyethylene fibers at the atomistic scale using MD chain pullout models. Inter-molecular shear transfer interactions are studied by integrating a cohesive zone model with a shear lag model [23]. Since molecular models are inherently nonlocal [24], a nonlocal shear lag model (NSLM) is developed for the first time, by extending the classical shear lag model (CSLM) of stress transfer between fiber and matrix in composites. Nonlocal theory incorporates long range interactions such that the material behavior at a point is influenced by the state of all points within the interaction zone in the body [25]. Integral type nonlocal models account for microstructural length scale information through a kernel function [24] and can serve to link the length scales. Nonlocal models also accurately capture wave dispersions due to material

2. MD modeling and simulation of inter-molecular interactions The foundational building block of UHMWPE fibers is the orthorhombic crystalline lattice [29] shown in Fig. 2. This lattice unit cell has dimensions a = 0.74 nm, b = 0.493 nm and c = 0.253 nm with mass density 1.0 g/cc [29]. It contains two full chains – one at the center and one at each of the four corners (one-fourth of the chain length). MD models are prepared by replicating this unit cell with the polymer chains axis (c- dimension) aligned in the X-direction. Five models with length about 10 nm, 15 nm, 20 nm and 25 and 50 nm are considered to verify length effects and non-local shear lag parameters discussed later. Table 1 describes the dimensions and total number of atoms of the four different models. All the model lengths are larger than the dislocation core size or de-bonded length reported in [9,22]. In the chain, there are both bonded (bond/angle/dihedral) interaction and non-bonded vdW interactions. However, in the transverse directions, chains interact with only non-bonded vdW interactions. Non-bonded interactions are much weaker than the bonded interactions. Therefore, polyethylene crystal or the UHMWPE fibers are anisotropic in nature – stiffer and stronger in the longitudinal direction and more compliant and weaker in the transverse directions. Inter-molecular interactions are determined through slab type shear-mode pullout simulations [9,30]. Fig. 3 shows schematic of loading and boundary conditions for the pullout simulations. During pullout, the middle two rows (along the thickness Y-direction) of chains across the entire width (Z-direction) of the model are pulled. Two rows of chains are considered to minimize through-thickness displacement gradient in the pullout chains. Each row has seven chains along the width direction (Z-direction) giving a total of fourteen pullout chains. Chains are pulled out in the X-direction by applying a velocity of 10 m/s [31] to its right end while keeping the other end free to move. The left end of the top and bottom parts are fixed while the right end of these parts are free to move as shown in Fig. 3. These boundary and loading conditions resemble fiber–matrix shear lag model in composites. Unlike 2

Computational Materials Science xxx (xxxx) xxxx

S.C. Chowdhury, et al.

(a) XY (ab) plane view

(b) 3D view

Fig. 2. Crystal unit cell of PE (Atom color: black – carbon, white – hydrogen).

all-carbon and hydrocarbon materials modeling [32,33]. AIREBO force field has bond formation and breaking capability. Fig. 5 shows the pullout force–displacement response for different length models. Pullout force is periodic with period approximately equal to 1c , where c is the lattice constant along the chain axis. The periodic behavior occurs due to the separation and proximity of the adjacent CH2 groups during pull-out. Force is zero or minimum at the equilibrium distance. When the CH2 groups move away from the equilibrium distance, +ve vdW force develops up to a certain displacement then it starts to reduce. We confirm this periodic stick –slip behavior of the pullout chains by monitoring the longitudinal displacement distribution at different pullout displacement shown in Fig. 6 for the 15 nm model. As shown in Fig. 6, at the beginning of simulations, there is no noticeable displacement. However, as pullout continues, displacement develops at the right end (i.e., pulling end) of the chains up to a certain length and then diminishes to zero at the left end. Length of this transition region (i.e., shear lag region) is about 10.0 nm which is close to the length reported by Coussens et al. [9]. Fig. 4 in their paper shows length of the gradient zone is about 8.3 nm (65C-atoms). Shear lag region, where the load transfer occurs, is the distance between the pulling end of the chains where the displacement develops and the point where the displacement diminishes to zero. Depending on the relative displacement between the pullout chain and the adjacent chain, load will be transferred through this region. Shear lag model type approximation has been reported in the literature [34] to study formation of solitons in interacting chains in aligned polyethylene fibers. The length over which tensile stress build up occurs as predicted by MD results are consistent with the soliton width reported in [34]. As the pullout displacement approaches 1/2c , the left end starts to exhibit displacement (red and cyan marker profiles in Fig. 6) and the pullout force becomes maximum

Table 1 Data for MD models. Model Designation

Dimensions (Lx × Ly × Lz), nm

Total Atoms

10 nm 15 nm 20 nm 25 nm 50 nm

10.14 × 7.89 × 5.18 15.20 × 7.89 × 5.18 20.27 × 7.89 × 5.18 25.34 × 7.89 × 5.18 50.68 × 7.89 × 5.18

53,760 80,640 107,520 134,400 268,800

Model Model Model Model Model

force controlled boundary conditions used in [9], we used velocity controlled pullout since it gives aligned pullout ends for all pullout chains. Before the pullout simulation, the model is first equilibrated to the stress free state with two step relaxation. In the first step of relaxation, the model is relaxed with NPT (isothermal-isobaric) ensemble considering periodic boundary conditions in all directions for 100 ps at T = 300 K, P = 1 atm with timestep 0.2 fs. Nose-Hoover thermostat and barostat are used to control temperature and pressure. In the second step of relaxation, the model is relaxed with NVT (canonical) ensemble. In this case, the model is periodic in the transverse and non-periodic in the longitudinal X-direction. Due to non-periodicity in the longitudinal X-direction direction, chain length becomes finite with chain ends. Since chain ends are not embedded as in [22], we did not hydrogenate the carbon atoms in the exposed chain ends which is consistent with [9]. Pullout simulation with NVT ensemble is carried out after the second stage of relaxation. Force and displacement at the right gripped end region of the pullout chains, and positions of atoms in all chains are recorded during the simulations for displacement distribution and load transfer analysis. Fig. 4 shows a 15 nm MD model with chains pulled out. As mentioned before, interatomic interactions are modeled with reactive force field AIREBO [28], which is widely accepted and used for

Fig. 3. Schematic of pullout model. 3

Computational Materials Science xxx (xxxx) xxxx

S.C. Chowdhury, et al.

Fig. 4. MD model with chains pulled out (15 nm model at 0.8 nm pullout displacement).

conformation where there is no displacement gradient (cross marker profile in Fig. 6). After that, displacement gradient again starts to build up (blue diamond marker profile in Fig. 6) and pullout force starts to pick up (second period in Fig. 5). Heo et al. [20] and Coussens et al. [9] also observed this stick–slip behavior in case of inter-fibrillar sliding and single chain pullout simulations. Neighboring adjacent chains also exhibit similar periodic behavior (not shown here). However, there is a difference in the displacement gradient between the pullout chains and the adjacent chains along the chain length (shown later in section 3) through which load is transferred. Displacement distribution of the longer length 50 nm model at different pullout displacements is shown in Fig. A1 in Appendix. For the five lengths considered, there is no noticeable difference in the peak force and there is no chain (or bond) breakage during the pullout. Peak force varies from 1.29 to 1.36 nN/chain with the corresponding pullout stress 6.7 – 7.0 GPa. This stress is in excellent agreement with the results obtained by Coussens et al. [9] and O’Connor et al. [20] from their simulations and experiment [19]. Peak pullout stress is lower than the stress required to break the chain which is about 17.0 GPa [16]. Since the peak pullout stress is limited by the non-bonded inter-chain (i.e., inter-molecular) interaction with the associated inter-molecular debonding and chain end slip, it could not be increased to the chain’s breakage stress unless inter-chains interaction increases. Higher hydrostatic pressure in transverse impact cases will increase this inter-chains interaction. We plan to study the effects of transverse impact on the chain pullout and inter-chain shear load transfer mechanism and will be reported elsewhere. While the chain slab is pulled out in the X-direction, displacement gradient occurs in the transverse Y-direction as shown in Fig. 7 due to inter-molecular interactions. Length of this displacement gradient region (i.e., zone of influence) on either side of the pullout slab is about 1.6 nm (which is in the same order of vdW cutoff distance 1.02 nm). As a first order approximation, we consider this distance as inter-molecular spacing h in the shear lag models discussed later.

Fig. 5. Pullout force–displacement response obtained from MD simulations for 10 nm, 15 nm, 20 nm, 25 nm and 50 nm models. (Force per chain is determined by normalizing the total pullout force with fourteen pullout chains.)

3. Nonlocal shear lag model In this section a nonlocal shear lag model is derived to describe the inherently nonlocal atomistic pullout response and to link the intermolecular interactions to higher length scales. It should be noted in MD simulations, chain ends are pulled at a velocity of 10 m/s inducing spatial variations in strain rate and strain rate dependent nonlinear behavior in shear directions. The average shear strain in the pullout simulations range from 4% to 7%. In this range, MD models predict nonlinear shear behavior [14]. Whereas shear lag models in this study

Fig. 6. Displacement distribution at different pullout displacement obtained from 15-nm MD model.

(red line in Fig. 5). After that, distribution profiles gradually shift upward with pronounced displacement at the left end indicating chain slippage. When the pullout displacement approaches 1c the entire pullout chain slips and the chains’ conformation goes back to the initial 4

Computational Materials Science xxx (xxxx) xxxx

S.C. Chowdhury, et al.

coupled Equation (3) for displacements

u1 = Ax + B + Ce u2 = Ax + B − Ce

Fig. 8. Unit cell.

assume linear elastic and time independent response. Classical shear lag model [11,12] utilizing unit cell in Fig. 8 is first compared to the MD predictions. The unit cell consists of two slabs of thickness t connected by inter-molecular spacing h over a length L. A prescribed displacement u2 (x = L) = εL is applied at the right end of the bottom slab 2 and the left end of the top slab 1 is constrained u1 (x = 0) = 0 in the axial didu (x = 0) = 0 and the right end of slab 1 rection. The left end of slab 2 2 dx du1 (x = 0) dx

= 0 are stress-free. Ef is the axial stiffness of the slabs and G is the shear stiffness of the inter-molecular region (xy). t is taken as equal to the width of two rows of chains being pulled in the MD simulations, Ef and G are computed from the MD simulations of UHMWPE in [14]. The equilibrium equation is: (1)

Shear stress in the inter-molecular region is given in terms of fibril displacements as:

τ=G

u2 − u1 h

σij, j = ρu¨ i , σij (x ) =

Lc2

(4)

2 x / Lc

∫V α (|x′ − x|) Cijkl εkl (x′) dV (x′), εij = 12 (ui,j + uj,i)

where σij is the stress tensor, Cijkl is the elastic modulus tensor, εij is the strain tensor, ui is the displacement vector and α (|x ′ − x|) is the nonlocal kernel function of Euclidean distance |x ′ − x|. The nonlocal stress gradient theory (1 − ℓ2∇2 ) σij = Cijkl εkl can be derived with a kernel function of the form α (|x ′ − x|) = (2π ℓ2)−1K 0 (|x ′ − x|/ℓ) where K0 is the modified Bessel function and ℓ is the nonlocal internal length scale parameter. The properties of the kernel are such that it reverts to the Dirac delta function in the limit of vanishing internal characteristic length leading to classical elasticity. The nonlocal stress gradient constitutive relation for shear stress is expressed as

(2)

d 2u1 = u2 − u1 dx 2

− De−

2 x / Lc

(5)

Combining Eqs. (1) and (2) and expressing the normal stresses in terms of displacement, the governing equations are:

− Lc2

2x Lc

+ De−

where A, B, C and D are constants determined using the boundary conditions in Fig. 8. Typical displacement, shear stress and normal stress distribution (normalized by the maximum values) in the unit cell predicted by CSLM is shown in Fig. 9. The displacements and stress distributions are symmetric about x = L/2 as expected. The transfer length over which the normal stress in the slabs build up and the interfacial shear stress decreases to almost zero. Slab aspect ratio is deL fined as, s = t . The characteristic length (Lc) is about 21.0 slab thickness for elastic properties used in Table 2. For a given slab thickness, the slab is long enough for the tensile stress build up at higher aspect ratio [35]. It should be noted that MD models predict oscillatory stick–slip type pullout response with multiple debond cycles over the model length. In order to describe the inter-molecular interactions using a cohesive zone traction separation behavior (discussed later in section 4), one debond cycle is considered in the shear lag analysis. Pullout displacements (u2 at x = L) at peak forces from the MD models are used as applied loading in the shear lag models. The displacement distributions at peak loads predicted by MD models are compared to the classical shear lag model in Fig. 10. The input properties obtained from MD models [14] are used in the CSLM are shown in Table 2. The displacement distributions predicted by CSLM, are in general higher than the MD predictions. Regardless of the t, h, and G values classical shear lag displacement distributions does not correlate well with MD predictions. The classical shear lag formulation does not consider any length scale effects and is considered a ‘local’ model. However, atomistic MD models inherently consider the influence of surrounding neighboring atoms, that is, nano-scale interactions. Due to the inherent nonlocal nature of MD solution, the CSLM predictions do not correlate well with MD results. Therefore, a nonlocal shear lag model is developed by extending the CSLM considering the unit cell shown in Fig. 8. While classical continuum models do not allow intrinsic size dependence, MD simulations beyond a few nanometers are computationally very expensive. Nonlocal continuum theory [36] is an integral approach that accounts for characteristic length scale effects where the stress at a material point is a function of strain within a finite volume surrounding that point through a non-local kernel. The non-local governing equations are:

Fig. 7. Variation of longitudinal displacement along transverse direction at x = 10nm in the 15 nm model.

dσ1 dσ τ =− 2 =− dx dx t

2 x / Lc

(3)

d 2u2 = u2 − u1 dx 2 E th

f where Lc = is the characteristic length over which load transfer G occurs. Lc = 5.20 nm for the properties shown in Table 2. Solving the

τ − ℓ2

Table 2 Model values [14]

d 2τ u − u1 =G 2 dx 2 h

(6)

Combining Eqs. (1) and (6), the classical SLM Eq. (3) for displacements are modified as

t(nm)

Ef(GPa)

G (GPa)

h (nm)

0.2465

257.0

3.81

1.60

d 4u d 2u1 ⎤ Lc2 ⎡ℓ2 41 − = u2 − u1 ⎢ dx ⎥ dx 2 ⎦ ⎣ 5

(7)

Computational Materials Science xxx (xxxx) xxxx

S.C. Chowdhury, et al.

Fig. 9. Classical shear lag model (a) slab displacement (b) slab stress distribution.

d 4u d 2u2 ⎤ Lc2 ⎡−ℓ2 42 + = u2 − u1 ⎢ ⎥ dx dx 2 ⎦ ⎣

∊2

Change of variables given by

U = u2 − u1; V = u2 + u1; X = x / Lc ∊ = ℓ/ Lc

d 4U d 2U − = −2U dX 4 dX 2

(9)

The general solution to the coupled fourth order differential equations is given in Eq. (10). Here, the nonlocal effect is introduced

(8)

Fig. 10. CSLM comparison of displacements (a) 10 nm (b) 15 nm (c) 20 nm (d) 25 nm. 6

Computational Materials Science xxx (xxxx) xxxx

S.C. Chowdhury, et al.

Fig. 11. NSLM (a) slab displacement (b) slab shear stress distribution.

Fig. 12. NSLM stress distribution (a) normal stress slab 1; (b) normal stress slab 2; and (c) effective stress.

7

Computational Materials Science xxx (xxxx) xxxx

S.C. Chowdhury, et al.

Fig. 13. NSLM comparison of displacements (a) 10 nm (b) 15 nm (c) 20 nm (d) 25 nm.

Fig. 14. Comparison of axial stresses (a) classical shear lag (b) nonlocal shear lag (average MD data at 0.2 nm displacement are obtained from twenty sets of data at 0.0002 nm interval around 0.2 nm displacement).

8

Computational Materials Science xxx (xxxx) xxxx

S.C. Chowdhury, et al.

Fig. 15. MD results (a) force and energy evolution (b) energy as a function of MD model length.

through ∊ = ℓ/ Lc and the constants corresponding to e x / ∊ Lc terms are taken as zero for boundedness.

u1 = Be−x / ∊ Lc +

Cx Lc

u2 = Be−x / ∊ Lc +

+ D + Ge

Cx Lc

2 x / Lc

+ D − Ge

2x Lc

+ He−

− He−

by the neighboring atoms during the pullout process. The nonlocal displacement distributions correlates much better to the MD predictions as shown in Fig. 13 for 10 nm, 15 nm, 20 nm and 25 nm models. The results are found to be insensitive to nonlocal parameter greater than the length of the model, l > L . Fig. 14 compares the axial stresses in the top and pull out slab predicted by NSLM to the MD virial stress calculations for the 25 nm model at an applied displacement of 0.20 nm near the complete chain slip. To calculate the virial stress, pullout slab and the adjacent slab are discretized into several regions along the X-direction. Dimensions of each region are 0.5 × 0.5 × 5.2 nm. The NSLM predictions correlate somewhat better to the MD results compared to the classical shear lag model. The differences in the results are attributed to the strain rate dependent nonlinear behavior not accounted in the shear lag models.

2 x / Lc 2 x / Lc

(10)

The constants B, C, D, G and H are determined using the four boundary conditions in Fig. 8. An additional boundary condition for the d2u

fifth constant is taken as 21 (x = 0) = 0 . Fig. 11 shows the effect of dx nonlocal parameter (∊) on the slab displacement and shear stress distributions in the unit cell predicted by NSLM for a 25 nm model. For clarity, displacement distributions corresponding to nonlocal parameter bounds l = L and l ≈ 0 are shown in Fig. 11(a). When the nonlocal parameter ∊ approaches zero, the length scale effects vanish leading to the classical solution as shown in Fig. 11(b). As l increases, the magnitude of displacement distribution becomes smaller, especially at the slab ends and are no longer symmetric. Consequently, the magnitude of shear stress distribution is lower at the left end and higher at the right end of the unit cell compared to the classical shear lag model. The results tend to converge as the nonlocal parameter approaches the length of the model. Fig. 12 shows the normal stress distribution in the slabs 1 and 2; and effective stress distributions in the unit cell predicted by NSLM for a 25 nm model. The effective stress in the unit cell is defined as the average of stresses in the fibrils σeff = (σ1 + σ2)/2 . Again, when the nonlocal parameter ∊ approaches zero, the length scale effects vanish leading to the classical solution as shown in Fig. 12. As ∊ increases, the normal stress distribution in fibril 1 becomes smaller, especially at the left fibril end that is constrained and are no longer symmetric. The normal stress distribution in fibril 2, is lower at the left end and higher at the right end that is pulled. The nonlocal effective stress distribution in the unit cell is nonlinear compared to a constant effective stress predicted by local shear lag model as shown in Fig. 12(c). These results are qualitatively consistent with nonlocal theory and models which regularizes stress fields around discontinuities such as crack tips [37] and holes [38]. The nonlocal parameter is determined by matching dispersion curves of nonlocal theory with dispersion curves of atomistic models when linking atomistic to continuum length scales [36]. The tensile stress build-up is found to occur over a length of approximately 40c (10 nm). Therefore, to describe MD length scale response, the nonlocal parameter, l = 40c is taken as the length of tensile stress build-up. This is reasonable, since this length of unit cell is engaged and is influenced

4. Discussion It is well known that fibril-scale structure and properties (fiber morphology) can affect the fiber-scale mechanical properties [1]. Fibrilscale finite element (FE) modeling approach [39] is attractive to explore the structure–property relationships of UHMWPE fibers. However, it is experimentally challenging to measure the nanoscale fibril properties and inter-fibrillar interactions needed as input for the fibrilscale FE models. The atomistic results from this study can serve as input in lieu of the experiments. Force-displacement of one debond cycle predicted by 15 nm MD model is shown in Fig. 15(a). Total energy (in units of aJ) absorbed in one cycle is obtained by integrating force–displacement. The total energy increases with increase in model length (Fig. 15(b)). This energy includes energy due to axial stretching, inter-molecular shearing within the chains and chain slipping. Energy associated with inter-molecular interactions can be represented using a cohesive zone model. A cohesive zone traction separation behavior for sliding Mode II can be derived by integrating cohesive zone model [40] with shear lag analysis [23] of MD pullout results. One debond cycle of the atomistic pullout response is considered to identify the cohesive zone model parameters. The constitutive response of the cohesive zone model is represented using a traction-separation behavior. The parameters that describe the traction-separation behavior are peak traction (τmax), relative displacement at damage onset (δII0 ), critical failure displacement (δIIF ) and fracture toughness (GIIC). Shear stress calculated using Equation (2) is compared in Fig. 16(a). The shear stresses predicted by the nonlocal shear lag model show a much better correlation to the atomistic model. The peak shear stress at the pullout end at an applied strain of 0.8% (strain at which complete 9

Computational Materials Science xxx (xxxx) xxxx

S.C. Chowdhury, et al.

Fig. 16. Inter-molecular interactions (a) shear stress (b) peak traction (c) traction-separation behavior.

the influence of key parameters such as fibril radius, fibril length, fibril stiffness and inter-fibrillar spacing on the fiber tensile properties leading to computational materials design. As shown in Fig. 17(a) (25 nm model), the more intimate the chains are the higher is the efficiency of load transfer. When the chains are far apart, there is less interactions and may lead to premature fiber failure. This loss of chain interactions is proposed a potential fiber failure mechanism in Werff et al. [43]. Similarly, the load transfer efficiency increases with higher inter-molecular shear stiffness as shown in Fig. 17(b). Fig. 17(c) shows the influence of slab thickness on the effective stress–strain. For a given length, the aspect ratio (L/t) increases with decrease in the slab thickness promoting more efficient tensile stress build up. Therefore, it seen that the structural factors including intermolecular spacing, chain aspect ratio and inter-molecular shear stiffness significantly affect the mechanical properties of the fiber. Fig. 18(a) compares the effective strain distribution in the unit cell predicted by local and nonlocal shear lag models for 25 nm model. The nonlocal effective strains are lower in magnitude compared to the local solution along the entire length. The nonlocal (long range) effects becomes important for inter-molecular spacing at the lower (nanometer)

pullout occurs in MD models) is plotted as a function of the model length shown in Fig. 16(b) and is found to plateau with increasing model length. A bilinear traction-separation is considered in this study as shown in Fig. 15(c). The peak traction τmax is taken as the plateau shear stress. The sliding displacement at damage onset is given by δII0 = τmax /(G / h) . Atomistic pullout displacements (Fig. 6) indicate complete debonding at a critical failure displacement δIIF = c . The area under the curve represents the fracture toughness in Mode II and is computed to be equal to approximately 17 mJ/m2. It should be noted that surface free energy during shearing for polyethylene chains is estimated to be 71 mJ/m2 in [41] and the activation energy for pure chain slip is reported to be 10–20 mJ/m2 [42]. The fracture energy estimated from MD pullout simulation is comparable to the activation energy for polyethylene chain slip in shearing reported in [40]. Therefore, similar to the slip-activation, chain pullout is governed by vdW interactions. Fibril effective tensile stress–strain behavior can be derived using the atomistically informed validated nonlocal shear lag model as shown in Fig. 17. Effective strain is defined as the average of the strains in the fibrils εeff = (u1 + u2)/2L . This approach would also allow quantifying 10

Computational Materials Science xxx (xxxx) xxxx

S.C. Chowdhury, et al.

Fig. 17. Nonlocal effective stress (a) effect of inter-molecular spacing (b) effect of inter- molecular shear stiffness (c) effect of slab thickness.

Fig. 18. Nonlocal effects (a) effective strain distribution (b) average effective strain ratio.

with increase in the nonlocal parameter ℓ .

length scales [44] as shown by the average effective strain ratio between nonlocal and local solution (εeff , nonlocal/ εeff , local ) in Fig. 18(b). That is, the material properties are size dependent. The small scale effects decrease with increase in the inter-molecular spacing h and increase 11

Computational Materials Science xxx (xxxx) xxxx

S.C. Chowdhury, et al.

0022. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the Army Research Laboratory or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation herein. The author SS wish to acknowledge the startup funding from the University of South Carolina. The author SS wish to acknowledge Dr. Valery Roy of the University of Delaware for useful discussions and help with solving the differential equations. The author SCC would like to thank Mr. Ethan A. Wise of the University of Delaware for his help in MD simulations data reduction. Data Availability The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study. Fig. A1. Displacement distribution at different pullout displacement obtained from 50 nm model.

References

5. Conclusions

[1] S. Sockalingam, S.C. Chowdhury, J.W. Gillespie, M. Keefe, Recent advances in modeling and experiments of Kevlar ballistic fibrils, fibers, yarns and flexible woven textile fabrics–a review, Text. Res. J. 40517516646039 (2016). [2] K. Krishnan, S. Sockalingam, S. Bansal, S.D. Rajan, Numerical simulation of ceramic composite armor subjected to ballistic impact, Compos. Part B Eng. 41 (2010), https://doi.org/10.1016/j.compositesb.2010.10.001. [3] S. Sockalingam, F.D. Thomas, D. Casem, J.W. Gillespie, T. Weerasooriya, Failure of Dyneema® SK76 single fiber under multiaxial transverse loading, 4051751879865, Text. Res. J. (2018), https://doi.org/10.1177/0040517518798653. [4] S. Sockalingam, J.W. Gillespie, M. Keefe, Influence of multiaxial loading on the failure of Kevlar KM2 single fiber, Text. Res. J. 88 (2018) 483–498, https://doi.org/ 10.1177/0040517516681961. [5] Y. Ohta, H. Murase, T. Hashimoto, Structural development of ultra-high strength polyethylene fibers: Transformation from kebabs to shishs through hot-drawing process of gel-spun fibers, J. Polym. Sci., Part B: Polym. Phys. 48 (2010) 1861–1872. [6] S. Gakkai, High-performance and specialty fibers: Concepts, technology and modern applications of man-made fibers for the future. Tokyo, Japan: Springer; 2016. doi:10.1007/978-4-431-55203-1. [7] P.B. McDaniel, S. Sockalingam, J.M. Deitzel, J.W. Gillespie, M. Keefe, T.A. Bogetti, et al., The effect of fiber meso/nanostructure on the transverse compression response of ballistic fibers, Compos. Part A Appl. Sci. Manuf. 94 (2017), https://doi. org/10.1016/j.compositesa.2016.12.003. [8] S. Sockalingam, D. Casem, T. Weerasooriya, P. McDaniel, J. Gillespie, Experimental Investigation of the High Strain Rate Transverse Compression Behavior of Ballistic Single Fibers, J. Dyn. Behav. Mater. 3 (2017) 474–484, https://doi.org/10.1007/ s40870-017-0126-2. [9] B. Coussens, H. Werff van der, Molecular modeling of chain pullout and scission for oriented polyethylene, 15th Int. Conf. Deform. Yield Fract. Polym., Kerkrade (2012) 11–14. [10] S.E. Szczesny, D.M. Elliott, Interfibrillar shear stress is the loading mechanism of collagen fibrils in tendon, Acta Biomater. 10 (2014) 2582–2590. [11] X. Wei, M. Naraghi, H.D. Espinosa, Optimal length scales emerging from shear load transfer in natural materials: application to carbon-based nanocomposite design, ACS Nano 6 (2012) 2333–2344. [12] H.L. Cox, The elasticity and strength of paper and other fibrous materials, J. Appl. Phys. 3 (1952) 72–79. [13] K.E. Strawhecker, E.J. Sandoz-Rosado, T.A. Stockdale, E.D. Laird, Interior morphology of high-performance polyethylene fibers revealed by modulus mapping, Polym (United Kingdom) 103 (2016) 224–232, https://doi.org/10.1016/j.polymer. 2016.09.062. [14] S.C. Chowdhury, J. Staniszewski, E.M. Martz, R.H. Ganesh, S. Sockalingam, B.Z. Haque, et al., A computational approach for linking molecular dynamics to finite element simulation of polymer chains in polyethylene fibers, Proc. Am. Soc. Compos. - 30th Tech. Conf. ACS (2015, 2015.). [15] H. Dong, Z. Wang, T.C. O’Connor, A. Azoug, M.O. Robbins, T.D. Nguyen, Micromechanical models for the stiffness and strength of UHMWPE macrofibrils, J. Mech. Phys. Solids 116 (2018) 70–98, https://doi.org/10.1016/j.jmps.2018.03. 015. [16] S.C. Chowdhury, S. Sockalingam, J.W. Gillespie, Molecular Dynamics Modeling of the Effect of Axial and Transverse Compression on the Residual Tensile Properties of Ballistic Fiber, Fibers (2017). [17] R.M. Elder, T.C. O’Connor, T.L. Chantawansri, Y.R. Sliozberg, T.W. Sirk, I.-C. Yeh, et al., Shock-wave propagation and reflection in semicrystalline polyethylene: a molecular-level investigation, Phys. Rev. Mater. 1 (2017) 43606, https://doi.org/ 10.1103/PhysRevMaterials. 1.043606. [18] S. Lee, G.C. Rutledge, Plastic Deformation of Semicrystalline Polyethylene by Molecular Simulation, Macromolecules 44 (2011) 3096–3108, https://doi.org/10. 1021/ma1026115. [19] J.M. Kim, R. Locker, G.C. Rutledge, Plastic Deformation of Semicrystalline Polyethylene under Extension, Compression, and Shear Using Molecular Dynamics

This paper investigates the inter-molecular interactions in tensile load transfer in ultrahigh molecular weight polyethylene (UHWMPE) single crystals at the atomistic scale. Molecular dynamics (MD) simulations of displacement controlled chain pullout are employed to study the load transfer mechanisms. The transfer of tensile load is found to follow a shear lag model and occur by means of van der Waals force dominated inter-molecular shear interactions. The tensile stress build up occurs over a length of approximately 40c, where c is the lattice constant along chain axis. Atomistic MD models inherently consider the influence of surrounding neighboring atoms. Therefore, a nonlocal shear lag model is developed for the first time by extending the classical shear lag model of stress transfer in composites. The nonlocal model predictions correlate somewhat better with the MD results compared to the classical shear lag formulation. The shear lag models may be improved by incorporating strain rate dependent and nonlinear shear which is a topic for future study. Combining the MD and shear lag model results, a bilinear cohesive traction-separation behavior is identified to describe the inter-molecular interactions with interface stiffness (2.38 GPa/nm), peak traction (0.14 GPa) and mode II fracture toughness (17 mJ/m2). This model is a first step towards developing a combined shear lag-cohesive zone finite element model for inter-fibrillar interactions and an integrated statistical shear lag model for developing a failure model of the fiber in the future. The model can also be extended to include dynamic loading with equilibrium equation. CRediT authorship contribution statement Sanjib C. Chowdhury: Formal analysis, Investigation, Methodology, Validation, Visualization, Writing - original draft. Subramani Sockalingam: Conceptualization, Formal analysis, Investigation, Methodology, Validation, Visualization, Writing - original draft. John W. Gillespie Jr.: Methodology, Supervision, Writing review & editing, Funding acquisition. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements Research was sponsored by the Army Research Laboratory and was accomplished under Cooperative Agreement Number W911NF-12-212

Computational Materials Science xxx (xxxx) xxxx

S.C. Chowdhury, et al.

[20]

[21]

[22]

[23]

[24]

[25] [26]

[27] [28]

[29] [30]

[31]

[32]

Simulation, Macromolecules 47 (2014) 2515–2528, https://doi.org/10.1021/ ma402297a. S.J. Heo, I. Jang, P.R. Barry, S.R. Phillpot, S.S. Perry, W.G. Sawyer, et al., Effect of the sliding orientation on the tribological properties of polyethylene in molecular dynamics simulations, J. Appl. Phys. (2008;103.), https://doi.org/10.1063/1. 2900884. W. Hoogsteen, H. Kormelink, G. Eshuis, G. ten Brinke, A.J. Pennings, Gel-spun polyethylene fibres, J. Mater. Sci. 23 (1988) 3467–3474, https://doi.org/10.1007/ BF00540480. T.C. O’Connor, M.O. Robbins, Chain Ends and the Ultimate Strength of Polyethylene Fibers, ACS Macro Lett. 5 (2016) 263–267, https://doi.org/10.1021/ acsmacrolett.5b00838. G. Guo, Y. Zhu, Cohesive-Shear-Lag Modeling of Interfacial Stress Transfer Between a Monolayer Graphene and a Polymer Substrate, J. Appl. Mech. 82 (2015) 31005, https://doi.org/10.1115/1.4029635. S. Ghosh, A. Kumar, V. Sundararaghavan, A.M. Waas, International Journal of Solids and Structures Non-local modeling of epoxy using an atomistically-informed kernel, Int. J. Solids Struct. 50 (2013) 2837–2845, https://doi.org/10.1016/j. ijsolstr.2013.04.025. A.C. Eringen, Nonlocal Continuum Field Theories Nonlocal Continuum Field Theories. n.d. T. Hui, C. Oskay, International Journal of Solids and Structures A nonlocal homogenization model for wave dispersion in dissipative composite materials, Int. J. Solids Struct. 50 (2013) 38–48, https://doi.org/10.1016/j.ijsolstr.2012.09.007. S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys. 117 (1995) 1–19. J. Stuart Steven, Alan B. Tutein, J.A. Harrison, A reactive potential for hydrocarbons with intermolecular interactions, J. Chem. Phys. 112 (2000) 6472–6486, https://doi.org/10.1063/1.481208. C.W. Bunn, The crystal structure of long-chain normal paraffin hydrocarbons. The “shape” of the CH2 group, Trans. Faraday Soc. 35 (1939) 482–491. S.C. Chowdhury, T. Okabe, Computer simulation of carbon nanotube pull-out from polymer by the molecular dynamics method, Compos. Part A Appl. Sci. Manuf. 38 (2007) 747–754, https://doi.org/10.1016/j.compositesa.2006.09.011. S.C. Chowdhury, J.W. Gillespie, Silica–silane coupling agent interphase properties using molecular dynamics simulations, J. Mater. Sci. 52 (2017) 12981–12998, https://doi.org/10.1007/s10853-017-1412-z. S.C. Chowdhury, B.Z. Haque, J.W. Gillespie, D.R. Hartman, Molecular simulations

[33]

[34]

[35] [36] [37] [38] [39]

[40]

[41]

[42]

[43]

[44]

13

of pristine and defective carbon nanotubes under monotonic and combined loading, Comput. Mater. Sci. 65 (2012) 133–143, https://doi.org/10.1016/j.commatsci. 2012.07.007. B.Z. Haque, S.C. Chowdhury, J.W. Gillespie, Molecular simulations of stress wave propagation and perforation of graphene sheets under transverse impact, Carbon N Y 102 (2016) 126–140, https://doi.org/10.1016/j.carbon.2016.02.033. A. Hammad, T.D. Swinburne, H. Hasan, S. Del Rosso, L. Iannucci, A.P. Sutton, Theory of the deformation of aligned polyethylene, Proc. R Soc. A Math. Phys. Eng. Sci. 471 (2015) 20150171, https://doi.org/10.1098/rspa.2015.0171. D. Hull, T.W. Clyne, An introduction to composite materials, Cambridge University Press, Second, 1996. A.C. Eringen, Nonlocal continuum field theories, Springer Science & Business Media (2002). S. Ghosh, A. Kumar, V. Sundararaghavan, A.M. Waas, Non-local modeling of epoxy using an atomistically-informed kernel, Int. J. Solids Struct. 50 (2013) 2837–2845. S.A. Silling, Origin and effect of nonlocality in a composite, J. Mech. Mater. Struct. 9 (2014) 245–258. J.M. Staniszewski, S. Sockalingam, T.A. Bogetti, J.W. Gillespie, Modeling the Fibrillation of Kevlar® KM2 Single Fibers Subjected to Transverse Compression, Fibers Polym. 19 (2018) 1479–1489, https://doi.org/10.1007/s12221-018-8127-x. S. Sockalingam, M. Dey, J.W. Gillespie Jr., M. Keefe, Finite element analysis of the microdroplet test method using cohesive zone model of the fiber/matrix interface, Compos. Part A Appl. Sci. Manuf. 56 (2014), https://doi.org/10.1016/j. compositesa.2013.10.021. In-Chul Yeh, Jan W. Andzelm, Calculations of free energy of surface interactions in crystalline polyethylene, J. Chem. Phys. 149 (1) (2018) 014701, https://doi.org/ 10.1063/1.5031026. A.T. Olsson, E. Schr, E. Bergvall, Ab initio and classical atomistic modelling of structure and defects in crystalline orthorhombic polyethylene : Twin boundaries, slip interfaces, and nature of barriers 121 (2017) 234–246, https://doi.org/10. 1016/j.polymer.2017.06.008. Werff, Van Der; Heisserer, U; Coussens, B; Stepanyan, R; Riedel, W; Laessig T., on the Ultimate Potential of High Strength, 30th Int. Symp. Ballist., 2017, p. 1863–1875. B. Arash, Q. Wang, A review on the application of nonlocal elastic models in modeling of carbon nanotubes and graphenes, Comput. Mater. Sci. 51 (2012) 303–313.