- Email: [email protected]

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 ﬁbers 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 inﬂuence of surrounding neighboring atoms. Therefore, a nonlocal shear lag continuum model is developed for the ﬁrst 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 identiﬁed to describe the inter-molecular interactions of the continuum with interface stiﬀness (2.38 GPa/nm), peak traction (0.14 GPa) and mode II fracture toughness (17 mJ/m2).

1. Introduction Polymeric ﬁbers such as Kevlar® KM2 and ultrahigh molecular weight polyethylene (UHMWPE) Dyneema® SK76 are used in soldier protection systems [1,2] due to their high speciﬁc tensile strength and tensile stiﬀness. During impact, ﬁbers are subjected to complex dynamic multiaxial loading [3,4]. Axial tensile speciﬁc toughness and longitudinal wave speed are important ﬁber properties contributing to the ballistic performance. Dyneema® SK76 ﬁbers possess a heterogeneous ﬁbrillar structure as shown in Fig. 1. The ﬁbrillar structure is comprised of a complex hierarchy of ﬁbril diameters which range from nanometer scale, to micrometer scale, up to millimeter scale that form under high draw ratios during processing [5,6]. The microﬁbrils 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 ﬁbrils with the long chain molecules aligned in the ﬁber direction. In order to design and synthesize next generation ﬁbers, 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-ﬁbrillar interaction [9] mechanisms are thought to play an important role in tensile load transfer within the ﬁber, which is similar to the collagen ﬁbrils [10,11]. These load transfer mechanisms are somewhat analogous to the classical shear lag [12] stress transfer at higher length scales that occurs between ﬁber and matrix in ﬁber-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 ﬁnite length chains, most of the eﬀective tensile properties of the crystal and ﬁbrils 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 ﬁber (microﬁbril 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 ﬁeld [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-ﬁbril tribological properties of crystalline polyethylene and they found that tribological properties are dependent on the sliding direction with respect to the ﬁbril 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 ﬁbers. 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 ﬁbers 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 identiﬁed that polyethylene yields by chain slip instead of bond breakage and the ineﬀective 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 ﬁbers, it is essential to understand the load transfer mechanism at diﬀerent 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 ﬁbers 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 ﬁrst time, by extending the classical shear lag model (CSLM) of stress transfer between ﬁber and matrix in composites. Nonlocal theory incorporates long range interactions such that the material behavior at a point is inﬂuenced 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 ﬁbers 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 eﬀects and non-local shear lag parameters discussed later. Table 1 describes the dimensions and total number of atoms of the four diﬀerent 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 ﬁbers are anisotropic in nature – stiﬀer 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 ﬁxed while the right end of these parts are free to move as shown in Fig. 3. These boundary and loading conditions resemble ﬁber–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 ﬁeld has bond formation and breaking capability. Fig. 5 shows the pullout force–displacement response for diﬀerent 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 conﬁrm this periodic stick –slip behavior of the pullout chains by monitoring the longitudinal displacement distribution at diﬀerent 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 ﬁbers. 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 proﬁles 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 ﬁrst equilibrated to the stress free state with two step relaxation. In the ﬁrst 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 ﬁnite 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 ﬁeld 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 proﬁle in Fig. 6). After that, displacement gradient again starts to build up (blue diamond marker proﬁle 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-ﬁbrillar sliding and single chain pullout simulations. Neighboring adjacent chains also exhibit similar periodic behavior (not shown here). However, there is a diﬀerence 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 diﬀerent pullout displacements is shown in Fig. A1 in Appendix. For the ﬁve lengths considered, there is no noticeable diﬀerence 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 eﬀects 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 inﬂuence) on either side of the pullout slab is about 1.6 nm (which is in the same order of vdW cutoﬀ distance 1.02 nm). As a ﬁrst 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 diﬀerent pullout displacement obtained from 15-nm MD model.

(red line in Fig. 5). After that, distribution proﬁles 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 ﬁrst 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 stiﬀness of the slabs and G is the shear stiﬀness 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 ﬁbril 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 modiﬁed 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 ﬁned 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 eﬀects and is considered a ‘local’ model. However, atomistic MD models inherently consider the inﬂuence 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 eﬀects where the stress at a material point is a function of strain within a ﬁnite 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 modiﬁed 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 diﬀerential equations is given in Eq. (10). Here, the nonlocal eﬀect 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) eﬀective 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 diﬀerences 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

ﬁfth constant is taken as 21 (x = 0) = 0 . Fig. 11 shows the eﬀect 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 eﬀects 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 eﬀective stress distributions in the unit cell predicted by NSLM for a 25 nm model. The eﬀective stress in the unit cell is deﬁned as the average of stresses in the ﬁbrils σeff = (σ1 + σ2)/2 . Again, when the nonlocal parameter ∊ approaches zero, the length scale eﬀects vanish leading to the classical solution as shown in Fig. 12. As ∊ increases, the normal stress distribution in ﬁbril 1 becomes smaller, especially at the left ﬁbril end that is constrained and are no longer symmetric. The normal stress distribution in ﬁbril 2, is lower at the left end and higher at the right end that is pulled. The nonlocal eﬀective stress distribution in the unit cell is nonlinear compared to a constant eﬀective 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 ﬁelds 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 inﬂuenced

4. Discussion It is well known that ﬁbril-scale structure and properties (ﬁber morphology) can aﬀect the ﬁber-scale mechanical properties [1]. Fibrilscale ﬁnite element (FE) modeling approach [39] is attractive to explore the structure–property relationships of UHMWPE ﬁbers. However, it is experimentally challenging to measure the nanoscale ﬁbril properties and inter-ﬁbrillar interactions needed as input for the ﬁbrilscale 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 inﬂuence of key parameters such as ﬁbril radius, ﬁbril length, ﬁbril stiﬀness and inter-ﬁbrillar spacing on the ﬁber 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 efﬁciency of load transfer. When the chains are far apart, there is less interactions and may lead to premature ﬁber failure. This loss of chain interactions is proposed a potential ﬁber failure mechanism in Werﬀ et al. [43]. Similarly, the load transfer eﬃciency increases with higher inter-molecular shear stiﬀness as shown in Fig. 17(b). Fig. 17(c) shows the inﬂuence of slab thickness on the eﬀective stress–strain. For a given length, the aspect ratio (L/t) increases with decrease in the slab thickness promoting more eﬃcient tensile stress build up. Therefore, it seen that the structural factors including intermolecular spacing, chain aspect ratio and inter-molecular shear stiﬀness signiﬁcantly aﬀect the mechanical properties of the ﬁber. Fig. 18(a) compares the eﬀective strain distribution in the unit cell predicted by local and nonlocal shear lag models for 25 nm model. The nonlocal eﬀective strains are lower in magnitude compared to the local solution along the entire length. The nonlocal (long range) eﬀects 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 eﬀective tensile stress–strain behavior can be derived using the atomistically informed validated nonlocal shear lag model as shown in Fig. 17. Eﬀective strain is deﬁned as the average of the strains in the ﬁbrils ε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 eﬀective stress (a) eﬀect of inter-molecular spacing (b) eﬀect of inter- molecular shear stiﬀness (c) eﬀect of slab thickness.

Fig. 18. Nonlocal eﬀects (a) eﬀective strain distribution (b) average eﬀective strain ratio.

with increase in the nonlocal parameter ℓ .

length scales [44] as shown by the average eﬀective 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 eﬀects 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 oﬃcial 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 diﬀerential 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 ﬁndings cannot be shared at this time as the data also forms part of an ongoing study. Fig. A1. Displacement distribution at diﬀerent 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 ﬁbrils, ﬁbers, yarns and ﬂexible 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 ﬁber under multiaxial transverse loading, 4051751879865, Text. Res. J. (2018), https://doi.org/10.1177/0040517518798653. [4] S. Sockalingam, J.W. Gillespie, M. Keefe, Inﬂuence of multiaxial loading on the failure of Kevlar KM2 single ﬁber, 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 ﬁbers: Transformation from kebabs to shishs through hot-drawing process of gel-spun ﬁbers, J. Polym. Sci., Part B: Polym. Phys. 48 (2010) 1861–1872. [6] S. Gakkai, High-performance and specialty ﬁbers: Concepts, technology and modern applications of man-made ﬁbers 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 eﬀect of ﬁber meso/nanostructure on the transverse compression response of ballistic ﬁbers, 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. Werﬀ 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, Interﬁbrillar shear stress is the loading mechanism of collagen ﬁbrils 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 ﬁbrous 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 ﬁbers 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 ﬁnite element simulation of polymer chains in polyethylene ﬁbers, 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 stiﬀness and strength of UHMWPE macroﬁbrils, 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 Eﬀect 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 reﬂection 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 inﬂuence of surrounding neighboring atoms. Therefore, a nonlocal shear lag model is developed for the ﬁrst 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 identiﬁed to describe the inter-molecular interactions with interface stiﬀness (2.38 GPa/nm), peak traction (0.14 GPa) and mode II fracture toughness (17 mJ/m2). This model is a ﬁrst step towards developing a combined shear lag-cohesive zone ﬁnite element model for inter-ﬁbrillar interactions and an integrated statistical shear lag model for developing a failure model of the ﬁber 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 ﬁnancial interests or personal relationships that could have appeared to inﬂuence 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., Eﬀect 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 ﬁbres, 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 paraﬃn 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 ﬁeld 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 eﬀect 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 ﬁber/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. Werﬀ, 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.

Copyright © 2021 COEK.INFO. All rights reserved.