On Surrogate Modeling of Static Response of AM Short-Fiber Thermoplastics Using Graph Neural Networks

# On Surrogate Modeling of Static Response of AM Short-Fiber Thermoplastics Using Graph Neural Networks
Source: [https://arxiv.org/html/2606.28996](https://arxiv.org/html/2606.28996)
Vipin KumarTrenton M\. RicksSuhasini GururajaSiddhartha SrivastavaGraduate Student, Department of Aerospace Engineering, Auburn University, Auburn, ALSr\. R&D Staff Member, Manufacturing Science Division, ORNL, Knoxville, TNResearch Aerospace Engineer, Multiscale and Multiphysics Modeling Branch, NASA Glenn Research Center, Cleveland, OHProfessor, Department of Aerospace Engineering, Auburn University, Auburn, ALAssistant Professor, Department of Aerospace Engineering, Auburn University, Auburn, AL

###### Abstract

Short\-fiber thermoplastic \(SFT\) composites are increasingly employed in lightweight aerospace and automotive structures owing to their favorable strength\-to\-weight ratio, high production rates, and recyclability\. Unlike continuous\-fiber systems, the mechanical response of SFTs is governed by complex mesoscale interactions among fiber orientation, spatial clustering, and manufacturing\-induced porosity\. These microstructural features exhibit significant spatial variability in additively manufactured components and strongly influence stiffness, damage initiation, and nonlinear deformation\. Although mesoscale finite element \(FE\) models can explicitly resolve such heterogeneity, their application to realistic three\-dimensional microstructures containing thousands of fibers remains computationally intractable\.

In this study, a data\-driven surrogate modeling framework is proposed to efficiently predict the mechanical behavior of additively manufactured, compression\-molded \(AM\-CM\) SFTs\. Three\-dimensional microstructures were reconstructed from high\-resolution micro\-computed tomography \(μ\\mu\-CT\) data and discretized into localized Voronoi\-based cells representing statistically distinct fiber\-interaction neighborhoods\. Each cell was homogenized via nonlinear FE simulations incorporating matrix damage\. The resulting stress\-strain responses were used to train a hybrid Graph Neural Network\-Long Short\-Term Memory \(GNN\-LSTM\) architecture capable of simultaneously encoding microstructural topology and capturing nonlinear, history\-dependent mechanical evolution\.

The trained surrogate accurately predicts the stiffness and stress\-strain behavior of unseen microstructures, with a coefficient of determinationR2≈0\.98R^\{2\}\\approx 0\.98relative to high\-fidelity FE simulations, while achieving more than two\-order\-of\-magnitude reduction in computational cost\. Coupling the framework with experimentally calibrated static damage laws demonstrates that fiber orientation, clustering, and porosity collectively govern the local effective stiffness matrix\. The proposed approach provides a physics\-informed, data\-efficient pathway to identify mechanically weak microstructural cells and accelerate digital\-twin development for SFT components\.

###### keywords:

Short\-fiber thermoplastics \(SFTs\) , Voronoi tessellation , Graph neural networks \(GNN\) , Long short\-term memory \(LSTM\) , Multiscale modeling , Surrogate modeling , Additive manufacturing

## 1Introduction

Short\-fiber thermoplastics \(SFTs\) are increasingly deployed in lightweight aerospace and automotive structures for their high production rates, recyclability, and competitive mechanical properties relative to traditional thermoset\-based composites\[[9](https://arxiv.org/html/2606.28996#bib.bib49),[2](https://arxiv.org/html/2606.28996#bib.bib43),[55](https://arxiv.org/html/2606.28996#bib.bib64)\]\. Recently, hybrid manufacturing that combines additive manufacturing \(AM\) with compression molding \(CM\) has emerged as a promising processing strategy\. This AM\-CM approach mitigates residual porosity, enhances inter\-bead bonding, and improves consolidation in pellet\-extrusion systems\. During AM, intense shear flow within the extrusion nozzle induces preferential fiber alignment along the deposition path\. The subsequent compression molding step, performed above the matrix melt temperature, reduces global porosity by approximately 50%nd strengthens fiber\-matrix interfacial bonding\[[37](https://arxiv.org/html/2606.28996#bib.bib63),[59](https://arxiv.org/html/2606.28996#bib.bib31),[17](https://arxiv.org/html/2606.28996#bib.bib32),[23](https://arxiv.org/html/2606.28996#bib.bib33),[67](https://arxiv.org/html/2606.28996#bib.bib34),[52](https://arxiv.org/html/2606.28996#bib.bib35),[53](https://arxiv.org/html/2606.28996#bib.bib36)\]\.

Despite these processing advantages, AM\-CM SFT components still exhibit pronounced microstructural heterogeneity\[[17](https://arxiv.org/html/2606.28996#bib.bib32),[6](https://arxiv.org/html/2606.28996#bib.bib44),[3](https://arxiv.org/html/2606.28996#bib.bib65),[64](https://arxiv.org/html/2606.28996#bib.bib66)\]\. Unlike continuous\-fiber laminates with layered construction, SFT microstructures exhibit spatial variability in fiber orientation, fiber clustering, and porosity due to complex flow kinematics during deposition and molding\[[75](https://arxiv.org/html/2606.28996#bib.bib42),[17](https://arxiv.org/html/2606.28996#bib.bib32)\]\. From a micromechanics perspective, these variations directly alter local load\-transfer efficiency\. In SFTs, stress redistributes through matrix shear \(shear lag\), making axial stress transfer highly sensitive to fiber proximity, alignment, and termination\. Regions of fiber clustering and fiber\-end neighborhoods introduce stress concentrations, amplify matrix shear stresses, and act as preferred sites for damage initiation under both monotonic and cyclic loading\[[30](https://arxiv.org/html/2606.28996#bib.bib46),[22](https://arxiv.org/html/2606.28996#bib.bib27),[66](https://arxiv.org/html/2606.28996#bib.bib29),[47](https://arxiv.org/html/2606.28996#bib.bib28)\]\.

In our previous work on AM\-CM chopped carbon\-fiber\-reinforced acrylonitrile butadiene styrene \(CF\-ABS\) SFTs, we confirmed a strong correlation between microstructural topology and mechanical performance\[[55](https://arxiv.org/html/2606.28996#bib.bib64)\]\. High\-resolution micro\-computed tomography \(μ\\mu\-CT\) and interrupted fatigue testing revealed that porosity nucleates and grows preferentially near fiber termini and within fiber clusters during cyclic loading\. Fiber ends disrupt matrix continuity and create stress concentrations, while dense fiber clusters restrict matrix flow, forming mechanically weak zones\. As the number of fatigue cycles increases, secondary microstructural evolution—including fiber\-matrix debonding, fiber rotation, and local rearrangement—modifies both the fiber length distribution \(FLD\) and the fiber orientation distribution \(FOD\)\. Quasi\-static tensile tests performed after fatigue interruption reveal reduced stiffness and an earlier onset of nonlinearity relative to the pristine state, confirming that the mechanical response of AM\-CM SFTs is primarily governed by the local topology of fiber interactions\.

Mean\-field homogenization approaches, including Mori–Tanaka and related inclusion\-based models, offer computational efficiency by spatially averaging reinforcement contributions within an equivalent homogeneous matrix\[[50](https://arxiv.org/html/2606.28996#bib.bib50),[58](https://arxiv.org/html/2606.28996#bib.bib2),[33](https://arxiv.org/html/2606.28996#bib.bib3),[10](https://arxiv.org/html/2606.28996#bib.bib4),[73](https://arxiv.org/html/2606.28996#bib.bib5)\]\. Because these formulations rely on spatial averaging, they cannot explicitly represent fiber\-fiber interactions, resin\-rich pockets, or localized damage hotspots that are directly observable inμ\\mu\-CT\-resolved SFT microstructures\[[49](https://arxiv.org/html/2606.28996#bib.bib26),[32](https://arxiv.org/html/2606.28996#bib.bib24),[4](https://arxiv.org/html/2606.28996#bib.bib25)\]\. While such models can provide reasonable estimates of global elastic properties, they remain fundamentally limited in capturing microstructure\-driven nonlinear response and progressive stiffness degradation\[[6](https://arxiv.org/html/2606.28996#bib.bib44),[69](https://arxiv.org/html/2606.28996#bib.bib23),[19](https://arxiv.org/html/2606.28996#bib.bib22)\]\.

At the other end of the spectrum, fully resolved three\-dimensional finite element \(FE\) models based onμ\\mu\-CT data can capture local stress concentrations and load\-transfer pathways\. However, the complexity of realistic SFT microstructures—which feature hundreds to thousands of fibers—necessitates extremely large meshes, resulting in prohibitive computational costs, particularly when accounting for nonlinear plasticity and damage\. This restricts their use for parametric studies, uncertainty quantification, and multiscale modeling\.

Surrogate modeling offers a middle ground between fine\-scale and coarse\-scale approaches by decomposing heterogeneous microstructures into smaller sub\-domains, hereafter referred to as cells, each representing a statistically homogeneous local fiber\-interaction neighborhood that retains key micromechanical features while reducing computational effort\[[72](https://arxiv.org/html/2606.28996#bib.bib21),[60](https://arxiv.org/html/2606.28996#bib.bib20),[45](https://arxiv.org/html/2606.28996#bib.bib19),[28](https://arxiv.org/html/2606.28996#bib.bib18),[35](https://arxiv.org/html/2606.28996#bib.bib16),[20](https://arxiv.org/html/2606.28996#bib.bib17),[61](https://arxiv.org/html/2606.28996#bib.bib48),[14](https://arxiv.org/html/2606.28996#bib.bib41),[16](https://arxiv.org/html/2606.28996#bib.bib37),[65](https://arxiv.org/html/2606.28996#bib.bib47),[26](https://arxiv.org/html/2606.28996#bib.bib45)\]\. However, when each cell still requires nonlinear FE homogenization, a surrogate does not fully close the computational gap\. Scalability remains limited, particularly as finer spatial resolution is required to capture microstructural variability\.

Whereas the surrogate strategies discussed above reduce cost through geometric decomposition while still relying on FE evaluation of each cell, data\-driven surrogate models based on machine learning aim to replace the FE solver itself with a learned mapping from microstructure to response\[[11](https://arxiv.org/html/2606.28996#bib.bib11),[62](https://arxiv.org/html/2606.28996#bib.bib12),[21](https://arxiv.org/html/2606.28996#bib.bib13),[34](https://arxiv.org/html/2606.28996#bib.bib14),[12](https://arxiv.org/html/2606.28996#bib.bib15)\]\. Image\-based convolutional neural networks and graph\-based learning methods have demonstrated promising capabilities for predicting effective properties or stress fields in simplified composite systems\[[24](https://arxiv.org/html/2606.28996#bib.bib38),[74](https://arxiv.org/html/2606.28996#bib.bib39),[16](https://arxiv.org/html/2606.28996#bib.bib37),[70](https://arxiv.org/html/2606.28996#bib.bib40)\]\. However, many existing frameworks rely on synthetic microstructures and do not explicitly represent the physical topology of fiber\-fiber interactions present in realμ\\mu\-CT data\. Furthermore, most existing frameworks are limited to state\-based property prediction\[[10](https://arxiv.org/html/2606.28996#bib.bib4),[76](https://arxiv.org/html/2606.28996#bib.bib10),[61](https://arxiv.org/html/2606.28996#bib.bib48)\], in which the response is assumed to depend only on the current configuration, and do not address history\-dependent behavior where the stress response is governed by the entire deformation history through evolving internal variables\[[27](https://arxiv.org/html/2606.28996#bib.bib1),[13](https://arxiv.org/html/2606.28996#bib.bib8),[31](https://arxiv.org/html/2606.28996#bib.bib9),[76](https://arxiv.org/html/2606.28996#bib.bib10),[48](https://arxiv.org/html/2606.28996#bib.bib68),[25](https://arxiv.org/html/2606.28996#bib.bib69),[40](https://arxiv.org/html/2606.28996#bib.bib70)\]\.

These limitations reveal a critical modeling gap: predictive frameworks are needed that \(i\) preserve experimentally observed microstructural topology, including fiber clustering, fiber\-end effects, and resin\-rich regions; \(ii\) capture history\-dependent nonlinear damage evolution under quasi\-static loading; and \(iii\) remain computationally tractable for coupon\-level analysis\.

In this work, we develop a data\-driven multiscale surrogate modeling framework that operates onμ\\mu\-CT\-resolved microstructural topology, characterized by fiber geometry and their connectivity through proximity\-dependent load\-transfer mechanisms\. The framework learns the mapping from these fiber statistics to the nonlinear, history\-dependent stress–strain response while enforcing physical consistency and enabling computationally efficient predictions\. Figure[1](https://arxiv.org/html/2606.28996#S1.F1)depicts the multiscale modeling workflow for AM\-SFTs\. The AM\-CM process produces a bead\-scale aligned preform characterized via high\-resolutionμ\\mu\-CT imaging\. The reconstructed microstructure is partitioned into mesoscale Voronoi cells that tessellate the domain into subdomains with statistically homogeneous fiber configurations, capturing local variations in clustering and mutual proximity\. Nonlinear FE homogenization is first performed to generate topology\-dependent constitutive data for each cell\. A hybrid Graph Neural Network–Long Short\-Term Memory \(GNN\-LSTM\) surrogate model is then trained to approximate the homogenization operator: fibers are treated as graph nodes, and load\-transfer pathways as weighted edges, thereby preserving the fiber\-to\-fiber interactions, while the LSTM captures path\-dependent plasticity and damage evolution in the matrix\. Individual Voronoi cell responses are subsequently assembled to predict coupon\-scale mechanical behavior while retaining proximity\-driven load\-transfer mechanisms\.

![Refer to caption](https://arxiv.org/html/2606.28996v1/Figure1.png)Figure 1:Overview of the multiscale modeling workflow\. \(a\) Pellet\-extrusion AM, followed by compression molding \(not shown\), produces a CF\-ABS preform\. High\-resolutionμ\\mu\-CT imaging reconstructs the 3D microstructure at the micrometer scale, revealing the matrix, fiber, and pore phases\. \(b\) The reconstructed microstructure is partitioned into mesoscale Voronoi cells using surrogate modeling to allow high\-fidelity coupon\-scale validation\. The workflow spans three length scales: microscale \[μ\\mum\] for high\-fidelity FE homogenization, mesoscale \[mm\] for surrogate\-based cell assembly, and coupon \[cm\] for experimental validation \(adapted from\[[55](https://arxiv.org/html/2606.28996#bib.bib64)\]\)\.
Note:The color map is used for visual differentiation of cells and does not represent fiber count or local fiber concentration\.The primary contribution of this work is the formulation of a topology\-preserving surrogate model that directly linksμ\\mu\-CT\-resolved fiber configurations to nonlinear history\-dependent constitutive response\. Unlike conventional approaches that rely on homogenized statistical descriptors, the proposed topology\-preserving framework explicitly retains local interaction mechanisms governing load transfer and damage evolution\. By replacing high\-fidelity nonlinear FE homogenization with a trained GNN–LSTM surrogate, the framework significantly reduces computational cost while maintaining the micromechanical fidelity required to predict stiffness degradation in AM\-CM SFT systems\.

## 2Materials and Methods

### 2\.1Additive manufacturing\-compression molding process

CF\-ABS panels containing 20 wt\.% chopped carbon fibers were fabricated using a two\-step process developed at Oak Ridge National Laboratory, combining pellet\-extrusion additive manufacturing \(AM\) with subsequent compression molding \(CM\), as described in prior studies\[[54](https://arxiv.org/html/2606.28996#bib.bib67),[55](https://arxiv.org/html/2606.28996#bib.bib64),[56](https://arxiv.org/html/2606.28996#bib.bib62)\]\(see Figure[1](https://arxiv.org/html/2606.28996#S1.F1)a\)\. Unlike conventional filament\-based fused deposition modeling \(FDM\), this approach utilizes a large\-format, pellet\-fed material extrusion system\. During AM, a 20 wt\.% CF\-ABS feedstock, comprising ABS pellets reinforced with chopped high\-tenacity carbon fibers \(Kaltex K20\-HTU\), was deposited using a KUKA robotic extrusion system\.

The extrusion was performed at a screw speed of 1200 rpm, a nozzle traverse speed of 280 mm/s, and a nozzle\-substrate gap of 9\.88 mm, producing single\-bead layers with a height of approximately 6–7 mm\. The extrusion system operated at a nozzle temperature of 523 K, with four internal heating zones progressively increasing from 363 K \(Zone 1\) to 498 K \(Zone 2\), 523 K \(Zone 3\), and 523 K at the nozzle\. This processing condition induces strong shear flow within the nozzle, promoting preferential fiber alignment along the deposition direction and resulting in a transversely anisotropic microstructure at the bead scale\. Fully dense \(100% infill\) additive preforms were produced with dimensions of approximately 450 mm × 300 mm × 6\.5 mm\. Subsequently, compression molding was performed above the ABS melt temperature \(523 K\) for four minutes under a controlled pressure of a 250\-ton hot press \(Carver Press, Model\#\\\#3895\.4NE10000\)\. This step enhanced inter\-bead fusion, reduced global porosity, and improved fiber\-matrix interfacial bonding, yielding consolidated panels with final dimensions of approximately 450 mm × 300 mm × 3\.5 mm\[[56](https://arxiv.org/html/2606.28996#bib.bib62),[55](https://arxiv.org/html/2606.28996#bib.bib64),[54](https://arxiv.org/html/2606.28996#bib.bib67)\]\.

Despite nominally identical processing parameters, the AM\-CM process inherently generates spatially variable microstructures\[[55](https://arxiv.org/html/2606.28996#bib.bib64),[9](https://arxiv.org/html/2606.28996#bib.bib49),[38](https://arxiv.org/html/2606.28996#bib.bib30)\]\. Variations in local flow kinematics, thermal gradients, and consolidation pressure lead to non\-uniform fiber clustering and residual porosity distributions\. While the fabricated panels meet macroscopic quality metrics, localized regions may exhibit high fiber clustering, which can serve as sites of high stress concentration\. From a micromechanics perspective, these features can promote early damage initiation\[[55](https://arxiv.org/html/2606.28996#bib.bib64)\]\.

### 2\.2Ultrasonic inspection, mechanical and microstructural characterization

To investigate the influence of inherent microstructural heterogeneity on mechanical response, three representative specimens were selected from a single AM\-CM panel following ultrasonic C\-scan inspection \(Figure[2](https://arxiv.org/html/2606.28996#S2.F2)a\)\. The selection criterion was based on spatially varying ultrasonic attenuation signatures, which correlate with porosity content and local density variations\.

Each specimen was subjected to an interrupted fatigue protocol to capture microstructural evolution across progressively increasing damage states\. Cyclic loading was applied in stepwise increments at a frequency of1515Hz and a stress ratio ofR=0\.1R=0\.1, with each loading block comprising 7500 cycles at a prescribed stress amplitude\. Starting from 30% of ultimate tensile strength\(σu​t​s\)\(\\sigma\_\{uts\}\), the applied stress was increased incrementally in steps of 10%σu​t​s\\sigma\_\{uts\}, up to specimen failure, as shown in Figure[2](https://arxiv.org/html/2606.28996#S2.F2)b\. At each stress block, the specimen was cyclically loaded for the full 7500 cycles, after which it was interrupted and fully unloaded before either proceeding to the next increment or undergoing characterization\. High\-resolutionμ\\mu\-CT scans were acquired at three key stages: the pristine state before any fatigue loading \(BT\), an intermediate stage below the fatigue limit at50%50\\%σu​t​s\\sigma\_\{uts\}\(BF\), and a higher\-load stage above the fatigue limit at65%65\\%σu​t​s\\sigma\_\{uts\}\(AF\) shown in Figure[2](https://arxiv.org/html/2606.28996#S2.F2)c\. At each interruption, the specimen wasμ\\mu\-CT scanned and subsequently subjected to quasi\-static tensile testing to quantify fatigue\-induced stiffness degradation before reloading for the next increment\.

![Refer to caption](https://arxiv.org/html/2606.28996v1/Figure2.png)Figure 2:Specimen selection strategy and interrupted fatigue experimental protocol\. \(a\) Ultrasonic C\-scan attenuation map of the AM\-CM panel used to select three specimens representing spatially distinct porosity levels; higher attenuation correlates with elevated local pore content\. \(b\) Stepwise cyclic loading protocol showing progressive stress block increments periodically interrupted to allow thermal equilibration before subsequent loading\. and c\) Threeμ\\mu\-CT imaging and quasi\-static tensile testing stages: before fatigue testing \(BT\), below the fatigue limit \(BF\), and above the fatigue limit \(AF\), enabling direct correlation between evolving microstructure and degrading mechanical response within each individual specimen \(adapted from\[[55](https://arxiv.org/html/2606.28996#bib.bib64)\]\)\.μ\\mu\-CT imaging was performed using a Zeiss Xradia 620 Versa system at Auburn University\. Specimens were mounted on a rotary stage perpendicular to the X\-ray source, with a source\-to\-sample distance of 28 mm and a detector distance of 36 mm\. Scans were conducted at 140 kV and 21 W using a 4×\\timesscintillator objective lens with 2×\\times2 binning and an air filter, yielding an isotropic voxel size of approximately 0\.6μ\\mum\. Projection images were reconstructed into three\-dimensional volumes using the system’s built\-in reconstruction software, producing approximately 1000 slices per specimen along the out\-of\-plane direction\.

Quantitative microstructural characterization was performed by extracting geometric descriptors of fibers and pores from the reconstructed volumes\. The fiber aspect ratio \(AR\) is defined as the ratio of fiber length to diameter \(lf/dfl\_\{f\}/d\_\{f\}\)\. Fiber orientation is described by an azimuthal angleϕ\\phifrom the x\-axis and a polar angleθ\\thetain theyy–zzplane; the same convention is applied to pores\. All descriptors were computed within 0\.2 mm sub\-volumes, a window size previously established to ensure convergence of constituent volume fractions and representative characterization of fiber and pore distributions\[[55](https://arxiv.org/html/2606.28996#bib.bib64)\]\. Descriptors were extracted using Dragonfly Pro software \(version 2022\.2\), with further details provided in section[3](https://arxiv.org/html/2606.28996#S3)\.

By repeatedly performingμ\\mu\-CT scanning and quasi\-static testing on the same specimen at each stage, this approach establishes a direct correlation between evolving microstructural features and the corresponding degradation in mechanical properties, while eliminating specimen\-to\-specimen variability\.

### 2\.3Microstructure\-resolved multiscale framework

Using microstructural geometric data obtained fromμ\\mu\-CT, a microstructure\-resolved multiscale framework has been developed to predict the quasi\-static mechanical response of AM\-CM SFTs\. Figure[3](https://arxiv.org/html/2606.28996#S2.F3)presents the overall framework\. To capture interaction physics while maintaining computational tractability at the coupon scale, discrete fiber\-matrix nonlinear FE homogenization is replaced by a physics\-informed surrogate model\. This surrogate efficiently learns the mapping between fiber interaction topology and the strain\-dependent stress response\. Later sections describe each of these aspects in detail\.

![Refer to caption](https://arxiv.org/html/2606.28996v1/Figure3.png)Figure 3:Microstructure reconstruction and mesoscale discretization workflow\. \(a\)μ\\mu\-CT imaging of an AM\-CM SFT specimen yields 3D\-volumetric data segmented into carbon fiber, ABS matrix, and pore phases, from which fiber length, orientation, and clustering statistics are extracted within 0\.2 mm sub\-volumes\. \(b\) A moving sphere traverses the reconstructed volume, grouping fibers whose centerlines intersect each sphere into localized interaction neighborhoods; representative points \(centroids of fiber\-sphere assemblies\) serve as seeds for spatial partitioning\. \(c\) Voronoi tessellation of the mesoscale cell produces a set of non\-overlapping polyhedral cells, each representing a mechanically distinct fiber\-interaction neighborhood for subsequent FE homogenization or surrogate evaluation \(adapted from\[[55](https://arxiv.org/html/2606.28996#bib.bib64)\]\)\.#### 2\.3\.1Extracting fiber features fromμ\\mu\-CT

The volumetric data from high\-resolutionμ\\mu\-CT scans were segmented into three phases—carbon fibers, ABS matrix, and pores—as shown in Figure[3](https://arxiv.org/html/2606.28996#S2.F3)a\. Fiber centerlines were then extracted from the segmented sub\-volume \(Figure[3](https://arxiv.org/html/2606.28996#S2.F3)b\), preserving key geometric features governing mechanical response, including fiber lengths, orientation, clustering, and spatial distribution\. In short\-fiber thermoplastics, stress transfer is governed by shear\-lag mechanisms, making the relative spatial arrangement of fibers critical\. In particular, regions of high fiber density \(clusters\) and low fiber density \(resin\-rich zones\) exhibit fundamentally different load transfer characteristics and damage evolution\. Accurately preserving these spatial heterogeneities is therefore essential for any physically meaningful microstructural representation\.

To capture these localized interactions, a moving\-sphere sampling procedure is employed \(Figure[3](https://arxiv.org/html/2606.28996#S2.F3)b\)\. A sphere of radiusRsR\_\{s\}traverses the reconstructed volume on a regular three\-dimensional Cartesian grid with uniform spacing equal toRsR\_\{s\}in all three directions, and fibers whose centerlines intersect the sphere are grouped into local neighborhoods\. The radiusRsR\_\{s\}thus defines an interaction length scale: a smallerRsR\_\{s\}captures highly localized fiber interactions, while a largerRsR\_\{s\}yields coarser neighborhoods that average over multiple clusters and matrix\-rich regions\. Moreover,RsR\_\{s\}is selected as a global parameter based on the mean fiber lengthlf,m​e​a​nl\_\{f,mean\}obtained from theμ\\mu\-CT data\. Specifically,Rs≳lf,m​e​a​nR\_\{s\}\\gtrsim l\_\{f,mean\}ensures that individual fibers are sufficiently resolved within each sampled neighborhood and that shear\-lag\-based load transfer along the fiber length is not artificially truncated\.

As the sphere traverses the microstructure, the chosen radiusRsR\_\{s\}generates a set of overlapping neighborhoods that capture local fiber interactions\. Due to this overlap, individual fibers may be included in multiple neighborhoods, ensuring continuity of interaction information across adjacent regions\. Each neighborhood is associated with a representative point, defined as the centroid of the combined set comprising the moving sphere center and the centroids of all fibers contained within the sphere\. These representative points provide a discrete sampling of the microstructure and form the basis for subsequent spatial partitioning\. A visualization of the moving\-sphere traversal across a representative microstructure is provided in Supplemental Material[Supplemental Material\-A\. Sphere\-Sliding Visualization](https://arxiv.org/html/2606.28996#Sx4)\.

#### 2\.3\.2Voronoi tessellation of mesoscale domain

The representative points obtained from the moving\-sphere sampling procedure serve as seeds to construct a Voronoi tessellation of the mesoscale domain\. This approach is preferred over alternative discretization methods, such as regular grids,kk\-means clustering, or Octree decomposition, because it preserves the intrinsic structural consistency of the underlying fiber distribution\. First, because the seeds are derived directly from the fiber sampling process, the resulting cells are inherently topology\-adaptive: fiber\-dense regions naturally produce smaller cells, while matrix\-rich regions yield larger cells, eliminating the need for user\-defined mesh scales\. Second, the tessellation provides a complete, non\-overlapping partition of the domain, ensuring volumetric integrity during homogenization\. Consequently, the mesoscale domain is decomposed into a set of non\-overlapping polyhedral cells \(Ωc\\Omega\_\{c\}\) that collectively span the entire domain \(Figure[3](https://arxiv.org/html/2606.28996#S2.F3)c\)\. Third, the equidistant nature of Voronoi boundaries defines a geometrically consistent local neighborhood\. This geometry is directly controlled by the sampling sphere radius \(RsR\_\{s\}\), which dictates the minimum spacing between seeds and thereby influences the spatial distribution and local fiber volume fraction\. For instance, a larger sampling radius \(Rs=180​μ​mR\_\{s\}=180~\\mu\\text\{m\}\) results in a coarser partition with fewer cells compared to a more refined sampling atRs=60​μ​mR\_\{s\}=60~\\mu\\text\{m\}\(see Figure[4](https://arxiv.org/html/2606.28996#S2.F4)\)\. The Python implementation of the moving\-sphere sampling and Voronoi tessellation described here is publicly available at[GitHub](https://github.com/AMPL-Gururaja)\([https://github\.com/AMPL\-Gururaja](https://github.com/AMPL-Gururaja)\)\. The repository includes scripts for sphere traversal, seed generation, Voronoi cell construction, STL export, and fiber volume fraction computation for each discrete cell\.

![Refer to caption](https://arxiv.org/html/2606.28996v1/Figure4.png)Figure 4:Dependence of Voronoi cell morphology on the moving\-sphere sampling radiusRsR\_\{s\}, illustrated for three representative values:RsR\_\{s\}= 60 µm \(left\),RsR\_\{s\}= 120 µm \(center\), andRsR\_\{s\}= 180 µm \(right\)\. AsRsR\_\{s\}increases, the number of seed points decreases, and individual cells grow larger\.
#### 2\.3\.3Discrete\-cell Finite\-Element Homogenization

Each discrete Voronoi cellΩc\\Omega\_\{c\}is treated as a repeating unit cell \(RUC\) with periodic boundary conditions \(PBCs\) over which effective properties are evaluated using FE simulations\. All FE simulations were performed using the commercially available software Abaqus 2022\. Automated Python scripting was employed to generate equivalent cuboidal RUCs from the irregular Voronoi polyhedra through the following mapping procedure\.

First, the bounding box of each Voronoi cell is computed from its vertex coordinates, and an equivalent cuboid is constructed such that its volume matches the Voronoi cell volume\. Second, the fiber centerline coordinates extracted fromμ\\mu\-CT—specifically the start and end points of each fiber segment clipped to the Voronoi cell—are directly transferred into the cuboidal RUC coordinate system by applying a rigid translation that maps the Voronoi centroid to the cuboid centroid, preserving the absolute fiber orientations and inter\-fiber spacing\. This mapping ensures two key microstructural properties are preserved: \(i\) the fiber volume fraction within the Voronoi cell, and \(ii\) the fiber orientation distribution within the cell\.

The fibers and matrix were modeled using three\-dimensional tetrahedral elements \(C3D4\), with PBCs applied on face, edge, and vertex pairs of the cuboidal RUC following the implementation detailed in our prior work\[[5](https://arxiv.org/html/2606.28996#bib.bib6),[55](https://arxiv.org/html/2606.28996#bib.bib64),[63](https://arxiv.org/html/2606.28996#bib.bib7)\]\. To extract the effective response of the RUC, six independent macroscopic strain cases \(three normal and three shear\) are applied sequentially in the form of far\-field displacement fields\. For each loading configuration, the stress and strain fields within the RUC are obtained from finite element simulations\. The effective \(homogenized\) stress response is determined through volume averaging of the local fields at integration points using Equation[1](https://arxiv.org/html/2606.28996#S2.E1):

𝝈¯=1V​∑e=1Ne∑g=1Nint𝝈g,e​Vg,e\\bar\{\\boldsymbol\{\\sigma\}\}=\\frac\{1\}\{V\}\\sum\_\{e=1\}^\{N\_\{e\}\}\\sum\_\{g=1\}^\{N\_\{\\mathrm\{int\}\}\}\\boldsymbol\{\\sigma\}^\{g,e\}\\,V^\{g,e\}\(1\)
where𝝈¯\\bar\{\\boldsymbol\{\\sigma\}\}is the homogenized second\-order stress tensor,NeN\_\{e\}is the total number of elements,NintN\_\{\\mathrm\{int\}\}is the number of integration points per element \(i\.e\. 1 in this case\),𝝈g,e\\boldsymbol\{\\sigma\}^\{g,e\}is the stress tensor at thegg\-th integration point in theee\-th element,Vg,eV^\{g,e\}is the volume \(weight\) of thegg\-th integration point in theee\-th element, andVVis the total volume of the RUC\.

The Voronoi cell\-level effective response is represented by a strain\-dependent constitutive operator \(tangent stiffness evolution\)ℂ​\(𝜺¯\)\\mathbb\{C\}\(\\bar\{\\boldsymbol\{\\varepsilon\}\}\)and the corresponding stress–strain curve\. Mesh convergence studies were conducted until the change in homogenized Young’s modulus between successive refinements fell below 2%\. A converged mesh element size of0\.5​μ​m0\.5~\\mu\\text\{m\}was used in the remainder of the study\. The effective stress–strain response obtained from FE homogenization of each discrete Voronoi cell is subsequently assembled across the full mesoscale Voronoi tessellation through volume\-weighted averaging to predict the coupon\-scale mechanical response, as described in detail in Section[2\.5](https://arxiv.org/html/2606.28996#S2.SS5)\.

#### 2\.3\.4Constitutive Modeling of Constituents

Carbon fibers are modeled as orthotropic linear\-elastic solids using the properties listed in Table[1](https://arxiv.org/html/2606.28996#S2.T1)\. Under the investigated quasi\-static loading regime, post\-mortem scanning electron microscopy \(SEM\) analysis reveals no evidence of fiber fracture\[[55](https://arxiv.org/html/2606.28996#bib.bib64)\]; stiffness degradation is therefore attributed primarily to matrix\-driven degradation mechanisms, and fiber damage has not been accounted for in the current analysis\. Additionally, the fiber\-matrix interface is assumed to be perfectly bonded\.

Table 1:Elastic properties of the carbon fiber and ABS matrix constituents used in the 20 wt\.% CF\-ABS discrete\-cell models\.ConstituentE11E22E33G12G13G23ν12\\nu\_\{12\}ν13\\nu\_\{13\}ν23\\nu\_\{23\}\(GPa\)\(GPa\)\(GPa\)\(GPa\)\(GPa\)\(GPa\)Carbon fibera23014149950\.250\.250\.30ABS matrixb2\.312\.312\.310\.860\.860\.860\.350\.350\.35
- •aCarbon fiber properties adopted from literature\[[71](https://arxiv.org/html/2606.28996#bib.bib54)\]\.
- •bABS matrix elastic modulus and Poisson’s ratio are determined experimentally\.

The baseline elastic properties of the ABS matrix phase were assumed to be isotropic and were experimentally determined in this study by quasi\-static uniaxial tensile testing in accordance with ASTM D638 Type I at 2 mm/min\. Standalone neat ABS specimens were manufactured using the same processing conditions as the composite matrix phase\. Testing was performed at room temperature under a displacement\-controlled loading rate to isolate the baseline elastic modulus \(E=2\.31E=2\.31GPa\) and Poisson’s ratio \(ν=0\.35\\nu=0\.35\), which are assigned as the initial undamaged properties in Table[1](https://arxiv.org/html/2606.28996#S2.T1)\.

Matrix damage: The ABS matrix was modeled as an isotropic elastoplastic material with progressive damage to capture the experimentally observed nonlinearity, near\-peak softening, and stiffness degradation\. All damage model calibrations were performed using true stress–true strain data to ensure consistency under finite deformation\. The initial elastic modulusEEwas obtained via linear regression within the small\-strain regime\. The global yield stressσY\\sigma\_\{Y\}was defined using the experimentally measured 0\.2% offset criterion, averaged across 3 specimens\.

Plastic hardening was represented using a two\-term Voce law\[[7](https://arxiv.org/html/2606.28996#bib.bib51),[18](https://arxiv.org/html/2606.28996#bib.bib53)\]:

σy​\(ep\)=σY\+R1​\(1−e−b1​ep\)\+R2​\(1−e−b2​ep\),\\sigma\_\{y\}\(e\_\{p\}\)=\\sigma\_\{Y\}\+R\_\{1\}\\left\(1\-e^\{\-b\_\{1\}e\_\{p\}\}\\right\)\+R\_\{2\}\\left\(1\-e^\{\-b\_\{2\}e\_\{p\}\}\\right\),\(2\)
whereepe\_\{p\}is the plastic strain andR1,b1,R2,b2R\_\{1\},b\_\{1\},R\_\{2\},b\_\{2\}are global hardening parameters identified by fitting the true stress–strain response up to peak stress while enforcing consistency with the prescribed yield stress \(see Figure[5](https://arxiv.org/html/2606.28996#S2.F5)a\)\. The two\-term form was adopted because ABS, as an amorphous thermoplastic, exhibits a characteristic saturating hardening response in which the hardening rate declines progressively with plastic strain\[[36](https://arxiv.org/html/2606.28996#bib.bib52)\]; a single exponential term is insufficient to simultaneously capture both the early nonlinear hardening and the near\-peak plateau, whereas the two\-term superposition provides the required flexibility without additional physical assumptions\.

![Refer to caption](https://arxiv.org/html/2606.28996v1/Figure5.png)Figure 5:Calibration of the ABS matrix elastoplastic damage constitutive model\. \(a\) Comparison of the calibrated two\-term Voce plasticity plus exponential damage model against true stress\(σ11\)\(\\sigma\_\{11\}\)–true strain\(ε11\)\(\\varepsilon\_\{11\}\)curves from three specimens, demonstrating excellent agreement through the elastic, nonlinear hardening, and post\-peak softening regimes\. \(b\) Evolution of the scalar damage variableDmD\_\{m\}with applied true strain\(ε11\)\(\\varepsilon\_\{11\}\), showing damage initiation near peak stress followed by progressive softening\. \(c\) Equivalent plastic strain\(εp\)\(\\varepsilon\_\{p\}\)at the onset of damage initiation \(matrix failure index\), illustrating the threshold beyond which crack\-band\-regularized softening is activated\.Progressive stiffness degradation observed near and beyond peak stress was captured using a phenomenological scalar continuum damage mechanics \(CDM\) framework\. Within this framework, the scalar damage variableDm∈\[0,1\)D\_\{m\}\\in\[0,1\)acts uniformly on the undamaged elastic stiffness tensorℂm\\mathbb\{C\}\_\{m\}, so that all normal and shear stiffness components are degraded by the same factor\(1−Dm\)\(1\-D\_\{m\}\)\. The nominal Cauchy stress is therefore:

𝝈m=\(1−Dm\)​ℂm:𝜺m\\boldsymbol\{\\sigma\}\_\{m\}=\(1\-D\_\{m\}\)\\,\\mathbb\{C\}\_\{m\}:\\boldsymbol\{\\varepsilon\}\_\{m\}\(3\)
This isotropic stiffness degradation representation is consistent with established CDM formulations for thermoplastic matrices\[[36](https://arxiv.org/html/2606.28996#bib.bib52)\]\.

Damage initiates onceep≥ep​0e\_\{p\}\\geq e\_\{p0\}and evolves according to an exponential softening law \(see Figure[5](https://arxiv.org/html/2606.28996#S2.F5)b and Figure[5](https://arxiv.org/html/2606.28996#S2.F5)c\):

Dm=1−exp⁡\(−⟨ep−ep​0⟩\+ep​f\),D\_\{m\}=1\-\\exp\\\!\\left\(\-\\frac\{\\langle e\_\{p\}\-e\_\{p0\}\\rangle\_\{\+\}\}\{e\_\{pf\}\}\\right\),\(4\)
where the Macaulay bracket⟨⋅⟩\+\\langle\\cdot\\rangle\_\{\+\}ensures damage activates only whenep\>ep​0e\_\{p\}\>e\_\{p0\}\. The damage initiation strainep​0e\_\{p0\}and the characteristic softening strainep​fe\_\{pf\}were both treated as phenomenological parameters identified directly from the post\-peak true stress–strain data of each specimen\. The complete set of calibrated parameters is summarized in Table[2](https://arxiv.org/html/2606.28996#S2.T2)\.

Table 2:Calibrated ABS elastoplastic and damage parameters identified from true stress–strain data\. Global plasticity parameters are common across all specimens; damage parameters are specimen\-specific owing to differing porosity levels\.
#### 2\.3\.5Effective Matrix Stiffness Incorporating Porosity

A two\-step homogenization strategy was employed to incorporate the effect of porosity on the matrix’s effective stiffness \(Figure[6](https://arxiv.org/html/2606.28996#S2.F6)\)\. First, a matrix\-only RUC containingμ\\mu\-CT\-resolved pores were numerically homogenized \(similar to the procedure described in[2\.3\.3](https://arxiv.org/html/2606.28996#S2.SS3.SSS3)\) to obtain an effective porous\-matrix stiffness tensorCmeffC^\{\\text\{eff\}\}\_\{m\}\.Second, this matrix stiffness tensor was assigned to all fiber\-matrix RUC models\. This procedure isolated pore\-induced matrix degradation while allowing the fiber topology to vary spatially between cells\[[1](https://arxiv.org/html/2606.28996#bib.bib61)\]\.

![Refer to caption](https://arxiv.org/html/2606.28996v1/Figure6.png)Figure 6:Figure 6:Two\-step micromechanical homogenization procedure for incorporating matrix porosity effects\. \(a\) Microstructure of SFTs consisting of carbon fibers and ABS matrix with microporosity, as reconstructed fromμ\\mu\-CT imaging\. \(a\) Step 1: a matrix\-only RUC withμ\\mu\-CT\-resolved pore morphology is homogenized under prescribed macroscopic strains to obtain the effective porous\-matrix stiffness tensorCmeffC^\{\\text\{eff\}\}\_\{m\}\(”Homogenized porous matrix”\), \(b\) Step 2:CmeffC^\{\\text\{eff\}\}\_\{m\}is assigned to all fiber\-matrix RUC used in the discrete\-cell FE simulations \(”Homogenized SFTs”\), yielding fatigue\-stage\-resolved constitutive data that reflect measured porosity levels without explicitly meshing pores in every cell\.

### 2\.4Surrogate Modeling Framework

A surrogate constitutive framework was introduced to efficiently approximate the nonlinear, history\-dependent homogenized response of each discrete Voronoi cell\. For each Voronoi cellΩc\\Omega\_\{c\}, the nonlinear finite\-element homogenization defines a path\-dependent constitutive operator in which the homogenized stress at pseudo\-timettdepends on the entire applied macroscopic strain history up tott:

𝝈¯\(Ωc\)​\(t\)=ℋ𝒯Ωc​\(\{𝜺¯​\(s\)\}s≤t\),\\bar\{\\boldsymbol\{\\sigma\}\}^\{\(\\Omega\_\{c\}\)\}\(t\)\\;=\\;\\mathcal\{H\}\_\{\\mathcal\{T\}\_\{\\Omega\_\{c\}\}\}\\\!\\bigl\(\\\{\\bar\{\\boldsymbol\{\\varepsilon\}\}\(s\)\\\}\_\{s\\leq t\}\\bigr\),\(5\)
where\{𝜺¯​\(s\)\}s≤t\\\{\\bar\{\\boldsymbol\{\\varepsilon\}\}\(s\)\\\}\_\{s\\leq t\}denotes the macroscopic strain history up to timett, andℋ𝒯Ωc\\mathcal\{H\}\_\{\\mathcal\{T\}\_\{\\Omega\_\{c\}\}\}is a constitutive operator parameterized by the local fiber\-interaction topology𝒯Ωc\\mathcal\{T\}\_\{\\Omega\_\{c\}\}, comprising fiber geometry, spacing, clustering, and orientation, withinΩc\\Omega\_\{c\}\. The explicit dependence on strain history reflects the path\-dependent nature of the response, which arises from evolving internal variables \(accumulated plastic strainepe\_\{p\}and scalar damageDmD\_\{m\}\) in the elastoplastic\-damage matrix described in Section[2\.3\.4](https://arxiv.org/html/2606.28996#S2.SS3.SSS4)\. The operator is therefore not reducible to a single\-valued function of current strain, and any surrogate intended to replace it must inherit this history\-dependent structure\. As the number of discrete cells in a mesoscale assembly grew to the order of10210^\{2\}–10310^\{3\}, repeated evaluation ofℋ𝒯Ωc\\mathcal\{H\}\_\{\\mathcal\{T\}\_\{\\Omega\_\{c\}\}\}through high\-fidelity nonlinear FE simulations became computationally prohibitive\. A surrogate operatorℋ^𝒯Ωc\\hat\{\\mathcal\{H\}\}\_\{\\mathcal\{T\}\_\{\\Omega\_\{c\}\}\}is therefore constructed such that

𝝈^\(Ωc\)​\(t\)=ℋ^𝒯Ωc​\(\{𝜺¯​\(s\)\}s≤t\)≈ℋ𝒯Ωc​\(\{𝜺¯​\(s\)\}s≤t\)\.\\hat\{\\boldsymbol\{\\sigma\}\}^\{\(\\Omega\_\{c\}\)\}\(t\)\\;=\\;\\hat\{\\mathcal\{H\}\}\_\{\\mathcal\{T\}\_\{\\Omega\_\{c\}\}\}\\\!\\bigl\(\\\{\\bar\{\\boldsymbol\{\\varepsilon\}\}\(s\)\\\}\_\{s\\leq t\}\\bigr\)\\;\\approx\\;\\mathcal\{H\}\_\{\\mathcal\{T\}\_\{\\Omega\_\{c\}\}\}\\\!\\bigl\(\\\{\\bar\{\\boldsymbol\{\\varepsilon\}\}\(s\)\\\}\_\{s\\leq t\}\\bigr\)\.\(6\)The surrogate was trained on FE\-generated stress–strain trajectories that satisfy equilibrium, the Hill–Mandel macrohomogeneity condition, and the two\-term Voce isotropic\-hardening elastoplastic\-damage constitutive law calibrated for the ABS matrix phase \(Section[2\.3\.4](https://arxiv.org/html/2606.28996#S2.SS3.SSS4)\)\. The construction ofℋ^𝒯Ωc\\hat\{\\mathcal\{H\}\}\_\{\\mathcal\{T\}\_\{\\Omega\_\{c\}\}\}was decomposed into two complementary components that mirror the structure of classical history\-dependent constitutive models: \(i\) a GNN that encodes the topology𝒯Ωc\\mathcal\{T\}\_\{\\Omega\_\{c\}\}into a fixed\-dimensional descriptor capturing the cell’s stiffness potential and matrix\-mediated load\-transfer pathways, and \(ii\) a LSTM network that propagates the strain history through learned hidden\- and cell\-state vectors that serve as surrogate internal variables encoding accumulated plasticity and damage\. The GNN therefore parametrizes the topology argument ofℋ^\\hat\{\\mathcal\{H\}\}, while the LSTM parametrizes the history argument\. These components are described in Sections[2\.4\.1](https://arxiv.org/html/2606.28996#S2.SS4.SSS1)and[2\.4\.2](https://arxiv.org/html/2606.28996#S2.SS4.SSS2), respectively, and combined into the integrated architecture in Section[2\.4\.3](https://arxiv.org/html/2606.28996#S2.SS4.SSS3)\.

#### 2\.4\.1Encoding Fiber Interactions via GNN

Each discrete Voronoi domainΩc\\Omega\_\{c\}was represented as an undirected weighted graph𝒢Ωc=\(𝒱Ωc,ℰΩc,𝒳Ωc\)\\mathcal\{G\}\_\{\\Omega\_\{c\}\}=\(\\mathcal\{V\}\_\{\\Omega\_\{c\}\},\\mathcal\{E\}\_\{\\Omega\_\{c\}\},\\mathcal\{X\}\_\{\\Omega\_\{c\}\}\), where𝒱Ωc\\mathcal\{V\}\_\{\\Omega\_\{c\}\}is the set of nodes representing individual fibers within the cell,ℰΩc\\mathcal\{E\}\_\{\\Omega\_\{c\}\}is the set of weighted edges encoding shear\-lag\-mediated load\-transfer pathways between fiber pairs, and𝒳Ωc∈ℝNf×8\\mathcal\{X\}\_\{\\Omega\_\{c\}\}\\in\\mathbb\{R\}^\{N\_\{f\}\\times 8\}is the node\-feature matrix whose rows store per\-fiber geometric attributes, withNf=\|𝒱Ωc\|N\_\{f\}=\|\\mathcal\{V\}\_\{\\Omega\_\{c\}\}\|denoting the number of fibers in the cell\. The microstructural input to the GNN for each fiber nodennconsists of its three\-dimensional start and end coordinates,𝐱n\(s\),𝐱n\(e\)∈ℝ3\\mathbf\{x\}\_\{n\}^\{\(s\)\},\\mathbf\{x\}\_\{n\}^\{\(e\)\}\\in\\mathbb\{R\}^\{3\}, extracted directly fromμ\\mu\-CT, together with the cell\-level fiber volume fractionVf,ΩcV\_\{f,\\Omega\_\{c\}\}obtained from phase segmentation of the reconstructed volume \(Figure[3](https://arxiv.org/html/2606.28996#S2.F3)\)\. From these primitives, the fiber lengthln=‖𝐱n\(e\)−𝐱n\(s\)‖l\_\{n\}=\\\|\\mathbf\{x\}\_\{n\}^\{\(e\)\}\-\\mathbf\{x\}\_\{n\}^\{\(s\)\}\\\|, orientation vector𝐭n=\(𝐱n\(e\)−𝐱n\(s\)\)/ln\\mathbf\{t\}\_\{n\}=\(\\mathbf\{x\}\_\{n\}^\{\(e\)\}\-\\mathbf\{x\}\_\{n\}^\{\(s\)\}\)/l\_\{n\}, and per\-fiber volume\-fraction contributionVf,n=Vf,Ωc​ln/∑klkV\_\{f,n\}=V\_\{f,\\Omega\_\{c\}\}\\,l\_\{n\}/\\sum\_\{k\}l\_\{k\}were computed for each Voronoi cell\. To accelerate learning, an 8\-dimensional feature vector is defined per fiber:

𝐱n=\[ln,Vf,n,𝐱n\(s\),𝐱n\(e\)\]∈ℝ8\\mathbf\{x\}\_\{n\}=\\left\[l\_\{n\},\\;V\_\{f,n\},\\;\\mathbf\{x\}\_\{n\}^\{\(s\)\},\\;\\mathbf\{x\}\_\{n\}^\{\(e\)\}\\right\]\\in\\mathbb\{R\}^\{8\}\(7\)
Figure[7](https://arxiv.org/html/2606.28996#S2.F7)depicts the graph representation for a discrete Voronoi cell\. Here, each fiber nodenis represented using two descriptors: \(i\) the local fiber volume\-fraction contributionVf,nV\_\{f,n\}and \(ii\) the spatial endpoint coordinates𝐱n\(s\)\\mathbf\{x\}\_\{n\}^\{\(s\)\}and𝐱n\(e\)\\mathbf\{x\}\_\{n\}^\{\(e\)\}, which together encode the fiber length and orientation\. The edges of𝒢Ωc\\mathcal\{G\}\_\{\\Omega\_\{c\}\}are assigned weightswp​qw\_\{pq\}that quantify the load\-transfer efficiency between fiber pairs \(see Figure[7](https://arxiv.org/html/2606.28996#S2.F7)\)\.

A note about edge weights: As previously mentioned, the edge weights capture the load\-transfer efficiency between adjacent fiber pairs via matrix\-mediated shear\. In short\-fiber composites, axial stress is transferred between adjacent fiber ends through interfacial shear stress transfer that occurs over a characteristic shear transfer lengthLcL\_\{c\}\. A neighboring fiber at end\-to\-end separationdp​qd\_\{pq\}therefore experiences a shear stress whose magnitude is proportional toexp⁡\(−dp​q/Lc\)\\exp\(\-d\_\{pq\}/L\_\{c\}\): fiber pairs withinLcL\_\{c\}interact strongly, while pairs separated by distances much greater thanLcL\_\{c\}are effectively mechanically decoupled\.LcL\_\{c\}was computed from Nairn’s \(1997\) shear\-lag formulation\[[51](https://arxiv.org/html/2606.28996#bib.bib59),[39](https://arxiv.org/html/2606.28996#bib.bib58),[15](https://arxiv.org/html/2606.28996#bib.bib60),[43](https://arxiv.org/html/2606.28996#bib.bib55),[44](https://arxiv.org/html/2606.28996#bib.bib56),[46](https://arxiv.org/html/2606.28996#bib.bib57)\], which extended the classical Cox result to account for the geometric confinement of the matrix between neighboring fibers:

Lc=rf​E11f​ln⁡\(R/rf\)2​π​G12mL\_\{c\}=r\_\{f\}\\sqrt\{\\frac\{E^\{f\}\_\{11\}\\,\\ln\(R/r\_\{f\}\)\}\{2\\pi\\,G^\{m\}\_\{12\}\}\}\(8\)
whererfr\_\{f\}is the fiber radius,E11fE^\{f\}\_\{11\}is the fiber axial Young’s modulus \(Table 1\),G12mG^\{m\}\_\{12\}is the matrix in\-plane shear modulus \(Table 1\), andR=rf​π/\(2​3​Vf,Ωc\)R=r\_\{f\}\\sqrt\{\\pi/\(2\\sqrt\{3\}\\,V\_\{f,\\Omega\_\{c\}\}\)\}is the mean matrix cylinder radius determined by local fiber packing geometry within the Voronoi cell\.

The pairwise interaction weight between fibersppandqqwas defined as:

wp​q=exp⁡\(−dp​qLc\)​\(0\.90​\|𝐭p⋅𝐭q\|\+0\.10\)w\_\{pq\}=\\exp\\\!\\left\(\-\\frac\{d\_\{pq\}\}\{L\_\{c\}\}\\right\)\\\!\\left\(0\.90\\,\\lvert\\mathbf\{t\}\_\{p\}\\cdot\\mathbf\{t\}\_\{q\}\\rvert\+0\.10\\right\)\(9\)
wheredp​qd\_\{pq\}was the minimum segment\-to\-segment distance between fibersppandqq, and𝐭p\\mathbf\{t\}\_\{p\},𝐭q\\mathbf\{t\}\_\{q\}are their unit direction vectors\. The exponential term encodes the attenuation of the shear\-mediated stress field with fiber–fiber separation: as fiber ends move apart, the driving stress available to transfer load into a neighboring fiber diminishes exponentially over the lengthLcL\_\{c\}\.

The orientation\-dependent edge weight is defined as\(0\.90​\|𝐭p⋅𝐭q\|\+0\.10\)\(0\.90\\lvert\\mathbf\{t\}\_\{p\}\\cdot\\mathbf\{t\}\_\{q\}\\rvert\+0\.10\), capturing both axial and residual interaction mechanisms\. The dominant term \(0\.90\) represents shear\-lag\-driven axial coupling, which scales with the cosine of the inter\-fiber angle\. It is strongest for parallel fibers, while highly misaligned fibers contribute minimally to axial stress redistribution\. The constant term \(0\.10\) accounts for a residual, orientation\-independent interaction arising from transverse stress concentration due to fiber–matrix elastic stiffness mismatch, ensuring nonzero coupling even for orthogonal fibers\. The0\.90/0\.100\.90/0\.10split was determined from the normalized peak interfacial shear stress ratio obtained in the two\-fiber interaction finite element simulations \(Appendix[Appendix A1\. Parametric Finite\-Element Validation of the Edge Weight Function](https://arxiv.org/html/2606.28996#Sx1)\)\. These simulations shown that the peak interfacial shear stress for orthogonal \(ϕ=90∘\\phi=90^\{\\circ\}\) fibers dropped to approximately 10% \(ranging between 8–22% depending on the inter\-fiber gap distance\) of the fully aligned \(ϕ=0∘\\phi=0^\{\\circ\}\) baseline configuration\. This baseline residual ratio served as the basis for using 0\.10, while the remaining normalized load\-bearing capacity \(0\.90\) was allocated to the orientation\-dependent cosine term\.

The edge weight function operates strictly at the pairwise level\. Higher\-order multi\-fiber interaction effects—such as load redistribution within fiber clusters, and cooperative load sharing—are not explicitly encoded in the edge weights but emerge through neighborhood aggregation in the GNN message\-passing operation \(Equation \([10](https://arxiv.org/html/2606.28996#S2.E10)\)\)\. The graph is undirectedwp​q=wq​pw\_\{pq\}=w\_\{qp\}and complete: every pair of fibers within a Voronoi cell is connected\. A lower bound of10−610^\{\-6\}was imposed on edge weights to maintain numerical stability of the degree\-normalized aggregation in Equation \([10](https://arxiv.org/html/2606.28996#S2.E10)\), without altering the graph topology\. The parametric finite\-element study \(Appendix[Appendix A1\. Parametric Finite\-Element Validation of the Edge Weight Function](https://arxiv.org/html/2606.28996#Sx1)\) further confirmed that both the cosine\-dependent axial coupling and the orientation\-independent residual are physically operative for fiber geometries and spacings consistent withμ\\mu\-CT observations\.

![Refer to caption](https://arxiv.org/html/2606.28996v1/Figure7.png)Figure 7:Graph construction and GNN encoding for a Voronoi cell\.Left: Fibers are represented as graph nodes, with node features given by local volume fractionVf​i​b​e​r,ΩcV\_\{fiber,\\Omega\_\{c\}\}and spatial endpoints𝐱n\(s\)\\mathbf\{x\}\_\{n\}^\{\(s\)\},𝐱n\(e\)\\mathbf\{x\}\_\{n\}^\{\(e\)\}\. Fiber–fiber interactions were modeled as weighted edges\(wp​q\)\(w\_\{pq\}\), computed from geometric proximity and directional alignment using a shear\-lag\-informed weight function \(Equation[9](https://arxiv.org/html/2606.28996#S2.E9)\)\.Center: The graph is processed through two stacked graph convolutional layers, where message passing and neighborhood aggregation update node embeddings to reflect local load redistribution\.Right: Masked global mean pooling aggregates the final node embeddings into a compact 16\-dimensional region\-level microstructural descriptorzΩc∈ℝ16z\_\{\\Omega\_\{c\}\}\\in\\mathbb\{R\}^\{16\}, representing the combined effect of fiber geometry and interaction topology on the Voronoi cell’s effective response\.Graph convolutional network \(GCN\) layer formulation: Node features were first projected from the 8\-dimensional input space to a 64\-dimensional embedding space via a linear layer with GELU activation\. Two stacked graph convolutional layers then updated node embeddings through degree\-normalized weighted aggregation with residual connections:

𝐡n\(l\+1\)=𝐡n\(l\)\+Dropout​\(LN​\(GELU​\(∑m∈𝒩​\(n\)wn​md^n​d^m​𝐖\(l\)​𝐡m\(l\)\)\)\)\\mathbf\{h\}\_\{n\}^\{\(l\+1\)\}=\\mathbf\{h\}\_\{n\}^\{\(l\)\}\+\\mathrm\{Dropout\}\\\!\\left\(\\mathrm\{LN\}\\\!\\left\(\\mathrm\{GELU\}\\\!\\left\(\\sum\_\{m\\in\\mathcal\{N\}\(n\)\}\\frac\{w\_\{nm\}\}\{\\sqrt\{\\hat\{d\}\_\{n\}\\,\\hat\{d\}\_\{m\}\}\}\\,\\mathbf\{W\}^\{\(l\)\}\\,\\mathbf\{h\}\_\{m\}^\{\(l\)\}\\right\)\\right\)\\right\)\(10\)
where𝐡m\(l\)\\mathbf\{h\}\_\{m\}^\{\(l\)\}is the embedding of neighboring fiber nodemmat layerll,wn​mw\_\{nm\}is the edge weight from Equation[9](https://arxiv.org/html/2606.28996#S2.E9),d^n=1\+∑jwn​j\\hat\{d\}\_\{n\}=1\+\\sum\_\{j\}w\_\{nj\}is the weighted node degree with self\-loop,𝐖\(l\)\\mathbf\{W\}^\{\(l\)\}is the trainable weight matrix at layerll, LN denotes layer normalization, and the residual connection𝐡n\(l\)\\mathbf\{h\}\_\{n\}^\{\(l\)\}stabilized gradient flow across layers\. The degree normalization factor1/d^n​d^m1/\\sqrt\{\\hat\{d\}\_\{n\}\\hat\{d\}\_\{m\}\}prevents embedding magnitudes from scaling with local fiber connectivity, ensuring that fibers in dense clusters were not artificially amplified relative to isolated fibers\. From a mechanical standpoint, this message\-passing operation approximated multi\-fiber load redistribution within the fiber network: edge weights encoded pairwise shear\-mediated stress\-transfer efficiency, neighborhood aggregation represented multi\-fiber load sharing among all neighboring fibers, and stacked layers iteratively propagated this redistribution across the interaction network—analogous to FE solver iterations\.

After two message\-passing layers, masked global mean pooling aggregated active node embeddings into a 64\-dimensional graph\-level representation, which was projected by a linear layer to a compact 16\-dimensional microstructural descriptor:

zΩc=𝐖pool​\(1Nf​∑n=1Nf𝐡n\(2\)\)∈ℝ16z\_\{\\Omega\_\{c\}\}=\\mathbf\{W\}\_\{\\mathrm\{pool\}\}\\\!\\left\(\\frac\{1\}\{N\_\{f\}\}\\sum\_\{n=1\}^\{N\_\{f\}\}\\mathbf\{h\}\_\{n\}^\{\(2\)\}\\right\)\\in\\mathbb\{R\}^\{16\}\(11\)
This embedding served as a compact, fixed\-length representation of the cell’s effective stiffness potential and interaction topology, capturing the combined influence of fiber geometry, spatial distribution, volume fraction, and interaction network structure on the cell’s effective mechanical response\. BecausezΩcz\_\{\\Omega\_\{c\}\}was derived exclusively from quantities computable fromμ\\mu\-CT fiber coordinates andVf,ΩcV\_\{f,\\Omega\_\{c\}\}, the GNN required no microstructural inputs beyond standard reconstruction outputs\.

#### 2\.4\.2Modeling History\-Dependent Nonlinear Evolution via LSTM

The homogenized stress response of each cell was a history\-dependent trajectory that reflected progressive accumulation of matrix plastic strainepe\_\{p\}and scalar damageDmD\_\{m\}under monotonic loading\. Because these internal variables evolved with loading history, the stress at any increment depended not only on the current strain but on the entire prior deformation path\. A feedforward network that treats each strain increment independently could not enforce this temporal consistency\. A LSTM network was therefore employed to process the loading history as a sequence, maintaining hidden state vectors𝐡t∈ℝ128\\mathbf\{h\}\_\{t\}\\in\\mathbb\{R\}^\{128\}and cell state vectors𝐜t∈ℝ128\\mathbf\{c\}\_\{t\}\\in\\mathbb\{R\}^\{128\}as surrogate internal state variables implicitly encoding accumulated plasticity and damage without explicitly tracking them\.

LSTM input features: The LSTM input at each strain incrementttwas a 21\-dimensional feature vector concatenated from two categories of inputs:

𝐪t=\[εtεpk,\(εtεpk\)2,\(εtεpk\)3,ln⁡\(1\+εtεpk\),Vf,Ωc\]⏟5 loading\-side features∥zΩc⏟16 D topology\\mathbf\{q\}\_\{t\}=\\underbrace\{\\left\[\\frac\{\\varepsilon\_\{t\}\}\{\\varepsilon\_\{\\mathrm\{pk\}\}\},\\;\\left\(\\frac\{\\varepsilon\_\{t\}\}\{\\varepsilon\_\{\\mathrm\{pk\}\}\}\\right\)^\{\\\!2\},\\;\\left\(\\frac\{\\varepsilon\_\{t\}\}\{\\varepsilon\_\{\\mathrm\{pk\}\}\}\\right\)^\{\\\!3\},\\;\\ln\\\!\\left\(1\+\\frac\{\\varepsilon\_\{t\}\}\{\\varepsilon\_\{\\mathrm\{pk\}\}\}\\right\),\\;V\_\{f,\\Omega\_\{c\}\}\\right\]\}\_\{\\text\{5 loading\-side features\}\}\\;\\Big\\\|\\;\\underbrace\{z\_\{\\Omega\_\{c\}\}\}\_\{\\text\{16 D topology\}\}\(12\)These input vector maintains a fixed size of 21 dimensions across all strain increments rather than growing with the sequence length\. The history of prior deformation is captured entirely through the step\-by\-step updates of the fixed\-dimensional internal state\. The first category comprised five loading\-side features derived analytically from the applied scalar strainεt\\varepsilon\_\{t\}and the fiber coordinate inputs used by the GNN\. The polynomial and logarithmic strain features, motivated by the hyperelastic model, enriched the nonlinear strain representation, enabling the LSTM to distinguish between the linear elastic, nonlinear hardening, and near\-peak softening regimes without requiring an explicit constitutive form\. Normalizing strain byεpk\\varepsilon\_\{\\mathrm\{pk\}\}mapped the stress peak of every cell to a universal reference strain of 1\.0, enabling the LSTM to learn a shape that was subsequently amplitude\-rescaled\. The cell\-level volume fractionVf,ΩcV\_\{f,\\Omega\_\{c\}\}directly conveyed the available reinforcement content\.

The LSTM hidden and cell state vectors evolved as:

\(𝐡t,𝐜t\)=LSTM​\(𝐪t,𝐡t−1,𝐜t−1\)\\left\(\\mathbf\{h\}\_\{t\},\\,\\mathbf\{c\}\_\{t\}\\right\)=\\mathrm\{LSTM\}\\\!\\left\(\\mathbf\{q\}\_\{t\},\\;\\mathbf\{h\}\_\{t\-1\},\\;\\mathbf\{c\}\_\{t\-1\}\\right\)\(13\)
with two stacked layers \(hidden dimension 128, inter\-layer dropout 0\.25\)\. The hidden and cell states serve as surrogate internal variables implicitly encoding accumulated plastic strain and damage history at each loading increment\. The GNN component performed geometric and topological encoding of the fiber interaction network, while the LSTM modeled the path\-dependent evolution of the constitutive response — a separation that mirrors classical constitutive modeling, in which microstructural descriptors govern stiffness capacity while evolution laws govern nonlinear progression\.

#### 2\.4\.3Integrated GNN–LSTM Architecture

The final surrogate mapping from fiber interaction topology to normalized stress trajectory was:

σ^t\(Ωc\)=Softplus​\(ℱout​\(\[𝐡t;zΩc\]\)\),𝐡t=ℱLSTM​\(𝐪t,𝐡t−1,𝐜t−1\),zΩc=ℱGNN​\(𝒢Ωc\)\\hat\{\\sigma\}\_\{t\}^\{\(\\Omega\_\{c\}\)\}=\\mathrm\{Softplus\}\\\!\\left\(\\mathcal\{F\}\_\{\\mathrm\{out\}\}\\\!\\left\(\\left\[\\mathbf\{h\}\_\{t\}\\;;\\;z\_\{\\Omega\_\{c\}\}\\right\]\\right\)\\right\),\\quad\\mathbf\{h\}\_\{t\}=\\mathcal\{F\}\_\{\\mathrm\{LSTM\}\}\\\!\\left\(\\mathbf\{q\}\_\{t\},\\,\\mathbf\{h\}\_\{t\-1\},\\,\\mathbf\{c\}\_\{t\-1\}\\right\),\\quad z\_\{\\Omega\_\{c\}\}=\\mathcal\{F\}\_\{\\mathrm\{GNN\}\}\\\!\\left\(\\mathcal\{G\}\_\{\\Omega\_\{c\}\}\\right\)\(14\)
whereℱGNN​\(𝒢Ωc\)\\mathcal\{F\}\_\{\\mathrm\{GNN\}\}\(\\mathcal\{G\}\_\{\\Omega\_\{c\}\}\)denoted the GNN encoder operating on the fiber interaction graph𝒢Ωc\\mathcal\{G\}\_\{\\Omega\_\{c\}\}of cellΩc\\Omega\_\{c\}, andℱout\\mathcal\{F\}\_\{\\mathrm\{out\}\}was a three\-layer residual output network\. Figure[8](https://arxiv.org/html/2606.28996#S2.F8)illustrates the complete architecture\.

The output network received a skip\-connected input formed by concatenating the LSTM hidden state𝐡t∈ℝ128\\mathbf\{h\}\_\{t\}\\in\\mathbb\{R\}^\{128\}with the GNN embeddingzΩc∈ℝ16z\_\{\\Omega\_\{c\}\}\\in\\mathbb\{R\}^\{16\}at every strain increment, producing a 144\-dimensional input to a three\-layer MLP \(144→64→32→1144\\to 64\\to 32\\to 1, GELU activations\)\. The skip connection ensured that the GNN\-encoded topology signal—which carried direct information about the cell’s peak stress potential through its fiber arrangement—remained accessible at every prediction step independently of LSTM hidden state depth\. Without this connection, the GNN amplitude signal would need to propagate through all LSTM gates across potentially 100 or more time steps, risking attenuation through gate saturation\. A Softplus nonlinearity, defined asf​\(x\)=ln⁡\(1\+ex\)f\(x\)=\\ln\(1\+e^\{x\}\), at the network output enforced the physical constraint that predicted stress remained strictly non\-negative throughout the trajectory\.

To recover stress in physical units, the surrogate output was trained against FE\-computed stress trajectories normalized by the cell\-level peak stress, with the peak stress predicted independently from the GNN embedding via a post\-hoc regression model trained on the same FE dataset\. This architecture enforced a physically meaningful factorization: the GNN captured topology\-dependent stiffness potential from the fiber interaction structure; the LSTM captured the path\-dependent evolution of the constitutive response\.

The peak stress amplitude was predicted using a Random Forest regressor \(300 trees, MultiOutputRegressor\) trained post\-hoc on the frozen GNN embeddingzΩcz\_\{\\Omega\_\{c\}\}concatenated with the cell\-level volume fractionVf,ΩcV\_\{f,\\Omega\_\{c\}\}, using the same 500\-cell dataset and 70–15–15 train–validation–test split\. Although the GNN auxiliary head also predicted peak stress during joint training, it exhibited generalization degradation on unseen cells due to its coupling with the evolving LSTM gradients during end\-to\-end optimization\. The post\-hoc Random Forest, applied to frozen embeddings after GNN–LSTM convergence, provided more robust amplitude generalization without requiring additional finite\-element simulations, achieving a coefficient of determinationR2=0\.98R^\{2\}=0\.98on the held\-out test set for peak stress prediction compared with the discrete FE model\.

![Refer to caption](https://arxiv.org/html/2606.28996v1/Figure8.png)Figure 8:Integrated GNN–LSTM architecture for physics\-informed, history\-dependent stress–strain prediction\. The GNN microstructure encoder \(left\) processes the fiber interaction graph𝒢Ωc\\mathcal\{G\}\_\{\\Omega\_\{c\}\}to produce a fixed 16\-dimensional cell embeddingzΩc∈ℝ16z\_\{\\Omega\_\{c\}\}\\in\\mathbb\{R\}^\{16\}\. At each strain incrementtt, this embedding is concatenated with four nonlinear strain features and the cell\-level volume fractionVf,ΩcV\_\{f,\\Omega\_\{c\}\}to form the 21\-dimensional input vector𝐪t∈ℝ21\\mathbf\{q\}\_\{t\}\\in\\mathbb\{R\}^\{21\}, which is fed sequentially into a two\-layer LSTM \(hidden dimension 128\)\. The LSTM hidden state𝐡t∈ℝ128\\mathbf\{h\}\_\{t\}\\in\\mathbb\{R\}^\{128\}evolves as a surrogate set of internal state variables encoding accumulated plasticity and damage history\. A skip connection concatenateszΩcz\_\{\\Omega\_\{c\}\}directly with𝐡t\\mathbf\{h\}\_\{t\}, producing a 144\-dimensional input to a three\-layer output MLP \(144→64→32→1144\\rightarrow 64\\rightarrow 32\\rightarrow 1, Softplus\), which yields the normalized predicted stressσ^t\(Ωc\)\\hat\{\\sigma\}\_\{t\}^\{\(\\Omega\_\{c\}\)\}at each increment\.
#### 2\.4\.4Physics\-Informed Regularization

Composite loss function: To ensure physically plausible surrogate predictions, a composite loss function was used alongside additional structural training constraints\. Model parameters were optimized by minimizing:

ℒ=ℒHuber\+0\.05​ℒmono\+0\.001​ℒsmooth\+0\.02​ℒzero\+0\.3​ℒaux\+0\.5​ℒpeak\\mathcal\{L\}=\\mathcal\{L\}\_\{\\mathrm\{Huber\}\}\+0\.05\\,\\mathcal\{L\}\_\{\\mathrm\{mono\}\}\+0\.001\\,\\mathcal\{L\}\_\{\\mathrm\{smooth\}\}\+0\.02\\,\\mathcal\{L\}\_\{\\mathrm\{zero\}\}\+0\.3\\,\\mathcal\{L\}\_\{\\mathrm\{aux\}\}\+0\.5\\,\\mathcal\{L\}\_\{\\mathrm\{peak\}\}\(15\)The scalar scaling coefficients were assigned to each loss component \(0\.05, 0\.001, 0\.02, 0\.3, and 0\.5\) and function as regularization hyperparameters\. These constants were determined empirically via a grid search on the validation subset to ensure stable convergence during joint training, balancing trajectory\-tracking accuracy with physical admissibility constraints \(e\.g\., preventing non\-physical oscillations without over\-penalizing post\-yield softening\)\.The primary termℒHuber\\mathcal\{L\}\_\{\\mathrm\{Huber\}\}was a strain\-weighted Huber loss \(δ=0\.1\\delta=0\.1\) computed exclusively over real,non\-padded data points:

ℒHuber=∑twt​ℋδ​\(σ^t−σtFE\)∑twt\\mathcal\{L\}\_\{\\mathrm\{Huber\}\}=\\frac\{\\displaystyle\\sum\_\{t\}w\_\{t\}\\,\\mathcal\{H\}\_\{\\delta\}\\\!\\left\(\\hat\{\\sigma\}\_\{t\}\-\\sigma\_\{t\}^\{\\mathrm\{FE\}\}\\right\)\}\{\\displaystyle\\sum\_\{t\}w\_\{t\}\}\(16\)
whereσ^t\\hat\{\\sigma\}\_\{t\}denotes the surrogate\-predicted normalized stress at strain incrementtt,σtFE\\sigma\_\{t\}^\{\\mathrm\{FE\}\}represents the corresponding high\-fidelity FE homogenization target stress, andwtw\_\{t\}is the strain\-dependent tracking weight matrix\. The tracking weights are defined explicitly as:

wt=1\+3​\[clip​\(εt/εpk−0\.80\.8,0,1\)\]2w\_\{t\}=1\+3\\\!\\left\[\\mathrm\{clip\}\\\!\\left\(\\frac\{\\varepsilon\_\{t\}/\\varepsilon\_\{\\mathrm\{pk\}\}\-0\.8\}\{0\.8\},\\;0,\\;1\\right\)\\right\]^\{\\\!2\}\(17\)whereεt\\varepsilon\_\{t\}is the applied macroscopic strain at steptt, andεpk\\varepsilon\_\{\\mathrm\{pk\}\}is the localized cell\-level strain corresponding to the peak stress point\. The termclip​\(x,a,b\)\\mathrm\{clip\}\(x,a,b\)represents a standard mathematical bounding operation defined to constrain the scaling window strictly betweenaaandbb\.

This custom weighting scheme is engineered to prioritize tracking accuracy along highly nonlinear paths\. It assigns a tracking penalty up to4×4\\timeshigher to data points approaching and exceeding the stress peak \(εt/εpk≥0\.8\\varepsilon\_\{t\}/\\varepsilon\_\{\\mathrm\{pk\}\}\\geq 0\.8\), where accurate prediction of the onset and progression of nonlinear matrix softening is most mechanically significant\. The Huber loss componentℋδ\\mathcal\{H\}\_\{\\delta\}is intentionally employed in preference to a standard mean\-squared error framework to reduce gradient sensitivity to extreme outlier stress fluctuations arising from highly localized damage events in individual cells, while retaining stable quadratic penalization for minor prediction errors\.

Pre\-peak monotonicity: Stress increments prior to the peak were regularized to prevent non\-physical oscillations:

ℒmono=1tpk−1​∑t=1tpk−1GELU​\(σ^t−1−σ^t\)\\mathcal\{L\}\_\{\\mathrm\{mono\}\}=\\frac\{1\}\{t\_\{\\mathrm\{pk\}\}\-1\}\\sum\_\{t=1\}^\{t\_\{\\mathrm\{pk\}\}\-1\}\\mathrm\{GELU\}\\\!\\left\(\\hat\{\\sigma\}\_\{t\-1\}\-\\hat\{\\sigma\}\_\{t\}\\right\)\(18\)
wheretpkt\_\{\\mathrm\{pk\}\}was the index of the true stress peak identified from the FE data for each cell\. The post\-peak regime was explicitly excluded from this constraint, consistent with the experimentally observed and FE\-computed post\-peak softening behavior of the ABS matrix following peak stress\.

Smoothness:Second\-order differences of the predicted stress trajectory were penalized to suppress artificial stiffness fluctuations:

ℒsmooth=1T−2​∑t=1T−2\(σ^t\+1−2​σ^t\+σ^t−1\)2\\mathcal\{L\}\_\{\\mathrm\{smooth\}\}=\\frac\{1\}\{T\-2\}\\sum\_\{t=1\}^\{T\-2\}\\left\(\\hat\{\\sigma\}\_\{t\+1\}\-2\\hat\{\\sigma\}\_\{t\}\+\\hat\{\\sigma\}\_\{t\-1\}\\right\)^\{\\\!2\}\(19\)whereTTrepresents the total sequence length \(total number of strain increments\) in the complete loading history\.\.

Elastic consistency:The predicted normalized stress at the first strain increment was penalized to enforce physically correct onset from zero stress:

ℒzero=σ^02\\mathcal\{L\}\_\{\\mathrm\{zero\}\}=\\hat\{\\sigma\}\_\{0\}^\{\\,2\}\(20\)whereσ^0\\hat\{\\sigma\}\_\{0\}is the surrogate\-predicted normalized stress at the initial strain increment \(t=0t=0\), ensuring that the predicted stress trajectory originates strictly from a state of zero stress\.

Auxiliary physics supervision:The auxiliary lossℒaux\\mathcal\{L\}\_\{\\mathrm\{aux\}\}penalized deviation of the GNN auxiliary head predictions\[E~,σ~pk,ε~pk,δ~\]\[\\tilde\{E\},\\tilde\{\\sigma\}\_\{\\mathrm\{pk\}\},\\tilde\{\\varepsilon\}\_\{\\mathrm\{pk\}\},\\tilde\{\\delta\}\]from their normalized FE\-computed targets, ensuring that the 16\-dimensional cell embedding encoded mechanically interpretable information—initial modulus, peak stress, peak strain, and post\-peak softening rate—throughout training rather than converging to a representation that satisfied the trajectory loss without physical meaning\.

Peak stress penalty: The term:

ℒpeak=‖σ^tpk−σtpkFE‖2\\mathcal\{L\}\_\{\\text\{peak\}\}=\\left\\\|\\hat\{\\mathbf\{\\sigma\}\}\_\{t\_\{\\mathrm\{pk\}\}\}\-\\mathbf\{\\sigma\}^\{\\text\{FE\}\}\_\{t\_\{\\mathrm\{pk\}\}\}\\right\\\|^\{2\}\(21\)applied concentrated penalization at the exact peak stress point, improving prediction accuracy at the mechanically critical transition between hardening and softening regimes independently of the weighted Huber loss\.

Optimizer and training schedule: The GNN and LSTM were trained jointly end\-to\-end using the AdamW optimizer \(initial learning rate5×10−45\\times 10^\{\-4\}, weight decay2×10−42\\times 10^\{\-4\}\)\[[42](https://arxiv.org/html/2606.28996#bib.bib73)\]with gradient norm clipping to 1\.0 to prevent unrealistically high gradients during early training epochs\. A ReduceLROnPlateau learning rate scheduler\[[8](https://arxiv.org/html/2606.28996#bib.bib74)\]halved the learning rate when the validation loss plateaued for 15 consecutive epochs, with a minimum learning rate of10−610^\{\-6\}\. Training ran for up to 800 epochs; early stopping was applied based on validation loss, with a patience of 57 epochs, to prevent overfitting\. The complete set of training hyperparameters is summarized in Appendix[Appendix A2\. GNN\-LSTM Training Configuration](https://arxiv.org/html/2606.28996#Sx2)\.

#### 2\.4\.5Error Propagation and Stability

Small\-cell\-level prediction errors with the GNN\-LSTM surrogate did not lead to any observed instabilities in the coupon\-scale simulations\. Because the coupon\-level response was obtained by volume averaging over many discrete cells, surrogate approximation errors were statistically mitigated through the averaging operation itself\. Robust macroscopic predictions were therefore achievable even in the presence of small cell\-level surrogate approximation errors, with the principal sources of coupon\-scale discrepancy attributable to structural idealizations in the Voronoi representation rather than to surrogate approximation errors\.

### 2\.5Mesoscale Assembly and Coupon\-Scale Prediction

Having obtained the effective stress–strain response of each discrete Voronoi cell—either through direct FE homogenization \(Section[2\.3\.3](https://arxiv.org/html/2606.28996#S2.SS3.SSS3)\) or through the GNN\-LSTM surrogate \(Section[2\.4](https://arxiv.org/html/2606.28996#S2.SS4)\)—the collection of cell responses defines a heterogeneous mesoscale constitutive field over the reconstructed volume\. The coupon\-scale response was obtained by assembling the Voronoi cell responses through volume\-weighted averaging, consistent with the Hill–Mandel homogenization principle:

𝚺​\(ϵ\)=1Vcoupon​∑ΨcVΨc​σ¯\(Ψc\)​\(ϵ\)\\mathbf\{\\Sigma\}\(\\mathbf\{\\epsilon\}\)=\\frac\{1\}\{V\_\{\\text\{coupon\}\}\}\\sum\_\{\\Psi\_\{c\}\}V\_\{\\Psi\_\{c\}\}\\bar\{\\mathbf\{\\sigma\}\}^\{\(\\Psi\_\{c\}\)\}\(\\mathbf\{\\epsilon\}\)\(22\)where𝚺\\mathbf\{\\Sigma\}is the macroscopic second\-order stress tensor of the coupon,VcouponV\_\{\\text\{coupon\}\}is the total volume of the coupon,Ψc\\Psi\_\{c\}denotes thecc\-th mesoscale Voronoi cell within the reconstructed assembly,VΨcV\_\{\\Psi\_\{c\}\}is the volume of that corresponding cell, andσ¯\(Ψc\)\\bar\{\\mathbf\{\\sigma\}\}^\{\(\\Psi\_\{c\}\)\}is the homogenized second\-order stress tensor predicted by the trained GNN\-LSTM surrogate predictions under the applied macroscopic strain tensorϵ\\epsilon\.

Because the macroscopic response is obtained by aggregation over many cells, moderate cell\-level surrogate errors are naturally stabilized through averaging\. This enables robust coupon\-scale prediction while preserving sensitivity to the spatial distribution of fiber interaction topology\.

## 3Results and Discussion

### 3\.1Microstructural Evolution under Cyclic Loading

μ\\mu\-CT analysis revealed that pore nucleation and growth were highly localized\. Cyclic loading induced measurable shifts in ensemble\-level microstructural descriptors, with pores forming preferentially in \(i\) fiber\-end neighborhoods and \(ii\) clustered fiber regions where local matrix ligaments are thin \(Figure[9](https://arxiv.org/html/2606.28996#S3.F9)a\)\. These locations coincide with regions of elevated shear and normal stress concentration as predicted by shear\-lag theory and inter\-fiber load\-transfer mechanics\. The preferential localization of porosity near fiber termini provided direct experimental evidence that stress\-transfer discontinuities at fiber ends serve as primary damage\-initiation sites\.

The nearest\-neighbor probability densityHp​\(r\)H\_\{p\}\(r\)provides complementary quantitative evidence of progressive fiber\-cluster dispersion accompanying pore growth \(Figure[9](https://arxiv.org/html/2606.28996#S3.F9)b\)\. Before testing,Hp​\(r\)H\_\{p\}\(r\)exhibited a single dominant peak nearr≈20r\\approx 20–22​μ​m22~\\mu\\text\{m\}, consistent with tight fiber clustering inherited from the manufacturing process\. Above the fatigue limit, however, a secondary peak emerged nearr≈45r\\approx 45–47​μ​m47~\\mu\\text\{m\}, producing a distinct bimodal distribution that signifies a clear signature of cluster dilation\.

While fiber rotation is geometrically constrained within a pristine, continuous polymer matrix, this observable shift in the FOD is a direct structural consequence of severe, localized matrix failure that accumulates during cyclic loading\. As fatigue cycles progress, a massive threefold increase in localized matrix porosity occurred preferentially within the thin matrix ligaments separating dense fiber clusters and around stress\-concentrating fiber ends \(see Table[3](https://arxiv.org/html/2606.28996#S3.T3)\)\. This progressive micro\-voiding, debonding, and matrix degradation effectively dismantled the local geometric constraints holding the fibers rigid\. Under macroscopic tensile strain, the fibers adjacent to these expanding void networks are physically tilted and reallocated into the newly vacated micro\-volumes\. This localized rearrangement and shifting of the fiber network under cyclic loading directly aligns with experimental observations of damage\-induced microstructural evolution reported in the literature\[[6](https://arxiv.org/html/2606.28996#bib.bib44),[3](https://arxiv.org/html/2606.28996#bib.bib65)\]\.

![Refer to caption](https://arxiv.org/html/2606.28996v1/Figure9.png)Figure 9:Fatigue\-driven microstructural damage in 20 wt\.% CF\-ABS, fabricated via AM\-CM at ORNL, was characterized usingμ\\mu\-CT\. \(a\) Representative cross\-sectional images at two damage states are shown: prior to fatigue loading \(Before test\) and after exceeding the fatigue limit \(Above\-fatigue limit\)\. Post\-fatigue scans reveal pronounced localized pore nucleation and growth, predominantly near fiber termini and within fiber\-clustering regions, indicating stress\-concentration\-mediated damage initiation and propagation\. \(b\) Nearest\-neighbor probability densityHp​\(r\)H\_\{p\}\(r\)tracked longitudinally for a representative sample \(Specimen 2\) for the Before test \(BT\) and Above fatigue limit \(AF\) stages, indicating progressive cluster dilation driven by pore growth within the inter\-fiber zones, resulting in a dispersed sub\-population of fibers with substantially increased nearest\-neighbor distances\.Table 3:Porosity volume fraction \(%\) evolution under cyclic loading\. Values represent the mean ± standard deviation \(adopted from\[[55](https://arxiv.org/html/2606.28996#bib.bib64)\]\)\.The extracted descriptors revealed a gradual yet systematic microstructural evolution under fatigue loading\. Figure[10](https://arxiv.org/html/2606.28996#S3.F10)shows shifts in orientation distributions and aspect\-ratio statistics, indicative of microstructural rearrangement and spatial reallocation of fiber\-interaction neighborhoods under cyclic loading\. Consistent with these observations, the porosity volume fraction increased progressively under cyclic loading for all three specimens \(Table[3](https://arxiv.org/html/2606.28996#S3.T3)\)\. The presence of pores reduces the effective load\-bearing capacity of the matrix and modifies the local stress state within the composite\. The increase is particularly pronounced for specimens loaded above the fatigue limit, where porosity growth exceeded the initial \(pre\-fatigue\) state by threefold\. Although these changes were moderate at the ensemble level, their spatial localization significantly amplifies the progression of localized damage\.

![Refer to caption](https://arxiv.org/html/2606.28996v1/Figure10.png)Figure 10:Statistical characterization of fiber morphology and its correlation with mechanical response across fatigue stages\.Left: Histograms of measured fiber aspect ratios at each fatigue stage \(BT, BF, AF\), revealing quantifiable shifts in fiber length distribution under cyclic loading\.Center: Histograms of the in\-plane fiber orientation angleϕ\\phimeasured from the x\-axis and angleθ\\thetain the y\-z plane\.Right: Experimentally measured axial stress\-strain responses in fiber\-dominated direction: before fatigue testing \(BT\), below\-fatigue limit \(BF\), and above fatigue limit \(AF\)\. Rows correspond to Specimen 1 \(top\), Specimen 2 \(middle\), and Specimen 3 \(bottom\), which span three distinct initial porosity levels\.When correlated with the macroscale mechanical response \(Figure[10](https://arxiv.org/html/2606.28996#S3.F10)\), the progressive increase in porosity and localized damage accumulation corresponded to only marginal reductions in the longitudinal elastic stiffness and a subtle, early onset of axial nonlinear deformation\. This highly stable axial behavior is consistent with the fiber\-dominated nature of short\-fiber architectures, where the longitudinal modulus is governed primarily by the high\-stiffness reinforcement phase and remains relatively insensitive to low\-to\-moderate variations in matrix void content within the1%1\\%–6%6\\%range\[[57](https://arxiv.org/html/2606.28996#bib.bib75)\]\. However, as fatigue\-induced porosity grows, it severely disrupts matrix continuity and degrades matrix\-mediated shear\-lag load transfer between adjacent fibers\. Consequently, while the axial modulus experiences minor degradation, this microstructural breakdown exerts a far more pronounced influence on the matrix\- and interface\-dominated mechanical properties, specifically the transverse normal and shear responses where voids act as localized stress concentrators\[[41](https://arxiv.org/html/2606.28996#bib.bib76),[51](https://arxiv.org/html/2606.28996#bib.bib59)\]\. Importantly, this distinct sensitivity establishes that the severity of anisotropic stiffness loss is dictated by the local topology of the fiber network and inter\-fiber matrix ligaments rather than a uniform volumetric loss of load\-bearing material\[[1](https://arxiv.org/html/2606.28996#bib.bib61),[55](https://arxiv.org/html/2606.28996#bib.bib64)\]\.

Taken together, these experimental observations establish a clear, path\-dependent microstructure–property relationship\. The localized geometric heterogeneity within the material defines the initial fiber interaction topology; cyclic loading systematically modifies this topology through localized porosity growth and microstructural rearrangement, and the evolved topology ultimately governs the subsequent quasi\-static mechanical response\. These experimentally characterized features form the direct physical basis for the mesoscale interaction\-based modeling framework utilized in the following section\.

### 3\.2Microstructure\-Dependent Stiffness from Finite Element Homogenization

Figure[11](https://arxiv.org/html/2606.28996#S3.F11)presents the local stress distributions for three discrete Voronoi cells subjected to identical far\-field loading conditions\. As shown in Figure[11](https://arxiv.org/html/2606.28996#S3.F11)a–c, these Voronoi cells were selected to represent contrasting microstructural topologies, varying significantly in their local fiber volume fractions, spatial alignments, clustering, and inter\-fiber spacing\. When subjected to simulated loading, these distinct microstructural descriptors produced markedly different mechanical responses, as reflected in the localized stress\-field heterogeneity across panels Figure[11](https://arxiv.org/html/2606.28996#S3.F11)d–f and quantified by the macroscopic stress–strain trajectories in Figure[12](https://arxiv.org/html/2606.28996#S3.F12)\.

![Refer to caption](https://arxiv.org/html/2606.28996v1/Figure11.png)Figure 11:Discrete\-cell multiscale modeling framework\. \(a\) Mesoscale Voronoi tessellation of the reconstructed AM\-CM microstructure\. \(b\) Three representative polyhedral discrete cells \(ID 1, ID 2, ID 3\) depicting different fiber topologies\. \(c\) Equivalent RUC FE models with explicitly resolved fiber geometries embedded in a homogenized porous matrix\. \(d–f\) Representative local stress fields under: fiber\-aligned normal loading in the x\-direction \(left\), transverse loading in the z\-direction \(center\), and axial shear in the x\-y plane \(right\)\. Color scales are independently normalized per loading case to highlight intra\-cell stress distribution patterns\.Under axial loading along the x\-axis \(σx​x\\sigma\_\{xx\}\), cells containing clusters of nearly aligned fibers \(such as ID 1\) exhibited pronounced stress localization along the fiber axes\. This pattern indicates highly efficient load transfer via axial fiber stretching, yielding a correspondingly high effective longitudinal stiffness\. Conversely, cells characterized by dispersed or obliquely oriented fibers \(such as ID 3\) exhibited more diffuse stress fields, in which a significantly larger fraction of the applied far\-field load was carried by the surrounding thermoplastic matrix rather than the reinforcement phase\. Because the continuous polymer matrix is substantially more compliant than the carbon fibers, these matrix\-dominated cells exhibit an earlier onset of localized plastic deformation\. The resulting variation in initial axial stiffness among the three representative cells examined here reaches approximately 8%–12%, highlighting the exceptional sensitivity of macroscale load\-transfer efficiency to localized fiber topology\.

Under in\-plane shear loading \(τx​y\\tau\_\{xy\}\), the stress concentrations arose predominantly within the thin matrix ligaments separating adjacent fibers rather than within the fibers themselves\. This shifts the structural burden, indicating that shear transfer across the microstructural domain was controlled primarily by matrix\-mediated interactions rather than direct fiber reinforcement\. Consequently, the axial shear stiffness was more strongly governed by local inter\-fiber spacing and matrix continuity than by absolute fiber orientation—a mechanistic distinction with critical implications for the anisotropic damage evolution observed at later, highly cycled fatigue stages\.

The cell\-level stress–strain responses shown in Figure[12](https://arxiv.org/html/2606.28996#S3.F12)a–c further illustrate the influence of local fiber topology on the resulting nonlinear constitutive behavior\. Cells whose fibers were preferentially aligned with the primary loading axis exhibited both a higher initial elastic modulus and a more gradual transition to nonlinear softening\. This observation is consistent with classical shear\-lag theory, where the axial stress\-transfer capacity and structural reinforcement efficiency scale directly with fiber orientation and aspect ratio\. In contrast, cells with fewer load\-aligned fibers exhibited lower initial stiffness and an accelerated onset of matrix\-dominated nonlinearity\.

![Refer to caption](https://arxiv.org/html/2606.28996v1/Figure12.png)Figure 12:Nonlinear mechanical responses for the chosen discrete cells \(ID 1, ID 2, ID 3\)\.Top row: predicted stress–strain curves under \(a\) fiber\-aligned normal loading \(σx​x\\sigma\_\{xx\}–εx​x\\varepsilon\_\{xx\}\), \(b\) transverse normal loading \(σz​z\\sigma\_\{zz\}–εz​z\\varepsilon\_\{zz\}\), and \(c\) axial shear loading \(τx​y\\tau\_\{xy\}–γx​y\\gamma\_\{xy\}\)\.Bottom row: corresponding normalized tangent stiffness degradation expressed as \(d\)Ex​x/Ex​x0E\_\{xx\}/E\_\{xx\}^\{0\}, \(e\)Ez​z/Ez​z0E\_\{zz\}/E\_\{zz\}^\{0\}, and \(f\)Gx​y/Gx​y0G\_\{xy\}/G\_\{xy\}^\{0\}as functions of applied strain\.The normalized tangent stiffness degradation curves in Figure[12](https://arxiv.org/html/2606.28996#S3.F12)d–f \(whereEi​j/Ei​j0E\_\{ij\}/E\_\{ij\}^\{0\}denotes the ratio of instantaneous to initial undamaged stiffness\) reveal three mechanistically significant trends\. First, cells with well\-connected fiber networks maintained their stiffness over a broader strain range, as load redistribution through multiple fiber–fiber interaction pathways remains operative even as matrix damage progresses\. Second, cells characterized by sparse or poorly connected fiber networks exhibited earlier stiffness reduction, consistent with an accelerated transition from fiber\- to matrix\-dominated load\-bearing\. Third, axial shear stiffness degraded more gradually than normal stiffness, reflecting the comparatively stable matrix\-controlled shear transfer mechanism\.

These observations bear close resemblance to connectivity\-driven stiffness transitions in disordered fiber networks, where the density and topology of load\-bearing links dictate the onset and rate of stiffness loss\[[29](https://arxiv.org/html/2606.28996#bib.bib71),[68](https://arxiv.org/html/2606.28996#bib.bib72)\]\. In the present context, locally connected fiber clusters enable load redistribution via redundant pathways, thereby delaying macroscopic stiffness degradation, whereas sparsely connected or resin\-rich cells offer limited redundancy and fail earlier via matrix plasticity and void growth\. These results establish that local fiber\-network connectivity—rather than volume fraction alone—governs the nonlinear evolution of stiffness in AM\-CM SFTs, and that the developed graph\-based surrogate modeling framework accurately captures this physical mechanism\.

### 3\.3Surrogate Prediction of Cell\-Level Constitutive Behavior

Direct nonlinear FE homogenization of every discrete cell is computationally prohibitive for mesoscale assemblies comprising hundreds to thousands of Voronoi cells\. To enable efficient evaluation of cell\-level constitutive behavior while preserving sensitivity to microstructural topology, the hybrid GNN\-LSTM surrogate model described in Section[2\.4](https://arxiv.org/html/2606.28996#S2.SS4)was trained on FE\-generated stress–strain curves\. The training dataset comprised stress–strain curves from 500 discrete Voronoi cells extracted from theμ\\mu\-CT microstructure, each evaluated under six independent loading cases \(three normal and three shear\), capturing a wide range of local fiber volume fractions, fiber orientations, and clustering configurations inherent to the AM\-CM microstructure\. The resulting dataset was partitioned using a stratified random shuffle \(seed 42\) into training \(70%\), validation \(15%\), and testing \(15%\)\.

Training convergence: The training history \(Figure[13](https://arxiv.org/html/2606.28996#S3.F13)\) exhibited two characteristic phases\. A rapid initial reduction in loss was associated with the surrogate learning the topology\-dependent elastic stiffness signatures encoded in the GNN embedding, and the peak stress predictions from the auxiliary head and Random Forest regressor\. This was followed by a more gradual convergence phase in which the LSTM captured the path\-dependent plasticity and post\-peak softening behavior governing the nonlinear damage regime\. The strain\-weighted Huber loss, which assigned up to4×4\\timeshigher weight near and beyond the stress peak, directed the majority of the learning capacity toward accurate representation of the mechanically critical nonlinear regime\. The absence of divergence between training and validation loss throughout both phases confirmed stable generalization without overfitting, supporting the adequacy of the 500\-cell dataset for the microstructural variability encountered\. Early stopping was triggered at convergence, well within the 800\-epoch maximum\.

![Refer to caption](https://arxiv.org/html/2606.28996v1/Figure13_1.png)Figure 13:Training and validation loss histories demonstrating convergence of the GNN–LSTM surrogate for discrete\-cell homogenization under normal loading in the fiber\-dominated direction \(σx​x\\sigma\_\{xx\}\)Predictive accuracy: A parity plot \(Figure[14](https://arxiv.org/html/2606.28996#S3.F14)a–f\) demonstrated near\-perfect agreement between surrogate\-predicted and FE\-computed stresses across the full stress range\. This level of accuracy confirmed that the learned graph representation encoded physically meaningful microstructural information, capturing the topology\-to\-response mapping across the full diversity of fiber arrangements present in theμ\\mu\-CT microstructure—aligned fiber clusters, mixed\-orientation neighborhoods, and resin\-rich regions—rather than relying on globally averaged descriptors that would obscure topology\-driven variability\.

![Refer to caption](https://arxiv.org/html/2606.28996v1/Figure13_2.png)Figure 14:Parity plots comparing the high\-fidelity FE homogenization benchmarks against the predictions generated by the hybrid GNN–LSTM surrogate framework\. Individual coefficient of determination \(R2R^\{2\}\) values are reported across the six independent stress components: \(a\) fiber\-dominated normalσx​x\\sigma\_\{xx\}, \(b\) transverse width normalσy​y\\sigma\_\{yy\}, \(c\) thickness normalσz​z\\sigma\_\{zz\}, \(d\) transverse shearτx​y\\tau\_\{xy\}, \(e\) axial shearτy​z\\tau\_\{yz\}, and \(f\) axial shearτz​x\\tau\_\{zx\}\.![Refer to caption](https://arxiv.org/html/2606.28996v1/Figure14.png)Figure 15:Representative comparison between nonlinear FE homogenization results \(True\) and GNN–LSTM surrogate predictions \(Predicted\) for three discrete Voronoi cells subjected to fiber\-dominated normal loading\. Each panel corresponds to a cell with a distinct local fiber topology: an aligned fiber cluster \(a\), a mixed\-orientation neighborhood \(b\), and a resin\-rich region \(c\)\.Representative comparisons between FE simulations and surrogate predictions across a range of fiber topologies are presented in Figure[15](https://arxiv.org/html/2606.28996#S3.F15)\. The surrogate reproduced both the initial elastic stiffness and the post\-yield nonlinear regime with high fidelity for cells characterized by aligned fiber clusters, mixed\-orientation neighborhoods, and resin\-rich regions, achievingR2=0\.99R^\{2\}=0\.99\. Because all cells shared identical constituent material parameters, the variability in predicted response arose exclusively from differences in the local fiber topology encoded by the GNN—fiber count, spatial arrangement, orientation, distribution, and pairwise interaction strengths captured through the shear\-lag edge weights\. This confirmed that the surrogate learned a topology\-sensitive constitutive mapping rather than a globally averaged response, and that the GNN embeddingzΩcz\_\{\\Omega\_\{c\}\}successfully encoded the mechanically distinguishing features of each cell’s fiber network\.

Computational efficiency: The trained surrogate predicted the complete normalized stress–strain trajectory of an individual cell in under one second on a single CPU, with subsequent de\-normalization via the Random Forest peak\-stress prediction adding negligible computational cost\. The total per\-cell prediction time therefore remained under one second, compared with several minutes required for the equivalent nonlinear FE homogenization running on 5 CPUs in Abaqus 2022\. This reduction in computational cost, by more than two orders of magnitude, enabled efficient evaluation of hundreds to thousands of discrete cells within a mesoscale Voronoi assembly and rendered the proposed multiscale framework tractable for fatigue\-stage\-resolved analyses across multiple specimens and damage states\. Furthermore, unlike FE homogenization where cells are coupled through boundary conditions, surrogate evaluation of individual Voronoi cells is completely independent once the macroscopic strain is localized to each cell, meaning all cells could in principle be evaluated simultaneously in a parallelized implementation, with the theoretical speedup scaling linearly with the number of available processors—an advantage that was not exploited in the present single\-CPU implementation but represents a promising avenue for near\-real\-time coupon\-scale predictions in future HPC deployments\.

### 3\.4Mesoscale Homogenization and Coupon\-Scale Validation

The ultimate purpose of the proposed multiscale framework is to establish a quantitative, predictive link between experimentally observed microstructural degradation and macroscopic mechanical response\. To evaluate this capability, mesoscale Voronoi assemblies were constructed directly fromμ\\mu\-CT reconstructions of fatigue\-tested specimens\. Three specimens, each representing a distinct interrupted fatigue stage \(Section[2](https://arxiv.org/html/2606.28996#S2)\), were selected to capture progressive damage states ranging from nominally undamaged to near\-failure\.

For each specimen and fatigue stage, the fiber centerline coordinates and cell\-level volume fractions extracted from theμ\\mu\-CT reconstruction were passed directly to the GNN\-LSTM surrogate to predict cell\-level constitutive responses, which were assembled into the coupon\-scale response\. No constitutive parameters were recalibrated at the coupon scale between fatigue stages; all observed changes in the predicted macroscopic response arose solely from the microstructure\-informed constitutive update propagated through the LSTM component of the surrogate, in which the cell parameters are uniformly applied across all time points\.

At each fatigue stage, the homogenization procedure incorporatedμ\\mu\-CT\-measured porosity through the two\-step porous\-matrix model described in Section[2\.3\.5](https://arxiv.org/html/2606.28996#S2.SS3.SSS5), which reduced the effective matrix stiffness assigned to all Voronoi cells\. The GNN\-LSTM then predicted the normalized stress trajectory for each cell based on itsμ\\mu\-CT fiber topology, and the Random Forest predicted the corresponding peak stress amplitude from the GNN embedding and local volume fraction\. As fatigue damage advanced, increasing porosity reduced the porous\-matrix stiffness, loss of fiber connectivity, and cluster dispersion were directly captured in lower GNN edge weights and reduced topology embeddings, and the Random Forest consequently predicted lower peak stresses for damaged cells, collectively degrading the assembled coupon response in a manner driven entirely by the evolvingμ\\mu\-CT microstructure rather than by any manually adjusted parameter\.

As illustrated in Appendix[Appendix A3\. Full Six\-Mode Constitutive Response of the Discrete Cell](https://arxiv.org/html/2606.28996#Sx3), the full six\-mode constitutive response is first obtained for each discrete Voronoi cell—either through FE homogenization or the GNN\-LSTM surrogate—and these cell\-level responses are then assembled through volume\-weighted averaging at the mesoscale Voronoi assembly level to produce the final coupon\-scale response\. Figure[16](https://arxiv.org/html/2606.28996#S3.F16)shows the longitudinal mesoscale stress–strain response comparison with experimental monotonic tensile measurements\. With increasing fatigue damage, both experiments and simulations showed a systematic reduction in elastic modulus and an earlier onset of nonlinear deformation\. These trends were reproduced by the multiscale model without any coupon\-scale recalibration, confirming that the observed macroscopic degradation was captured entirely through the microstructure\-informed constitutive update propagated fromμ\\mu\-CT observations through the surrogate\.

The model reproduced the progressive softening of the stress–strain response across the three fatigue conditions\. Quantitatively, predicted elastic moduli agreed with experimental measurements within 5–7% across all three specimens and fatigue stages, confirming that the surrogate captured the correct sensitivity to microstructural degradation across the full range of damage studied, as shown in Figure[16](https://arxiv.org/html/2606.28996#S3.F16)\.

![Refer to caption](https://arxiv.org/html/2606.28996v1/Figure15.png)Figure 16:Comparison of experimentally measured and multiscale\-predicted stress–strain responses obtained from the hybrid FE\-surrogate framework, in which GNN\-LSTM surrogate predictions of discrete Voronoi cell constitutive responses are assembled through mesoscale periodic boundary conditions and volume\-weighted averaging to predict the coupon\-scale response, for all three specimens at each fatigue stage \(BT: before test, BF: below fatigue limit, AF: above fatigue limit\)\. Top \(a\), middle \(b\), and bottom \(c\) panels correspond to Specimen 1, Specimen 2, and Specimen 3, respectively\.The discrepancies between simulation and experiment were attributed primarily to physical mechanisms not represented in the current constitutive framework, rather than to surrogate approximation error\. The matrix model captures elastoplastic damage and porosity\-induced stiffness reduction, but does not explicitly resolve fiber\-matrix interfacial debonding or discrete matrix cracking\. Both become increasingly active at advanced stages of fatigue, causing further loss of stiffness in highly damaged regions\. The surrogate itself contributed negligibly to this error, as confirmed by comparing coupon\-scale predictions from the GNN\-LSTM surrogate with those obtained from direct FE homogenization of the same Voronoi assemblies\. These agree closely across all fatigue stages, isolating the source of the experiment\-simulation gap to the constitutive model rather than to the surrogate model\.

Taken together, these results demonstrated that the proposed FE\-Voronoi\-GNN framework successfully propagated microstructural information across three length scales—from the explicitly resolved fiber level through the topology\-sensitive surrogate to the macroscopic coupon response—using onlyμ\\mu\-CT fiber coordinates and volume fractions as inputs, with no specialized characterization beyond standard reconstruction\. By directly linking manufacturing\-induced and fatigue\-driven microstructural evolution to measurable changes in mechanical performance, the framework provided a physically grounded and computationally efficient basis for assessing the structural condition and remaining fatigue life of AM\-CM short\-fiber thermoplastic composite components\.

It is worth noting that while a mean\-field RVE approach using average microstructural descriptors may yield comparable elastic modulus predictions at the coupon scale, the present framework offers capabilities that such an approach fundamentally cannot provide: \(i\) spatially resolved identification of mechanically weak cells governed by local fiber clustering and sparse connectivity, \(ii\) direct propagation of fatigue\-stage\-specific microstructural evolution without constitutive recalibration, and \(iii\) capture of topology\-driven variability in nonlinear damage evolution that is entirely lost under mean\-field averaging\.

## 4Conclusions and Future Work

This study established a microstructure\-resolved multiscale surrogate modeling framework to predict the nonlinear mechanical response of SFT composites manufactured via AM\-CM\. By combining interaction\-topology\-based discrete homogenization with a physics\-informed GNN\-LSTM surrogate, the framework provided an efficient, physically interpretable pathway for linking local fiber topology directly to macroscopic structural performance\.

The key findings are summarized below\.

- •Microstructure\-controlled mechanics in discrete cells:High\-resolutionμ\\mu\-CT reconstructions demonstrate that the AM\-CM microstructure can be decomposed into localized Voronoi cells representing neighborhoods of interacting fibers\. Even when global fiber volume fraction and constituent properties remain unchanged, these cells exhibit substantial variability in stiffness and nonlinear response arising solely from differences in fiber clustering, alignment, and spacing\. The mechanical behavior of SFTs is therefore governed not only by volume fraction but also by the topology of local fiber interaction networks that control load\-transfer pathways between fibers and matrix\. Nonlinear FE homogenization showed that cells containing aligned or highly connected fiber clusters sustained load primarily through axial fiber stretching, yielding higher stiffness and delayed degradation, whereas cells with dispersed or sparsely connected configurations relied more on matrix\-mediated load transfer, producing earlier nonlinear deformation and accelerated stiffness reduction\.
- •Physics\-informed surrogate modeling:A hybrid GNN\-LSTM architecture was developed and trained to approximate the nonlinear constitutive response of discrete cells\. The surrogate required 3D fiber centerline endpoints and cell\-level volume fractions as inputs, both of which are available directly from standardμ\\mu\-CT reconstruction\. By representing the fiber interaction network as a graph with physics\-informed edge weights that encode both axial shear\-lag coupling and the orientation\-independent transverse stress\-concentration mechanism, as confirmed by parametric two\-fiber FEA, the GNN captured topology\-dependent stiffness\. The LSTM learned the path\-dependent normalized constitutive behavior; a skip connection from the GNN embedding to the output network ensured topology\-encoded amplitude information directly modulated stress predictions at every loading increment\. This two\-stage prediction pipeline reproduced FE\-predicted stress–strain responses withR2≈0\.98R^\{2\}\\approx 0\.98on de\-normalized physical stress trajectories while reducing per\-cell evaluation time by more than two orders of magnitude, enabling near\-real\-time prediction over hundreds to thousands of Voronoi cells\.
- •Multiscale propagation of microstructural degradation\.When integrated within mesoscale Voronoi assemblies derived directly fromμ\\mu\-CT reconstructions at each fatigue stage, the surrogate\-enabled framework accurately predicted coupon\-level tensile responses across multiple damage states and three specimens without recalibrating any constitutive parameters\. Predicted elastic moduli agreed with experiment within 5–7%, and area under the stress–strain curve within 10%, confirming that microstructure\-driven degradation mechanisms were correctly propagated across scales through the surrogate alone, driven entirely by the evolvingμ\\mu\-CT microstructure\.

Collectively, these results established that the mechanical response of AM\-CM SFTs was fundamentally governed by the topology of local fiber interaction networks rather than by global microstructural averages alone\. The proposed FE–Voronoi–GNN framework established a physically grounded bridge between microstructural evolution and macroscopic mechanical performance, deployable on anyμ\\mu\-CT dataset from which fiber centerlines and volume fractions could be extracted, with no specialized inputs beyond standard reconstruction\.

Beyond the specific CF\-ABS system studied, the framework provides a scalable foundation for microstructure\-aware digital twins of discontinuous fiber composites\. The ability to predict constitutive response rapidly and directly from minimal microstructural inputs opens opportunities for accelerated materials design, process optimization, and structural performance assessment in additively manufactured composite systems\.

Future work will focus on: \(i\) incorporating explicit fiber\-matrix interfacial debonding mechanisms to capture intra\-cell failure modes not currently represented; \(ii\) embedding thermodynamic admissibility constraints directly within the GNN\-LSTM architecture to enforce energetic consistency for more generalized global loading histories, including cyclic fatigue loading and non\-proportional multiaxial loading paths at the coupon scale, beyond the monotonic quasi\-static regime studied here; and \(iii\) extending the framework to other thermoplastic composite systems and processing routes, including injection molding and multi\-material additive manufacturing\.

## 5Acknowledgments

This research was supported by the NASA 80NSSC24M0154 Grant\. This research was also supported in part by the U\.S\. Department of Energy Office of Energy Efficiency and Renewable Energy \(EERE\), Advanced Materials and Manufacturing Technologies Office \(AMMTO\)\.

The authors would like to express their sincere gratitude to Brandon L\. Hearley at the NASA Glenn Research Center for his valuable review, insights, and assistance in the preparation and writing of this manuscript\.

The authors utilized AI\-based tools to assist with proofreading and refining the manuscript’s clarity\. Specifically, Claude and Grammarly were employed to improve grammar, sentence flow, and readability\. No generative AI tools were used for data analysis, interpretation of results, or drawing scientific conclusions\. The responsibility for the content and final approval of the manuscript rests entirely with the authors\.

## References

- \[1\]N\. Abhilash and G\. Suhasini\(2020\)Effect of chemical vapor infiltration induced matrix porosity on the mechanical behavior of ceramic matrix minicomposites\.ASCE\-ASME J Risk and Uncert in Engrg Sys Part B Mech Engrg6,pp\. 041005\.External Links:[Document](https://dx.doi.org/10.1115/1.4047465)Cited by:[§2\.3\.5](https://arxiv.org/html/2606.28996#S2.SS3.SSS5.p1.2),[§3\.1](https://arxiv.org/html/2606.28996#S3.SS1.p5.2)\.
- \[2\]J\. H\. S\. Almeida, A\. Miettinen, F\. Léonard, B\. G\. Falzon, and P\. J\. Withers\(2025\)Microstructure and damage evolution in short carbon fibre 3d\-printed composites during tensile straining\.Composites Part B: Engineering292,pp\. 112073\.External Links:ISSN 1359\-8368,[Document](https://dx.doi.org/10.1016/j.compositesb.2024.112073)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p1.1)\.
- \[3\]S\. F\. Almeida and Z\. S\. Neto\(1994\)Effect of void content on the strength of composite laminates\.Composite Structures28,pp\. 139–148\.External Links:ISSN 0263\-8223,[Document](https://dx.doi.org/10.1016/0263-8223%2894%2990044-2)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p2.1),[§3\.1](https://arxiv.org/html/2606.28996#S3.SS1.p3.1)\.
- \[4\]L\. Anqi, J\.J\.C\. Remmers, A\.W\.J\. V\. Dommelen, J\.T\. Massart, and G\.D\.M\. Geers\(2025\)A mean\-field homogenization model for fiber reinforced composite materials in large deformation with plasticity\.International Journal of Solids and Structures310,pp\. 113200\.External Links:ISSN 0020\-7683,[Document](https://dx.doi.org/10.1016/j.ijsolstr.2024.113200)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p4.1)\.
- \[5\]T\. Arabatti, N\. K\. Parambil, and S\. Gururaja\(2016\)Micromechanical modeling of damage development in polymer composites\.Advanced Composites Letters25,pp\. 096369351602500301\.External Links:[Document](https://dx.doi.org/10.1177/096369351602500301)Cited by:[§2\.3\.3](https://arxiv.org/html/2606.28996#S2.SS3.SSS3.p3.1)\.
- \[6\]M\.F\. Arif, N\. Saintier, F\. Meraghni, J\. Fitoussi, Y\. Chemisky, and G\. Robert\(2014\)Multiscale fatigue damage characterization in short glass fiber reinforced polyamide\-66\.Composites Part B: Engineering61,pp\. 55–65\(en\)\.External Links:ISSN 13598368,[Document](https://dx.doi.org/10.1016/j.compositesb.2014.01.019)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p2.1),[§1](https://arxiv.org/html/2606.28996#S1.p4.1),[§3\.1](https://arxiv.org/html/2606.28996#S3.SS1.p3.1)\.
- \[7\]F\. Bandinelli, L\. Peroni, and A\. Morena\(2023\)Elasto\-plastic mechanical modeling of fused deposition 3d printing materials\.Polymers15,pp\. 234\.External Links:[Document](https://dx.doi.org/10.3390/polym15010234)Cited by:[§2\.3\.4](https://arxiv.org/html/2606.28996#S2.SS3.SSS4.p4.1)\.
- \[8\]Y\. Bengio\(2012\)Practical recommendations for gradient\-based training of deep architectures\.Vol\.7700,Springer Berlin Heidelberg\.External Links:[Document](https://dx.doi.org/10.1007/978-3-642-35289-8%5F26)Cited by:[§2\.4\.4](https://arxiv.org/html/2606.28996#S2.SS4.SSS4.p15.3)\.
- \[9\]K\. M\. M\. Billah, V\. Kumar, N\. Rathod, A\. Phadatare, D\. Nuttall, T\. Smith, U\. Vaidya, S\. Kim, V\. Kunc, and A\. A\. Hassen\(2024\)Experimental investigation of short fiber reinforced thermoplastic composite joining in an integrated additive manufacturing and compression molding system\.Additive Manufacturing; Advanced Materials Manufacturing; Biomanufacturing; Life Cycle Engineering1,pp\. V001T02A013\.External Links:[Document](https://dx.doi.org/10.1115/MSEC2024-125490)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p1.1),[§2\.1](https://arxiv.org/html/2606.28996#S2.SS1.p3.1)\.
- \[10\]K\. Breuer and M\. Stommel\(2021\)Prediction of short fiber composite properties by an artificial neural network trained on an rve database\.Fibers9,pp\. 8\.External Links:ISSN 2079\-6439,[Document](https://dx.doi.org/10.3390/fib9020008)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p4.1),[§1](https://arxiv.org/html/2606.28996#S1.p7.1)\.
- \[11\]S\. Cassola, M\. Duhovic, T\. Schmidt, and D\. May\(2022\)Machine learning for polymer composites process simulation – a review\.Composites Part B: Engineering246,pp\. 110208\.External Links:[Document](https://dx.doi.org/10.1016/j.compositesb.2022.110208)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p7.1)\.
- \[12\]C\. Chen and G\. X\. Gu\(2019\)Machine learning for composite materials\.MRS Communications9,pp\. 556–566\.External Links:[Document](https://dx.doi.org/10.1557/mrc.2019.32)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p7.1)\.
- \[13\]Q\. Chen, R\. Jia, and S\. Pang\(2021\)Deep long short\-term memory neural network for accelerated elastoplastic analysis of heterogeneous materials: an integrated data\-driven surrogate approach\.Composite Structures264,pp\. 113688\.External Links:[Document](https://dx.doi.org/10.1016/j.compstruct.2021.113688)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p7.1)\.
- \[14\]F\. Chinesta, P\. Ladeveze, and E\. Cueto\(2011\)A Short Review on Model Order Reduction Based on Proper Generalized Decomposition\.Archives of Computational Methods in Engineering18,pp\. 395–404\.External Links:ISSN 1886\-1784,[Document](https://dx.doi.org/10.1007/s11831-011-9064-7)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p6.1)\.
- \[15\]H\. L\. Cox\(1952\-03\)The elasticity and strength of paper and other fibrous materials\.British Journal of Applied Physics3,pp\. 72\.External Links:[Document](https://dx.doi.org/10.1088/0508-3443/3/3/302)Cited by:[§2\.4\.1](https://arxiv.org/html/2606.28996#S2.SS4.SSS1.p5.6)\.
- \[16\]M\. Dai, M\. F\. Demirel, Y\. Liang, and J\. Hu\(2021\)Graph neural networks for an accurate and interpretable prediction of the properties of polycrystalline materials\.Computational Materials7,pp\. 103\(en\)\.External Links:ISSN 2057\-3960,[Document](https://dx.doi.org/10.1038/s41524-021-00574-w)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p6.1),[§1](https://arxiv.org/html/2606.28996#S1.p7.1)\.
- \[17\]M\. Donal\(1997\)Short fiber orientation and its effects on the properties of thermoplastic composite materials\.\.Polymer\-Plastics Technology and Engineering8,pp\. 101–54\(en\)\.External Links:[Document](https://dx.doi.org/10.1080/03602557708545033)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p1.1),[§1](https://arxiv.org/html/2606.28996#S1.p2.1)\.
- \[18\]V\. E\.\(1948\)The relationship between stress and strain for homogeneous deformation\.Journal of the Institute of Metals74,pp\. 537–562\.Cited by:[§2\.3\.4](https://arxiv.org/html/2606.28996#S2.SS3.SSS4.p4.1)\.
- \[19\]A\. Estefani and L\. Távara\(2024\)Numerical multiscale analysis of 3D printed short fiber composites parts: Filament anisotropy and toolpath effects\.Engineering Reports6,pp\. 12799\(en\)\.External Links:ISSN 2577\-8196, 2577\-8196,[Document](https://dx.doi.org/10.1002/eng2.12799)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p4.1)\.
- \[20\]B\. Eyri, O\. Gul, N\. G\. Karsli, and T\. Yilmaz\(2025\)Modeling short fiber reinforced polymer matrix composite materials using material designer\.Polymer Composites46,pp\. 10350–10360\(en\)\.External Links:ISSN 0272\-8397, 1548\-0569,[Document](https://dx.doi.org/10.1002/pc.29625)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p6.1)\.
- \[21\]S\. Ferdousi, Z\. Demchuk, W\. Choi, R\. C\. Advincula, and Y\. Jiang\(2025\)A deep learning and finite element approach for exploration of inverse structure–property designs of lightweight hybrid composites\.Composite Structures365,pp\. 119179\.External Links:[Document](https://dx.doi.org/10.1016/j.compstruct.2025.119179)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p7.1)\.
- \[22\]S\. Y\. Fu, B\. Lauke, Y\. H\. Zhang, and Y\.\-W\. Mai\(2005\)On the post\-mortem fracture surface morphology of short fiber reinforced thermoplastics\.Composites Part A: Applied Science and Manufacturing36,pp\. 987–994\.External Links:ISSN 1359\-835X,[Document](https://dx.doi.org/10.1016/j.compositesa.2004.11.005)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p2.1)\.
- \[23\]U\. Gandhi, D\. B\. Sebastian, V\. Kunc, and Y\. Song\(2016\)Method to measure orientation of discontinuous fiber embedded in the polymer matrix from computerized tomography scan data\.Journal of Thermoplastic Composite Materials29,pp\. 1696–1709\(en\)\.External Links:ISSN 0892\-7057, 1530\-7980,[Document](https://dx.doi.org/10.1177/0892705715584411)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p1.1)\.
- \[24\]J\. Ge, X\. Chao, H\. Hu, W\. Tian, W\. Li, and L\. Qi\(2024\)Numerical Assessment of Effective Elastic Properties of Needled Carbon/Carbon Composites Based on a Multiscale Method\.Journal of Carbon Research10,pp\. 85\(en\)\.External Links:ISSN 2311\-5629,[Document](https://dx.doi.org/10.3390/c10030085)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p7.1)\.
- \[25\]J\. Ghaboussi, J\. H\. Garrett, and X\. Wu\(1991\)Knowledge‐based modeling of material behavior with neural networks\.Journal of Engineering Mechanics117,pp\. 132–153\.External Links:[Document](https://dx.doi.org/10.1061/%28ASCE%290733-9399%281991%29117%3A1%28132%29)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p7.1)\.
- \[26\]D\. E\. Gulmez, J\. Maldonado, K\. Masania, J\. Sinke, and C\. Dransfeld\(2023\)Quantification of structural response and edge orientation of Chopped Tape Thermoplastic Composites in net\-shaped specimens\.Composite Structures321,pp\. 117302\(en\)\.External Links:ISSN 02638223,[Document](https://dx.doi.org/10.1016/j.compstruct.2023.117302)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p6.1)\.
- \[27\]J\. Gunst, E\. Vaghefi, E\. Mirkoohi, G\. Deshpande, B\. Prorok, and S\. Gururaja\(2025\)Leveraging machine learning to quantify the effect of porosity on mode i delamination in ceramic matrix composites\.Composite Structures372,pp\. 119542\.External Links:[Document](https://dx.doi.org/10.1016/j.compstruct.2025.119542)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p7.1)\.
- \[28\]M\. H\. Hamza, L\. Borkowski, and A\. Chattopadhyay\(2025\)Physics\-constrained machine learning surrogate model for time\-dependent behavior of ceramic matrix composites\.Composites Part B: Engineering307,pp\. 112825\.External Links:ISSN 1359\-8368,[Document](https://dx.doi.org/10.1016/j.compositesb.2025.112825)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p6.1)\.
- \[29\]D\. A\. Head, A\. J\. Levine, and F\. C\. MacKintosh\(2003\)Deformation of cross\-linked semiflexible polymer networks\.Phys\. Rev\. Lett\.91,pp\. 108102\.External Links:[Document](https://dx.doi.org/10.1103/PhysRevLett.91.108102)Cited by:[§3\.2](https://arxiv.org/html/2606.28996#S3.SS2.p6.1)\.
- \[30\]J\. Horst and J\. Spoormaker\(1997\)Fatigue fracture mechanism and fractography of short\-glassfibre\-reinforced polyamide 6\.Journal of Materials Science32,pp\. 3641–3651\.External Links:[Document](https://dx.doi.org/10.1023/A%3A1018634530869)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p2.1)\.
- \[31\]L\. Jiang, Y\. Zhou, and F\. Jin\(2022\)Design of short fiber\-reinforced thermoplastic composites: a review\.Polymer Composites43,pp\. 4835–4847\.External Links:[Document](https://dx.doi.org/10.1002/pc.26817)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p7.1)\.
- \[32\]J\. M\. Kaiser and M\. Stommel\(2012\)An extended mean\-field homogenization model to predict the strength of short\-fibre polymer composites\(Website\)Note:Accessed: Dec\. 29, 2025\. \[Online\]\. Available:[https://journals\.ub\.ovgu\.de/index\.php/techmech/article/view/725](https://journals.ub.ovgu.de/index.php/techmech/article/view/725),Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p4.1)\.
- \[33\]S\. Kammoun, I\. Doghri, L\. Brassart, and L\. Delannay\(2015\)Micromechanical modeling of the progressive failure in short glass–fiber reinforced thermoplastics – first pseudo\-grain damage model\.Composites Part A: Applied Science and Manufacturing73,pp\. 166–175\.External Links:ISSN 1359\-835X,[Document](https://dx.doi.org/10.1016/j.compositesa.2015.02.017)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p4.1)\.
- \[34\]F\. Kibrete, T\. Trzepieciński, H\. S\. Gebremedhen, and D\. E\. Woldemichael\(2023\)Artificial intelligence in predicting mechanical properties of composite materials\.Journal of Composites Science7,pp\. 364\.External Links:[Document](https://dx.doi.org/10.3390/jcs7090364)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p7.1)\.
- \[35\]J\. Köbler, M\. Schneider, F\. Ospald, H\. Andrä, and R\. Müller\(2018\)Fiber orientation interpolation for the multiscale analysis of short fiber reinforced composite parts\.Computational Mechanics61,pp\. 729–750\(en\)\.External Links:ISSN 0178\-7675, 1432\-0924,[Document](https://dx.doi.org/10.1007/s00466-017-1478-0)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p6.1)\.
- \[36\]A\. Krairi and I\. Doghri\(2014\)A thermodynamically\-based constitutive model for thermoplastic polymers coupling viscoelasticity, viscoplasticity and ductile damage\.International Journal of Plasticity60,pp\. 163–181\.External Links:[Document](https://dx.doi.org/10.1016/j.ijplas.2014.04.010)Cited by:[§2\.3\.4](https://arxiv.org/html/2606.28996#S2.SS3.SSS4.p6.2),[§2\.3\.4](https://arxiv.org/html/2606.28996#S2.SS3.SSS4.p8.1)\.
- \[37\]V\. Kumar, S\. P\. Alwekar, V\. Kunc, E\. Cakmak, V\. Kishore, T\. Smith, J\. Lindahl, U\. Vaidya, C\. Blue, M\. Theodore, S\. Kim, and A\. A\. Hassen\(2021\)High\-performance molded composites using additively manufactured preforms with controlled fiber and pore morphology\.Additive Manufacturing37,pp\. 101733\.External Links:ISSN 2214\-8604,[Document](https://dx.doi.org/10.1016/j.addma.2020.101733)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p1.1)\.
- \[38\]V\. Kumar, A\. A\. Hassen, V\. Kunc, D\. Nuttall, A\. Kircaliali, S\. I\. Talabi, J\. Reynolds, J\. Vaughan, P\. Mhatre, T\. Smith, B\. Šeta, J\. Spangenberg, and C\. Blue\(2026\)Genesis of a novel high\-rate composite manufacturing process using large\-scale additive manufacturing – compression molding \(am\-cm\) system: possibilities and limitations∗\.Composites Part B: Engineering309,pp\. 113101\.External Links:ISSN 1359\-8368,[Document](https://dx.doi.org/10.1016/j.compositesb.2025.113101)Cited by:[§2\.1](https://arxiv.org/html/2606.28996#S2.SS1.p3.1)\.
- \[39\]C\. M\. Landis, M\. A\. McGlockton, and R\. M\. McMeeking\(1999\)An improved shear lag model for broken fibers in composite materials\.Journal of Composite Materials33,pp\. 667–680\.External Links:[Document](https://dx.doi.org/10.1177/002199839903300704)Cited by:[§2\.4\.1](https://arxiv.org/html/2606.28996#S2.SS4.SSS1.p5.6)\.
- \[40\]B\. Liu, E\. Ocegueda, M\. Trautner, A\. M\. Stuart, and K\. Bhattacharya\(2023\)Learning macroscopic internal variables and history dependence from microscopic models\.Journal of the Mechanics and Physics of Solids178,pp\. 105329\.External Links:ISSN 0022\-5096,[Document](https://dx.doi.org/10.1016/j.jmps.2023.105329)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p7.1)\.
- \[41\]Z\. Liu, Y\. Lei, X\. Zhang, Z\. Kang, and J\. Zhang\(2022\)Effect mechanism and simulation of voids on hygrothermal performances of composites\.Polymers14\.External Links:ISSN 2073\-4360,[Document](https://dx.doi.org/10.3390/polym14050901)Cited by:[§3\.1](https://arxiv.org/html/2606.28996#S3.SS1.p5.2)\.
- \[42\]I\. Loshchilov and F\. Hutter\(2019\)Decoupled weight decay regularization\.External Links:1711\.05101Cited by:[§2\.4\.4](https://arxiv.org/html/2606.28996#S2.SS4.SSS4.p15.3)\.
- \[43\]D\.B\. Marshall, B\.N\. Cox, and A\.G\. Evans\(1985\)The mechanics of matrix cracking in brittle\-matrix fiber composites\.\.Acta Metallurgica33,pp\. 2013–2021\.External Links:[Document](https://dx.doi.org/10.1016/0001-6160%2885%2990124-5)Cited by:[§2\.4\.1](https://arxiv.org/html/2606.28996#S2.SS4.SSS1.p5.6)\.
- \[44\]D\.B\. Marshall and B\.N\. Cox\(1987\)Tensile fracture of brittle matrix composites: influence of fiber strength\.Acta Metallurgica35,pp\. 2607–2619\.External Links:ISSN 0001\-6160,[Document](https://dx.doi.org/10.1016/0001-6160%2887%2990260-4)Cited by:[§2\.4\.1](https://arxiv.org/html/2606.28996#S2.SS4.SSS1.p5.6)\.
- \[45\]N\. Mentges, H\. Çelik, C\. Hopmann, M\. Fagerström, and S\. M\. Mirkhalaf\(2023\)Micromechanical modelling of short fibre composites considering fibre length distributions\.Composites Part B: Engineering264,pp\. 110868\.External Links:ISSN 1359\-8368,[Document](https://dx.doi.org/10.1016/j.compositesb.2023.110868)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p6.1)\.
- \[46\]L\. Mishnaevsky and P\. Brøndsted\(2009\)Micromechanical modeling of damage and fracture of unidirectional fiber reinforced composites: a review\.Computational Materials Science44,pp\. 1351–1359\.External Links:ISSN 0927\-0256,[Document](https://dx.doi.org/10.1016/j.commatsci.2008.09.004)Cited by:[§2\.4\.1](https://arxiv.org/html/2606.28996#S2.SS4.SSS1.p5.6)\.
- \[47\]S\. Mortazavian and A\. Fatemi\(2015\)Fatigue behavior and modeling of short fiber reinforced polymer composites: a literature review\.International Journal of Fatigue70,pp\. 297–321\.External Links:ISSN 0142\-1123,[Document](https://dx.doi.org/10.1016/j.ijfatigue.2014.10.005)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p2.1)\.
- \[48\]M\. Mozaffar, R\. Bostanabad, W\. Chen, K\. Ehmann, J\. Cao, and M\. A\. Bessa\(2019\)Deep learning predicts path\-dependent plasticity\.Proceedings of the National Academy of Sciences116,pp\. 26414–26420\.External Links:[Document](https://dx.doi.org/10.1073/pnas.1911815116)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p7.1)\.
- \[49\]V\. Müller, M\. Kabel, H\. Andrä, and T\. Böhlke\(2015\)Homogenization of linear elastic properties of short\-fiber reinforced composites – A comparison of mean field and voxel\-based methods\.International Journal of Solids and Structures67\-68,pp\. 56–70\.External Links:ISSN 0020\-7683,[Document](https://dx.doi.org/10.1016/j.ijsolstr.2015.02.030)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p4.1)\.
- \[50\]A\. Nagaraja, A\. BR, C\. Ramesh, D\. Ghosh, and S\. Gururaja\(2024\)Multiscale modeling for mechanical property estimation of polymer matrix composites with nano\- and micro\-scale reinforcements\.Journal of Reinforced Plastics and Composites45,pp\. 1001–1021\.External Links:[Document](https://dx.doi.org/10.1177/07316844241291718)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p4.1)\.
- \[51\]J\. A\. Nairn\(1997\)On the use of shear\-lag methods for analysis of stress transfer in unidirectional composites\.Mechanics of Materials26,pp\. 63–80\.External Links:ISSN 0167\-6636,[Document](https://dx.doi.org/10.1016/S0167-6636%2897%2900023-9)Cited by:[§2\.4\.1](https://arxiv.org/html/2606.28996#S2.SS4.SSS1.p5.6),[§3\.1](https://arxiv.org/html/2606.28996#S3.SS1.p5.2)\.
- \[52\]P\. Parandoush and D\. Lin\(2017\)A review on additive manufacturing of polymer\-fiber composites\.Composite Structures182,pp\. 36–53\.External Links:ISSN 0263\-8223,[Document](https://dx.doi.org/10.1016/j.compstruct.2017.08.088)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p1.1)\.
- \[53\]C\.H\. Park and W\.I\. Lee\(2012\)3 \- compression molding in polymer matrix composites\.InManufacturing Techniques for Polymer Matrix Composites \(PMCs\),Woodhead Publishing Series in Composites Science and Engineering,pp\. 47–94\.External Links:ISBN 978\-0\-85709\-067\-6,[Document](https://dx.doi.org/10.1533/9780857096258.1.47)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p1.1)\.
- \[54\]P\. Pathak, S\. Gururaja, V\. Kumar, D\. Nuttall, A\. Mahmoudi, M\.M\. Khonsari, and U\. Vaidya\(2025\)Examining infrared thermography based approaches to rapid fatigue characterization of additively manufactured compression molded short fiber thermoplastic composites\.Composites Structure351,pp\. 118610\.External Links:ISSN 0263\-8223,[Document](https://dx.doi.org/10.1016/j.compstruct.2024.118610)Cited by:[§2\.1](https://arxiv.org/html/2606.28996#S2.SS1.p1.1),[§2\.1](https://arxiv.org/html/2606.28996#S2.SS1.p2.1)\.
- \[55\]P\. Pathak and S\. Gururaja\(2025\)Porosity evolution in additively manufactured compression molded short fiber thermoplastics under cyclic loading: insights from micro\-computed tomography and infrared thermography\.Composites Part A: Applied Science and Manufacturing194,pp\. 108881\.External Links:ISSN 1359\-835X,[Document](https://dx.doi.org/10.1016/j.compositesa.2025.108881)Cited by:[Figure 1](https://arxiv.org/html/2606.28996#S1.F1),[Figure 1](https://arxiv.org/html/2606.28996#S1.F1.4.2),[§1](https://arxiv.org/html/2606.28996#S1.p1.1),[§1](https://arxiv.org/html/2606.28996#S1.p3.1),[Figure 2](https://arxiv.org/html/2606.28996#S2.F2),[Figure 2](https://arxiv.org/html/2606.28996#S2.F2.2.1),[Figure 3](https://arxiv.org/html/2606.28996#S2.F3),[Figure 3](https://arxiv.org/html/2606.28996#S2.F3.2.1),[§2\.1](https://arxiv.org/html/2606.28996#S2.SS1.p1.1),[§2\.1](https://arxiv.org/html/2606.28996#S2.SS1.p2.1),[§2\.1](https://arxiv.org/html/2606.28996#S2.SS1.p3.1),[§2\.2](https://arxiv.org/html/2606.28996#S2.SS2.p4.5),[§2\.3\.3](https://arxiv.org/html/2606.28996#S2.SS3.SSS3.p3.1),[§2\.3\.4](https://arxiv.org/html/2606.28996#S2.SS3.SSS4.p1.1),[§3\.1](https://arxiv.org/html/2606.28996#S3.SS1.p5.2),[Table 3](https://arxiv.org/html/2606.28996#S3.T3),[Table 3](https://arxiv.org/html/2606.28996#S3.T3.12.2)\.
- \[56\]P\. Pathak, K\. F\. Bristy, V\. Kumar, K\. M\. M\., and S\. Gururaja\(2026\)On the influence of cyclic loading frequency on fatigue limit of short‐fiber thermoplastics\.Fatigue and Fracture of Engineering Materials and Structures45,pp\. 1001–1021\.External Links:[Document](https://dx.doi.org/10.1111/ffe.70241)Cited by:[§2\.1](https://arxiv.org/html/2606.28996#S2.SS1.p1.1),[§2\.1](https://arxiv.org/html/2606.28996#S2.SS1.p2.1)\.
- \[57\]D\. K\. Pokkalla, V\. Kumar, E\. Jo, A\. A\. Hassen, E\. Cakmak, S\. Alwekar, V\. Kunc, U\. Vaidya, H\. Baid, and S\. Kim\(2022\)Anisotropic mechanical properties of polymer composites from a hybrid additive manufacturing\-compression molding process using x\-ray computer tomography\.InNondestructive Characterization and Monitoring of Advanced Materials, Aerospace, Civil Infrastructure, and Transportation XVI,pp\. 46\.External Links:[Document](https://dx.doi.org/10.1117/12.2614620)Cited by:[§3\.1](https://arxiv.org/html/2606.28996#S3.SS1.p5.2)\.
- \[58\]J\.R\. Pothnis, S\. Gururaja, and D\. Kalyanasundaram\(2019\)Development and characterization of electric field directed preferentially aligned cnt nanocomposites\.Mechanics of Advanced Materials and Structures26,pp\. 35–41\.External Links:[Document](https://dx.doi.org/10.1080/15376494.2018.1534165)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p4.1)\.
- \[59\]F\. Rezaei, R\. Yunus, N\. A\. Ibrahim, and E\. S\. Mahdi\(2008\)Development of Short\-Carbon\-Fiber\-Reinforced Polypropylene Composite for Car Bonnet\.Polymer\-Plastics Technology and Engineering47,pp\. 351–357\(en\)\.External Links:ISSN 0360\-2559, 1525\-6111,[Document](https://dx.doi.org/10.1080/03602550801897323)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p1.1)\.
- \[60\]L\. Schraa, C\. Rodricks, G\. Kalinka, K\. Roetsch, C\. Scheffler, A\. Sambale, K\. Uhlig, M\. Stommel, and V\. Trappe\(2025\)Characterisation and modelling of the fibre\-matrix interface of short fibre reinforced thermoplastics using the push\-out technique\.Composites Part B: Engineering297,pp\. 112317\.External Links:ISSN 1359\-8368,[Document](https://dx.doi.org/10.1016/j.compositesb.2025.112317)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p6.1)\.
- \[61\]A\. Senthilnathan, V\. Saseendran, P\. Acar, N\. Yamamoto, and V\. Sundararaghavan\(2024\)Comparison and validation of stochastic microstructure characterization and reconstruction: machine learning vs\. deep learning methodologies\.Acta Materialia278,pp\. 120220\(en\)\.External Links:ISSN 1359–6454,[Document](https://dx.doi.org/10.1016/j.actamat.2024.120220)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p6.1),[§1](https://arxiv.org/html/2606.28996#S1.p7.1)\.
- \[62\]H\. Sharma, G\. Arora, M\. K\. Singh, V\. Ayyappan, P\. Bhowmik, S\. M\. Rangappa, and S\. Siengchin\(2025\)Review of machine learning approaches for predicting mechanical behavior of composite materials\.Discover Applied Sciences7,pp\. 1238\.External Links:[Document](https://dx.doi.org/10.1007/s42452-025-07616-8)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p7.1)\.
- \[63\]L\. Shuguang and W\. Anchana\(2004\)Unit cells for micromechanical analyses of particle\-reinforced composites\.Mechanics of Materials36,pp\. 543–572\.External Links:ISSN 0167\-6636,[Document](https://dx.doi.org/10.1016/S0167-6636%2803%2900062-0)Cited by:[§2\.3\.3](https://arxiv.org/html/2606.28996#S2.SS3.SSS3.p3.1)\.
- \[64\]A\. Stamopoulos, K\. Tserpes, P\. Prucha, and D\. Vavrik\(2016\)Evaluation of porosity effects on the mechanical properties of carbon fiber\-reinforced plastic unidirectional laminates by x\-ray computed tomography and mechanical testing\.Journal of Composite Materials50,pp\. 2087–98\.External Links:ISSN 0021\-9983,[Document](https://dx.doi.org/10.1177/0021998315602049)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p2.1)\.
- \[65\]J\. Summerscales, F\. J\. Guild, N\. R\. L\. Pearce, and P\. M\. Russell\(2001\)Voronoi cells, fractal dimensions and fibre composites\.Journal of Microscopy201,pp\. 153–162\(en\)\.External Links:ISSN 0022\-2720, 1365\-2818,[Document](https://dx.doi.org/10.1046/j.1365-2818.2001.00841.x)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p6.1)\.
- \[66\]S\. Tamboura, M\.A\. Laribi, J\. Fitoussi, M\. Shirinbayan, R\. T\. Bi, A\. Tcharkhtchi, and H\. B\. Dali\(2020\)Damage and fatigue life prediction of short fiber reinforced composites submitted to variable temperature loading: application to sheet molding compound composites\.International Journal of Fatigue138,pp\. 105676\.External Links:ISSN 0142\-1123,[Document](https://dx.doi.org/10.1016/j.ijfatigue.2020.105676)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p2.1)\.
- \[67\]H\. L\. Tekinalp, V\. Kunc, G\. M\. Velez\-Garcia, C\. E\. Duty, L\. J\. Love, A\. K\. Naskar, C\. A\. Blue, and S\. Ozcan\(2014\)Highly oriented carbon fiber–polymer composites via additive manufacturing\.Composites Science and Technology105,pp\. 144–150\.External Links:ISSN 0266\-3538,[Document](https://dx.doi.org/10.1016/j.compscitech.2014.10.009)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p1.1)\.
- \[68\]S\. Torquato\(2002\)Random Heterogeneous Materials: Microstructure and Macroscopic Properties\.Vol\.55,Springer Nature\.Note:Publication Title: Applied Mechanics Reviews \- APPL MECH REVExternal Links:[Document](https://dx.doi.org/10.1115/1.1483342)Cited by:[§3\.2](https://arxiv.org/html/2606.28996#S3.SS2.p6.1)\.
- \[69\]K\. Uhlig, M\. Tosch, L\. Bittrich, A\. Leipprand, S\. Dey, A\. Spickenheuer, and G\. Heinrich\(2016\)Meso\-scaled finite element analysis of fiber reinforced plastics made by Tailored Fiber Placement\.Composite Structures143,pp\. 53–62\.External Links:ISSN 0263\-8223,[Document](https://dx.doi.org/10.1016/j.compstruct.2016.01.049)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p4.1)\.
- \[70\]Y\. Vitulyova, T\. Babenko, K\. Kolesnikova, N\. Kiktev, and O\. Abramkina\(2025\)A Hybrid Approach Using Graph Neural Networks and LSTM for Attack Vector Reconstruction\.Computers14,pp\. 301\(en\)\.External Links:ISSN 2073\-431X,[Document](https://dx.doi.org/10.3390/computers14080301)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p7.1)\.
- \[71\]L\. Wang, Z\. Wang, S\.M\. Dong, W\. Zhang, and Y\. Wang\(2013\)Finite element simulation of stress distribution and development of cf/sic ceramic–matrix composite coated with single layer sic coating during thermal shock\.Composites Part B: Engineering51,pp\. 204–214\.External Links:[Document](https://dx.doi.org/10.1016/j.compositesb.2013.03.028)Cited by:[1st item](https://arxiv.org/html/2606.28996#S2.I1.i1.p1.1)\.
- \[72\]L\. Wu, L\. Adam, and L\. Noels\(2021\)Micro\-mechanics and data\-driven based reduced order models for multi\-scale analyses of woven composites\.Composite Structures270,pp\. 114058\.External Links:ISSN 0263\-8223,[Document](https://dx.doi.org/10.1016/j.compstruct.2021.114058)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p6.1)\.
- \[73\]H\. Xu, M\. Kuczynska, N\. Schafet, F\. Welschinger, and J\. Hohe\(2022\)FE\-based damage modeling approach for short fiber reinforced thermoplastics under quasi\-static load coupling anisotropic viscoplasticity and matrix degradation\.Journal of Composite Materials56,pp\. 3113–3125\.External Links:[Document](https://dx.doi.org/10.1177/00219983221109327)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p4.1)\.
- \[74\]M\. Yacouti and M\. Shakiba\(2025\)Integrated convolutional and graph neural networks for predicting mechanical fields in composite microstructures\.Composites Part A: Applied Science and Manufacturing190,pp\. 108618\(en\)\.External Links:ISSN 1359835X,[Document](https://dx.doi.org/10.1016/j.compositesa.2024.108618)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p7.1)\.
- \[75\]S\. Yu, J\. Y\. Hwang, and S\. H\. Hong\(2020\)3D microstructural characterization and mechanical properties determination of short basalt fiber\-reinforced polyamide 6,6 composites\.Composites Part B: Engineering187,pp\. 107839\(en\)\.External Links:ISSN 13598368,[Document](https://dx.doi.org/10.1016/j.compositesb.2020.107839)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p2.1)\.
- \[76\]Y\. Zhao, Z\. Chen, and X\. Jian\(2023\)A high\-generalizability machine learning framework for analyzing the homogenized properties of short fiber\-reinforced polymer composites\.Polymers15,pp\. 3962\.External Links:[Document](https://dx.doi.org/10.3390/polym15193962)Cited by:[§1](https://arxiv.org/html/2606.28996#S1.p7.1)\.

## Appendix A1\. Parametric Finite\-Element Validation of the Edge Weight Function

The edge weight function \(Equation[9](https://arxiv.org/html/2606.28996#S2.E9)\) embeds two physically motivated assumptions regarding pairwise fiber–fiber load transfer: \(i\) matrix\-mediated stress transfer between proximal fiber ends decays with increasing fiber–fiber separation, and \(ii\) load\-transfer efficiency depends on relative fiber orientation through a dominant cosine\-dependent axial shear\-lag coupling and a residual orientation\-independent transverse stress concentration\. To confirm the physical basis of both terms in the CF\-ABS system at geometries consistent withμ\\mu\-CT observations, a systematic parametric finite\-element study was conducted using a two\-fiber RUC model\. Because Equation[9](https://arxiv.org/html/2606.28996#S2.E9)governs physics at the pairwise level \(multi\-fiber effects are handled through GNN message\-passing aggregation\), the two\-fiber RUC represents the minimal and most direct unit for validating its functional form\.

##### Two\-fiber RUC model geometry

The model consists of two cylindrical carbon fibers of equal lengthLf​1=Lf​2=50​μL\_\{f1\}=L\_\{f2\}=50~\\mum and radiusrf=4​μr\_\{f\}=4~\\mum embedded in an ABS matrix cuboid \(Figure[17](https://arxiv.org/html/2606.28996#Sx1.F17)\)\. The first fiber is fixed in all analyses; the second is varied such that the end\-to\-end separationLg​a​pL\_\{gap\}ranges from 20 to 70μ\\mum in six increments, and the relative orientation angleϕ\\phitakes values of0∘0^\{\\circ\},30∘30^\{\\circ\},60∘60^\{\\circ\}, and90∘90^\{\\circ\}, yielding 24 parametric cases\. All cases were analyzed under uniaxial tensile strainε0=4%\\varepsilon\_\{0\}=4\\%applied along the primary fiber longitudinal direction, consistent with the homogenization protocol of Section 2\.3\.3\. For angled configurations, RUC dimensions were adjusted to accommodate the tilted geometry while keeping fiber dimensions fixed, resulting in orientation\-dependent fiber volume fractions:Vf=7\.86%V\_\{f\}=7\.86\\%at0∘0^\{\\circ\},6\.60%6\.60\\%at30∘30^\{\\circ\},5\.60%5\.60\\%at60∘60^\{\\circ\}, and4\.83%4\.83\\%at90∘90^\{\\circ\}\.

![Refer to caption](https://arxiv.org/html/2606.28996v1/FigureA2_1.png)Figure 17:Two\-fiber RUC model of CF\-ABS used for parametric shear\-lag studies\.Left:ϕ=0∘\\phi=0^\{\\circ\}with fiber lengthsLf​1,Lf​2L\_\{f1\},L\_\{f2\}and end\-to\-end separationLgapL\_\{\\text\{gap\}\}labeled\.Right:ϕ=30∘\\phi=30^\{\\circ\}case\.
##### Interfacial shear stress extraction

For each parametric case, the interfacial shear stress resultant:

τ=τx​y2\+τy​z2\+τz​x2\\tau=\\sqrt\{\\tau\_\{xy\}^\{2\}\+\\tau\_\{yz\}^\{2\}\+\\tau\_\{zx\}^\{2\}\}\(23\)was extracted at the first matrix integration point immediately adjacent to the fiber–matrix interface along the full fiber axis fromz=0z=0toz=Lfz=L\_\{f\}\. The peak interfacial shear stressτp​e​a​k\\tau\_\{peak\}was identified and normalized by localVfV\_\{f\}to enable consistent comparison across orientations \(Figure[18](https://arxiv.org/html/2606.28996#Sx1.F18)b\)\.

##### Validation of the distance\-decay term

![Refer to caption](https://arxiv.org/html/2606.28996v1/FigureA2_2.png)Figure 18:Shear\-lag response of two\-fiber CF\-ABS RVEs across all 24 parametric cases \(four orientations×\\timessix gap distances\)\. \(a\) Absolute interfacial shear stressτ​\(z\)\\tau\(z\)profiles extracted along the matrix–fiber interface path\. Each set of curves corresponds to one fiber orientation \(ϕ=0∘\\phi=0^\{\\circ\},30∘30^\{\\circ\},60∘60^\{\\circ\},90∘90^\{\\circ\}\), with individual curves representing gap distancesLgap=20L\_\{\\mathrm\{gap\}\}=20–7070mm\. \(b\)VfV\_\{f\}\-normalised peak interfacial shear stressτpeak/Vf\\tau\_\{\\mathrm\{peak\}\}/V\_\{f\}\[MPa per%Vf\\%V\_\{f\}\] as a three\-dimensional surface over theLgapL\_\{\\mathrm\{gap\}\}–ϕ\\phiparameter space\. Normalization by local fiber volume fraction removes the geometricVfV\_\{f\}variation introduced by the tilted\-fiber box geometry\. \(c\) Fiber length efficiencyηl\\eta\_\{l\}estimated from the data\-driven exponential decay of the normalisedτ\\tauprofile, shown as a three\-dimensional surface over theLgapL\_\{\\mathrm\{gap\}\}–ϕ\\phiparameter space\.Figure[18](https://arxiv.org/html/2606.28996#Sx1.F18)a presentsτ\\tauprofiles for all 24 cases\. Across all four orientations, interfacial shear concentrated at both fiber termini, confirming that end\-to\-end proximity governs pairwise load transfer regardless of orientation\. At0∘0^\{\\circ\}and30∘30^\{\\circ\},τp​e​a​k\\tau\_\{peak\}decreased monotonically with increasingLg​a​pL\_\{gap\}\(15–22% reduction at0∘0^\{\\circ\}between 20 and 70μ\\mum\), directly confirming the exponential decay termexp⁡\(−dp​q/Lc\)\\exp\(\-d\_\{pq\}/L\_\{c\}\)in Equation[9](https://arxiv.org/html/2606.28996#S2.E9)\. At60∘60^\{\\circ\}and90∘90^\{\\circ\}, the gap effect was negligible \(less than 8%\), confirming that the distance\-decay mechanism is physically operative only when fibers participate in axial shear\-lag\-governed load transfer—precisely the condition captured by the orientation\-dependent component\.

##### Validation of the two\-component orientation factor

Figure[18](https://arxiv.org/html/2606.28996#Sx1.F18)b reveals two physically distinct orientation regimes motivating the\(0\.9​\|𝐭p⋅𝐭q\|\+0\.10\)\(0\.9\|\\mathbf\{t\}\_\{p\}\\cdot\\mathbf\{t\}\_\{q\}\|\+0\.10\)structure in Equation 9\. TheVfV\_\{f\}\-normalized peak stress decreased systematically with misalignment, from 25–32 MPa per %VfV\_\{f\}at0∘0^\{\\circ\}to 2–7 MPa per %VfV\_\{f\}at90∘90^\{\\circ\}—an 8–10\-fold reduction that validates the dominant cosine\-dependent axial shear\-lag term\. Critically, even at90∘90^\{\\circ\}a non\-negligible residual of 2–7 MPa per %VfV\_\{f\}persisted, arising from transverse stress concentration due to fiber–matrix stiffness mismatch rather than axial shear\-lag\. This residual was gap\-insensitive \(less than 8% variation across the fullLg​a​pL\_\{gap\}range\), confirming its orientation\-independent character\. The observed residual fraction of 0\.08–0\.22 relative to the0∘0^\{\\circ\}maximum directly yields the constant 0\.10 in Equation[9](https://arxiv.org/html/2606.28996#S2.E9), with the complementary 0\.9 assigned to the cosine\-dependent term\.

##### Fiber length efficiency

The fiber length efficiencyηl\\eta\_\{l\}, presented as a three\-dimensional surface in Figure[18](https://arxiv.org/html/2606.28996#Sx1.F18)c, quantifies the fraction of theoretical fiber stiffness contribution activated over the actual fiber length as directly observed from the two\-fiber FEA results\. Values ofηl=0\.72\\eta\_\{l\}=0\.72–0\.750\.75for0∘0^\{\\circ\}fiber pairs confirm that, at the fiber lengths and separations present in theμ\\mu\-CT microstructure, pairwise fiber–fiber interaction is strongly active and a large fraction of available reinforcement efficiency is realized through end\-to\-end load transfer\. For30∘30^\{\\circ\}fibers,ηl=0\.20\\eta\_\{l\}=0\.20–0\.380\.38, consistent with the intermediate\-orientation behavior described above\.

## Appendix A2\. GNN\-LSTM Training Configuration

Table[4](https://arxiv.org/html/2606.28996#Sx2.T4)consolidates the hyperparameters used for GNN\-LSTM optimization \(described in Section[2\.4](https://arxiv.org/html/2606.28996#S2.SS4)\)\.

Table 4:Training hyperparameters used for GNN–LSTM model optimization\.Two implementation choices warrant brief justification\. A batch size of 1 was used because Voronoi cells contain variable numbers of fibers, which would otherwise require variable\-size padding for batched graph processing; sequential processing avoids this overhead while gradient norm clipping maintains stability\. The 500\-cell dataset was partitioned using stratified random shuffling rather than a sequential split, since a sorted split would risk grouping all low\-fiber\-count or high\-VfV\_\{f\}cells into a single subset, producing an out\-of\-distribution validation set and overoptimistic generalization estimates\.

## Appendix A3\. Full Six\-Mode Constitutive Response of the Discrete Cell

To verify that the discrete\-cell FE homogenization produces physically consistent constitutive responses under all six independent loading modes, Figure[19](https://arxiv.org/html/2606.28996#Sx3.F19)presents the full von Mises stress fields for a representative Voronoi cell subjected to three normal and three shear loading cases\. Under normal loading in theyy\- andzz\-directions \(Figure[19](https://arxiv.org/html/2606.28996#Sx3.F19)\(a–b\), stress is distributed relatively uniformly through the matrix, consistent with matrix\-dominated transverse response and the comparatively small contribution of off\-axis fibers to transverse stiffness\. Under normal loading in thexx\-direction \(fiber\-dominated, Figure[19](https://arxiv.org/html/2606.28996#Sx3.F19)c\), pronounced stress localization occurs along fiber axes, confirming efficient axial load transfer\.

Under shear loading \(Figure[19](https://arxiv.org/html/2606.28996#Sx3.F19)\(d–f\)\), stress concentrations developed in matrix ligaments between adjacent fibers, consistent with shear\-lag\-governed inter\-fiber coupling\. These include one transverse shear configuration acting within the cross\-sectionalzz–yyplane, and two axial shear configurations acting across the longitudinalyy–zzandxx–zzfiber planes\. These results confirm that the discrete\-cell FE implementation correctly captures the anisotropic, microstructure\-sensitive constitutive behavior that the GNN\-LSTM surrogate is trained to reproduce\.

![Refer to caption](https://arxiv.org/html/2606.28996v1/FigureA4_1.png)Figure 19:Full six\-mode constitutive response of a representative discrete Voronoi cell from finite\-element homogenization\. Von Mises stress fields are shown under: \(a\) normal loading in thexx\-direction \(σx​x\\sigma\_\{xx\}\); \(b\) transverse loading in theyy\-direction \(σy​y\\sigma\_\{yy\}\); \(c\) transverse loading in thezz\-direction \(σ33\\sigma\_\{33\}\); \(d\) axial shear in thexx–yyplane \(τ12\\tau\_\{12\}\); \(e\) axial shear in theyy–zzplane \(τ23\\tau\_\{23\}\); \(f\) axial shear in thexx–zzplane \(τ13\\tau\_\{13\}\)\. The contrasting stress patterns under normal versus shear loading confirm that the cell’s constitutive response is strongly anisotropic and topology\-sensitive\.
## Supplemental Material\-A\. Sphere\-Sliding Visualization