Multiscale Modeling of PEEK Using Reactive Molecular Dynamics Modeling and Micromechanics NASA / TM — 2018-219873 May 2018

Polyether ether ketone (PEEK) is a high-performance, semi-crystalline thermoplastic that is used in a wide range of engineering applications, including some structural components of aircraft. The design of new PEEK-based materials requires a precise understanding of the multiscale structure and behavior of semi-crystalline PEEK. Molecular Dynamics (MD) modeling can efficiently predict bulk-level properties of single phase polymers, and micromechanics can be used to homogenize those phases based on the overall polymer microstructure. In this study, MD modeling was used to predict the mechanical properties of the amorphous and crystalline phases of PEEK. The hierarchical microstructure of PEEK, which combines the aforementioned phases, was modeled using a multiscale modeling approach facilitated by NASA’s MSGMC. The bulk mechanical properties of semi-crystalline PEEK predicted using MD modeling and MSGMC agree well with vendor data, thus validating the multiscale modeling approach.


NASA STI Program . . . in Profile
Since its founding, NASA has been dedicated to the advancement of aeronautics and space science. The NASA Scientific and Technical Information (STI) Program plays a key part in helping NASA maintain this important role.
The NASA STI Program operates under the auspices of the Agency Chief Information Officer. It collects, organizes, provides for archiving, and disseminates NASA's STI. The NASA STI Program provides access to the NASA Technical Report Server-Registered (NTRS Reg) and NASA Technical Report Server-Public (NTRS) thus providing one of the largest collections of aeronautical and space science STI in the world. Results are published in both non-NASA channels and by NASA in the NASA STI Report Series, which includes the following report types: • TECHNICAL PUBLICATION. Reports of completed research or a major significant phase of research that present the results of NASA programs and include extensive data or theoretical analysis. Includes compilations of significant scientific and technical data and information deemed to be of continuing reference value. NASA counter-part of peer-reviewed formal professional papers, but has less stringent limitations on manuscript length and extent of graphic presentations.
• TECHNICAL MEMORANDUM. Scientific and technical findings that are preliminary or of specialized interest, e.g., "quick-release" reports, working papers, and bibliographies that contain minimal annotation. Does not contain extensive analysis.
• CONTRACTOR REPORT. Scientific and technical findings by NASA-sponsored contractors and grantees.
• CONFERENCE PUBLICATION. Collected papers from scientific and technical conferences, symposia, seminars, or other meetings sponsored or co-sponsored by NASA.
• SPECIAL PUBLICATION. Scientific, technical, or historical information from NASA programs, projects, and missions, often concerned with subjects having substantial public interest.
• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientific and technical material pertinent to NASA's mission.
For more information about the NASA STI program, see the following: •

Introduction
Polymer matrix composites are increasingly used in place of metals in the aerospace industry due to their comparable mechanical properties and lightweight nature, which increases fuel efficiency and reduces emissions. Polyether ether ketone (PEEK) is a high-performance semi-crystalline thermoplastic ideal for aerospace applications because of its excellent thermo-mechanical properties and resistance to aging. To design improved composite matrix materials incorporating PEEK, a thorough understanding of the molecular structure and thermo-mechanical behavior of PEEK is imperative.
Limited research has been performed to understand the influence of the molecular structure of PEEK on bulk material properties. Talbott et al. demonstrated experimentally that Young's modulus, shear modulus, tensile yield strength, and shear yield strength increase as the relative crystalline content of PEEK increases (Ref. 1). Talbott et al. also determined by x-ray scans that the crystalline structure deteriorates significantly (~10 percent drop in crystallinity) when the applied compression loading increases beyond the yield point (Ref. 1). Despite these initial efforts, more research is needed to provide a better understanding of the effect of molecular structure on bulk thermo-mechanical properties of PEEK. Computational simulation techniques can be used effectively to establish structure-property relationships for PEEK. Specifically, molecular dynamics (MD) simulation techniques are well-suited for predicting bulk level properties of single phase polymers based on molecular structure (Refs. 2 and 3). However, MD simulation techniques alone cannot predict the bulk level properties of semi-crystalline materials because the characteristic length scale of these materials (microns) is orders of magnitude larger than those that can be efficiently simulated with fully atomistic MD models (nanometers). To effectively simulate semi-crystalline PEEK materials, a multi-scale approach can be used in which the amorphous and crystalline phases are modeled separately using MD, and micromechanics is subsequently used to obtain the bulk properties of PEEK using repeating unit cells (RUCs) that represent the spatial arrangement of these phases at larger length scales.
The objective of this study is to demonstrate that MD modeling and Multi-Scale Generalized Method of Cells (MSGMC) (Ref. 4) can be used together to accurately predict the mechanical properties of semicrystalline PEEK. This paper presents the methodology for establishing well-equilibrated amorphous and crystalline PEEK MD models with the reactive force field ReaxFF (Ref. 5), as well as the procedure for predicting the mechanical response of each phase at the molecular level. It also presents the methodology for modeling the semi-crystalline structure of PEEK in the MSGMC module of MAC/GMC using the simulated molecular response as input. Following this, the resulting predicted mechanical properties are compared to experiment for modeling validation.

Hierarchical Structure of Semi-Crystalline PEEK
Multiscale modeling of PEEK is of particular interest since the crystalline phase contains a hierarchical microstructure with various substructures spanning multiple length scales. There are four different structures within the crystalline phase (from largest to smallest in characteristic length): the root network (see Figure 1), the spherulite, the lamella stacks, and the granular crystal block (GCB) (Ref. Inside the GCB are primary crystal blocks bonded together through a secondary structure; Wang et al. stated that the secondary structure is less ordered than the primary crystal blocks. Thus, for modeling purposes herein, the secondary structure has been treated as amorphous PEEK. The primary crystal blocks are 25 nm across and the whole GCB is 120 nm across. Hudson et al. observed lamellar bundles 100 to 200 nm wide in a PEEK/poly(ether imide) blend within spherulites, lending credence to the existence of the GCB structure reported by Wang et al. (Ref. 8). The lamella stacks structure is composed of many GCB's stacked side-by-side. The length of a lamella stack (inside the spherulite) is estimated to be 5 µm, the width (into the page in Figure 2) is estimated to be 0.5 µm, and the thickness is estimated to be 0.12 µm (one GCB thick). The spherulite is composed of many lamella stacks oriented radially, and the core is denser than the outer spherulite region. The diameter of this structure is estimated to be 5 µm. The root network is composed of many spherulites distributed throughout the network, each with their own network of narrow, fibrillary lamellar stacks extending outward from the spherulite (Ref.

Molecular Dynamics Modeling
This section details the MD modeling of amorphous PEEK and crystalline PEEK separately, and then discusses the results of each phase together.

MD Modeling of Amorphous PEEK
The PEEK monomer molecular structure is shown in Figure 3. In addition to the 22 atoms shown in the figure, the hydrogen atoms were also included, for a total of 34 atoms per monomer. Five periodic amorphous PEEK models were created in LAMMPS (Ref. 10), each containing 125 monomers sparsely spaced in the cubic simulation box with side lengths of 100.3 Å, resulting in a mass density of 0.059 g/cc. The fixed-bond OPLS all-atom force field (Ref. 11) was used for the initial set-up of the model. The velocity field of each model at the time of creation was initialized with a different seed value, resulting in five independent model samples. Each model simulation box was densified (using the 'fix deform' LAMMPS command) from a side length of 100.3 to 40 Å over 1 ns followed by densification from 40 to 35.84 Å over 6 ns for a total of 7 ns. The final mass density was 1.30 g/cc, which compares well with the experimentally-observed density of 1.30 g/cc (Ref. 12). During densification, the temperature was set at 300 K using the NVT ensemble with 1 fs time steps. Mass density profiles of each sample were smooth indicating good distribution of PEEK monomers throughout the simulation box.
Simulated polymerization took place immediately after densification over 200 ps with a time step of 0.1 fs and a temperature of 300 K (NVT ensemble). Bonds were created between the carbon (20) and oxygen (14) atoms in Figure 3 with the 'fix bond/create' command every 10 ps with a 50 percent probability and a bond cutoff distance of 7 Å. The final average molecular weight was 2233 ± 287 g/mol.
Once the MD models were fully polymerized, the force field was changed from fixed-bond to ReaxFF for predicting thermo-mechanical properties. To prepare the MD models for the transition to ReaxFF, the temperature of each model was lowered from 300 to 0.3 K (NVT, 1 fs time step) over 2.5 ns. This step was performed using the fixed-bond OPLS force field. The ReaxFF parameter set developed by Liu et al. (Ref. 13) was used for all ReaxFF simulations because of its ability to accurately predict mechanical properties of polymers (Ref. 3). Once the force field was switched to ReaxFF, the temperature of each model was increased from 0.3 to 300 K (NVT, 0.1 fs time step) over 100 ps. Each sample was equilibrated for 2 ns at 300 K and 1 atm using the NPT ensemble with a 0.1 fs time step. Mass density profiles were smooth indicating negligible local residual stresses in the models. The average mass density was 1.30 ± 0.01 g/cc, which compares well with the experimentally-observed mass density of 1.30 g/cc (Ref. 12). The structure of one of the equilibrated MD models is shown in Figure 4.

MD Modeling of Crystalline PEEK
The crystalline PEEK unit cell is orthorhombic with the dimensions a = 7.75 ± 0.01 Å, b = 5.89 ± 0.01 Å, and c = 9.883 ± 0.005 Å. It contains two molecules per unit cell, and has an average molecular weight of 192.2 g/mol and a calculated density of 1.415 g/cc. The unit cell used in this work is shown in Figure 5. Five replicate periodic unit cells were created in LAMMPS, each containing two PEEK monomers (for a total of 68 atoms) in a rectangular prism simulation box with side lengths of 6 Å (x), 6.8 Å (y), and 14.2806 Å (z), resulting in a mass density of 1.643 g/cc. The fixed-bond OPLS all-atom force field was used to initialize the unit cell in LAMMPS; the velocity field of each model at the time of creation was initialized with a different seed value, resulting in five independent model samples. After velocity initialization, each sample's energy was minimized with LAMMPS' built-in minimize function using the following parameters: the conjugate gradient algorithm (Polak-Ribiere version), an energy tolerance of 1×10 -4 (unitless), a force tolerance of 1×10 -6 (kcal/mol-Å), 10,000 maximum iterations, and 1,000,000 maximum force/energy evaluations. The minimized structure was then run at 300 K while holding the volume constant for 3 ns. The time step was 0.1 fs.
Each sample was then equilibrated for 12 ns while setting the temperature at 300 K and pressure at 1 atm using the NPT ensemble with 0.1 fs time steps. Subsequently, each sample was transitioned to ReaxFF for predicting mechanical properties. To prepare the MD models for the transition to ReaxFF, the temperature of each model was lowered from 300 to 1 K over 14 ns while holding the pressure at 1 atm (NPT, 0.1 fs time step). As with the amorphous PEEK, the ReaxFF parameter set developed by Liu et al. (Ref. 13) was used for all ReaxFF simulations. After transition of the force field to ReaxFF, the temperature of each model was increased from 1 to 300 K over 2 ns while holding the pressure at 1 atm (NPT, 0.1 fs time step). Each sample was then equilibrated for 2 ns at 300 K and 1 atm using the NPT ensemble with a 0.1 fs time step. The average mass density was 1.38 ± 0.01 g/cc, which is a little lower than the calculated unit cell density of 1.415 (Ref. 14), but higher than the experimental density of PEEK (1.30 g/cc) (Ref. 12). The structure of one of the equilibrated MD models is shown in Figure 6. The dimensions of the equilibrated MD models were x = 6.36 ± 0.02 Å, y = 7.21 ± 0.02 Å, and z = 15.15 ± 0.05 Å.

Results of MD Modeling
The five relaxed amorphous and crystalline MD samples were subjected to uniaxial tensile strains along the x-, y-, and z-axes for a total of 30 tensile straining simulations (15 amorphous and 15 crystalline). The 'fix deform' command was used to apply a 2×10 8 s -1 strain rate over 1 ns resulting in a 20 percent engineering strain. The NPT ensemble was used to specify a temperature of 300 K and a 1 atm pressure in the transverse directions allowing for Poisson contractions. The time step was set to 0.1 fs. The five relaxed crystalline MD samples were subjected to shearing simulations with respect to the xy-, xz-, and yz-planes after switching to a triclinic simulation box, minimizing, and then equilibrating for 2 ns at 300 K and 1 atm (NPT, 0.1 fs time step). The 'fix deform' command was used to apply a 2×10 8 s -1 strain rate over 500 ps resulting in a 10 percent engineering strain. The NPT ensemble was used to specify a temperature of 300 K and a 1 atm pressure affecting the non-shearing planes. The time step was set to 0.1 fs. The amorphous MD samples were not sheared because their shear modulus could be calculated from the results of the uniaxial simulations by assuming isotropy.
The stress-strain curves from the axial tensile straining simulations of amorphous PEEK were analyzed to give Young's modulus and Poisson's ratio. The "Segmented" package (Ref. 15) in the R statistical language was used to fit a bilinear regression model to the stress-strain curves. Because the predicted stress-strain responses were highly nonlinear and inconsistent with each other beyond the observed yield point (see, for example, Figure 7(a)), only the regions immediately near the observed yield point (Figure 7(b)) were used for locating the exact location of the break point via the R analysis. The determination of the exact breakpoints in the R analysis was based on an initial estimate of the breakpoints in this region of each stress-strain curve. The slope of the first line of the bilinear fit was assumed to be the Young's modulus. The Poisson's ratios were calculated as the average ratio of the transverse to axial strains within the initial slope regime of the bilinear fit. The average predicted amorphous Young's modulus was 3620 ± 943 MPa, where the uncertainty is the standard deviation over the five replicate models. The average simulated amorphous Poisson's ratio was calculated to be 0.39 ± 0.15. Shear modulus was calculated from the Young's modulus and Poisson's ratio; the average predicted shear modulus was 1.28 ± 0.28 GPa. The uncertainty of the shear modulus was calculated as the standard deviation of three shear moduli calculated from the upper uncertainty, average, and lower uncertainty of the predicted Young's modulus and Poisson's ratio. It is important to note that the modeled amorphous PEEK was assumed to be isotropic.
Unlike the isotropic amorphous model, the crystalline PEEK model is orthotropic with nine independent elastic moduli. The crystalline Young's moduli and shear moduli were calculated in a similar manner as the amorphous Young's modulus. Table 1 contains the predicted average amorphous and crystalline elastic moduli. To the authors' knowledge, mechanical properties of 100 percent crystalline PEEK have never been published in the open literature. Though the crystalline PEEK predictions in Table 1 cannot be compared with experimental data, they do exhibit the expected type of anisotropy. That is, the E1 value, calculated from the direction aligned vertically in Figure 6, is significantly higher than the other two axial moduli (140 GPa vs. 4.87 and 7.18 GPa). This behavior can be explained by the presence of a continuous line of covalent bonds in the x1 direction, whereas the behavior along the x2-and x3-axes are dominated by van der Waals bonds. Table 1 shows that the E2 and E3 moduli of the crystalline PEEK are higher than the E modulus of the amorphous PEEK. At first, one would think that the tangled chains of amorphous PEEK would exhibit stronger van der Waals bonding than the crystalline PEEK. However, it is important to recall that a crystal is a highly ordered structure with an optimal packing arrangement. This optimal packing arrangement gives rise to very strong van der Waals bonding hence the E2 and E3 moduli of the crystalline PEEK. Note also that, while the predicted value of ν12 is high, it is within the allowable range for orthotropic materials given the other properties (see Ref. 16), and a positive definite stiffness matrix results.

MSGMC Micromechanics Modeling
The generalized method of cells (GMC) micromechanics theory is an efficient, semi-analytical method that provides the homogenized, nonlinear constitutive response of a heterogeneous material. Its foundations for single scale analysis, along with validation of its results, are well-established in the literature (c.f. Ref. 17). The GMC method considers the material microstructure, on a given length scale, to be periodic. Demonstrative RUCs are shown at various length scales in Figure 8. The RUC is discretized into Nα by Nβ by Nγ subcells, each of which may contain a distinct material. However, as indicated in Figure 8, the unique feature of MSGMC is that the materials occupying the subcells on a given length scale may also be heterogeneous materials, represented by unique RUCs. A given analysis may consist of an arbitrary number of explicit length scales and associated homogenizations, denoted by k (see Figure 8). The highest length scale considered is denoted as level 0, whereas, the current length scale under consideration is length scale i. The GMC theory assumes a first-order displacement field in the subcells at a given scale, resulting in uniform stresses and strains in each subcell. Assuming infinitesimal strains, for the elastic case considered herein, the constitutive equation for the subcells at level i is given by, where σi is the stress tensor, Ci is the stiffness tensor, and εi is the strain tensor. The superscript (αi βi γi) denotes the particular subcell at level i, and αi, βi, and γi indicate the position of the subcell along the x1-, x2-, and x3-directions, respectively. Note that the method is not limited to the linear elastic case considered here, and previous work (Refs. 18 to 20, and 17) have included nonlinearity due to damage and inelasticity. Satisfaction of displacement and traction continuity between subcells in an average (integral) sense, and imposition of periodicity conditions along the repeating unit cell boundaries enable the establishment of a system of linear algebraic equations, which can be solved to determine the strain concentration tensor Note that, within a multiscale analysis, all terms in Equation (2) depend on the location of the level i RUC within the hierarchy of MSGMC. That is, referring to Figure 8, the strain in a given subcell at level k depends on the path taken down the length scales from level 0. Substituting Equation (2) into Equation (1), The RUC average stress tensor is given by, And recognizing that the effective elastic constitutive equation at level i is given by, Equations (5) and (6) indicate that the effective stiffness tensor, * i C , at level i is given by, In MSGMC, the scales are linked by equilibrating the homogenized RUC average stress, strain, and stiffness tensors at Level i to the local stress, strain, and stiffness tensor of a given subcell at Level i-1 (with appropriate transformation to account for the potential coordinate system change from scale to scale). That is, where 2 i T and 4 i T are the appropriate second and fourth order coordinate transformation tensors, respectively. Hence, it is clear that starting with the lowest scale (k) RUC (see Figure 8), whose subcells contain only monolithic materials, the effective stiffness tensor can be calculated using the standard GMC method. This stiffness tensor (after appropriate coordinate transformation) then represents the homogenized material occupying one of the subcells within an RUC at the next higher length scale. Given the transformed effective stiffness tensors of all subcells at this next higher length scale, the effective stiffness tensor of the RUC at this level can be determined. This stiffness tensor can then be transformed and passed along to the next higher length scale, and the process repeats until the highest length scale considered (0) is reached.
As an example, for an MSGMC analysis considering three length scales (0, 1, and 2), the overall effective stiffness tensor can be written using Equations (7) and (8) where a contracted notation has been used for the triple summation at each scale. Note that in Equation (9), the superscript on the bracketed terms indicates that all variables within the brackets are a function of the subcell indices from the next higher length scale (including lower scale dimensions and subcell indices). The intent of this notation is to fully define the location of a subcell at a given scale as one progresses down the length scales. For example, using this notation, the effective stiffness tensor at level 2, from Equation (8), can be written as, as there are distinct * 2 C values for every level 1 subcell, while there are distinct level 1 RUCs present within each level 0 subcell.
Converse to this multiscale homogenization procedure, MSGMC can perform multiscale localization of the stress and strain tensors. While not employed herein because only effective properties are presented for the PEEK material considered, the multiscale localization is needed for inclusion of nonlinearity from damage and inelasticity (see Refs. 19,20,and 17). For the three length scale example, the local strain tensor in an arbitrary lowest scale (level 2) subcell can be written using Equations (2) and (8) as, Again, the superscript on the bracketed terms indicates that all variables within the brackets are a function of the subcell indices from the next higher length scale (including subcell indices). The stress tensor for any subcell at any length scale can be similarly determined through localization, or by simply using the strain tensor, along with the constitutive Equation (1), at the appropriate length scale. Because of its ability to handle multiple length scales in a single analysis, MSGMC is ideal for multiscale modeling of materials such as PEEK that exhibit identifiable microstructures across multiple length scales.

Multiscale Modeling Framework for PEEK
The multiscale micromechanics modeling strategy used for the PEEK microstructure is shown in Figure 9, wherein four length scales (Levels 0 -3) were explicitly considered. Superscripts indicate the associated level and subscripts indicate the specific axis. Previous work incorporated MD with GMC in a hierarchical multiscale framework across three scales to model Epon 862 epoxy with graphene nanoplatelets (GNP) and a carbon fiber composite with GNP/Epon 862 across 4 scales (Refs. 21 and 22. For this study, the crystalline and amorphous regions of Level 3 were assigned the properties in Table 1, calculated using MD. The elastic properties of the GCB, lamella stacks, and spherulite were predicted using MSGMC, as implemented in the MAC/GMC software package developed at the NASA Glenn Research Center (Refs. 23 and 24). In this work, the root network was not modeled, rather, one spherulite structure with surrounding amorphous material was modeled to represent the PEEK structure. A spherulite with 30 percent crystallinity (the typical experimentally observed crystallinity (Ref. 8) was considered to be representative of bulk-level semi-crystalline PEEK.
The PEEK crystallite, at Level 2, was structured as a 5x5 doubly-periodic square RUC with four subcells composed of crystalline PEEK and the balance composed of amorphous PEEK. This structure was based on the four central primary crystals of the granular crystal block (GCB) shown in Figure 2. The dimensions of each row and each column were set such that the crystalline volume fraction was 86.2 percent. This volume fraction was chosen based on the 86.2 percent volume fraction reported for the crystals in this region by Wang et al (Ref. 7).
As shown in Figure 9, the Level 1 lamella stack structure was modeled as a 2x2 doubly-periodic, rectangular RUC with three subcells of amorphous PEEK and the fourth cell composed of the crystallite. The specific width and height of this RUC are dependent upon the crystalline volume fraction; when the crystalline volume fraction changes, the width and height of the RUC change resulting in a change in the aspect ratio. Thus, the aspect ratio (width of RUC divided by height) ranged from 1.47 to 1.73. This structure is henceforth referred to as the plank structure. The 2 1 x crystallite direction was aligned with the 1 3 x plank direction. The dimensions of this structure were calculated automatically based on a given volume fraction of crystallite. Note that the axis superscripts refer to the length scale levels in Figure 9.
The spherulite structure was modeled as a 6x6x6 triply periodic cubic RUC. The spherulite subcells contained plank RUCs with varying volume fractions of crystallite; the volume fraction decreased as the distance from the spherulite center increased to capture the spreading of the planks with distance from the center of the spherulite (see Fig. 9). The range of crystallite volume fractions in the subcells is shown in Table 2. The crystalline volume fraction of all of the subcells averaged together was 30 percent, which is the typical measured crystallinity of PEEK (Ref. 8); thus the spherulite was assumed to be representative of semi-crystalline PEEK. The orientation of planks in spherulite subcells was determined from a vector from the center of the spherulite to the center of each subcell, as shown in Figure 10. Once the plank axial ( 1 1 x ) direction was determined, two random orthogonal vectors, which were orthogonal to the plank axial direction, were established to orient the remaining two plank axes ( 1 2 x and 1 3 x ). A program written in the Python programming language was used to calculate the dimensions, crystalline volume fractions, and orientation vectors of each subcell of the spherulite, and to write the MAC/GMC input file. The two random, orthogonal vectors changed each time the Python program was executed, resulting in slight deviation from isotropy in the mechanical properties of the spherulite. Thus, the results from 50 different runs were averaged to obtain the mechanical properties of the spherulite.

Results and Discussion
To obtain an approximation of the deviation in the predictions, three different input files were run in MAC/GMC: one with the upper standard deviation values from Table 1, one with the average values from  Table 1, and one with the lower standard deviation values from Table 1. All three files (each run 50 times) were used to obtain a range of values for the predicted mechanical properties of semi-crystalline PEEK, shown in Table 3. As can be seen in Table 3, the mechanical properties are nearly isotropic. The experimental data used for comparison to the predictions is the vendor data given by Solvay for their KT-880 NT PEEK (Ref. 12). The predicted elastic properties of semi-crystalline PEEK (30 percent crystalline) are shown alongside the vendor data and relative error in Table 4. The predicted Young's modulus of 4.03 ± 1.16 GPa compares well with the vendor value of 3.70 GPa (crosshead rate of 50 mm/min). The predicted Poisson's ratio of 0.40 ± 0.10 compares well with the vendor data of 0.38. While the vendor did not provide a shear modulus, assuming their PEEK material to be isotropic results in a vendor shear modulus of 1.45 GPa, which compares well with the predicted shear modulus of 1.48 ± 0.33 GPa. As such, the multiscale modeling framework coupling MD with MSGMC has been validated for true prediction of the elastic properties of PEEK without the use of any experimental data as input or calibration of any model parameters.

Conclusion
MD modeling was used to predict the mechanical properties of the amorphous and crystalline phases of PEEK. The hierarchical microstructure of PEEK was modeled using a multiscale modeling approach facilitated by MSGMC incorporating the properties of amorphous and crystalline PEEK. The bulk mechanical properties of semi-crystalline PEEK predicted using MD modeling and MSGMC agree well with vendor data, thus validating the multiscale modeling approach. It should be noted that no experimental input data or calibration of model parameters were used to arrive at the predictions. This validated multiscale modeling approach can be used improve the properties of PEEK through parametric studies; one such study could assess the influence of crystallinity on mechanical properties. This approach can inform the design process to better tailor PEEK materials to meet specific performance requirements.