A Neural Network Transformer Model for Composite Microstructure
  Homogenization
Heterogeneity and uncertainty in a composite microstructure lead to either
computational bottlenecks if modeled rigorously or to solution inaccuracies in
the stress field and failure predictions if approximated. Although methods
suitable for analyzing arbitrary and non-linear microstructures exist, their
computational cost makes them impractical to use in large-scale structural
analysis. Surrogate models or Reduced Order Models (ROMs) commonly enhance
efficiencies but are typically calibrated with a single microstructure.
Homogenization methods, such as the Mori-Tanaka method, offer rapid
homogenization for a wide range of constituent properties. However, simplifying
assumptions, like stress and strain averaging in phases, render the
consideration of both deterministic and stochastic variations in microstructure
infeasible. This paper illustrates a transformer neural network architecture
that captures the knowledge of various microstructures and constituents,
enabling it to function as a computationally efficient homogenization surrogate
model. Given an image or an abstraction of an arbitrary composite
microstructure of linearly elastic fibers in an elastoplastic matrix, the
transformer network predicts the history-dependent, non-linear, and homogenized
stress-strain response. Two methods for encoding microstructure features were
tested: calculating two-point statistics using Principal Component Analysis
(PCA) for dimensionality reduction and employing an autoencoder with a
Convolutional Neural Network (CNN). Both methods accurately predict the
homogenized material response. The developed transformer neural network offers
an efficient means for microstructure-to-property translation, generalizable
and extendable to a variety of microstructures. The paper describes the network
architecture, training and testing data generation, and performance under
cycling and random loadings.
Emil Pitz, Kishore Pochiraju
# A Neural Network Transformer Model for Composite Microstructure Homogenization

###### Abstract

Heterogeneity and uncertainty in a composite microstructure lead to either computational bottlenecks if modeled rigorously, or to solution inaccuracies in the stress field and failure predictions if approximated. Although methods suitable for analyzing arbitrary and non-linear microstructures exist, their computational cost makes them impractical to use in large-scale structural analysis. Surrogate models or Reduced Order Models (ROMs), commonly enhance efficiencies, but they are typically calibrated with a single microstructure. Homogenization methods, such as the Mori-Tanaka method, offer rapid homogenization for a wide range of constituent properties. However, simplifying assumptions, like stress and strain averaging in phases, render the consideration of both deterministic and stochastic variations in microstructure infeasible.

This paper illustrates a transformer neural network architecture that captures the knowledge of various microstructures and constituents, enabling it to function as a computationally efficient homogenization surrogate model. Given an image or an abstraction of an arbitrary composite microstructure of linearly elastic fibers in an elastoplastic matrix, the transformer network predicts the history-dependent, non-linear, and homogenized stress-strain response. Two methods were tested that encode features of the microstructure. The first method calculates two-point statistics of the microstructure and uses Principal Component Analysis (PCA) for dimensionality reduction. The second method uses an autoencoder with a Convolutional Neural Network (CNN). Both microstructure encoding methods accurately predict the homogenized material response. The developed transformer neural network offers an efficient means for non-linear and history-dependent microstructure-to-property translation, which is generalizable and extendable to a wide variety of microstructures. The paper describes the network architecture, training and testing data generation and the performance of the transformer network under cycling and random loadings.

Transformer Microstructure Homogenization Surrogate Model History-dependence Microstructure encoding

## 1 Introduction

Multiple studies have shown that Neural Networks (NNs) can replace constitutive laws of homogeneous materials for various material behaviors, such as plastic behavior [1, 2], plastic material under cyclic loading [3], hyperelastic materials [4], or for modeling traction-separation behavior [5]. Extensions of NN models that predict the orthotropic linear elastic properties of 2D and 3D-microcells have been shown in [6, 7, 8] using CNNs. A surrogate model for the anisotropic electrical response of graphene/polymer nanocomposites was presented in [9].

To model the non-linear response of heterogeneous materials,NNs surrogates of a non-linear strain energy function were presented in [10, 11]. In [12, 13], localized Mises stress distributions were generated using CNNs for a given
monotonic uniaxial tensile loading at a prescribed loading rate. The approaches in [12, 13] were limited to monotonic loading and history dependence was ignored. The loading history dependence can be significant for the homogenized material level to find the stress state for irreversible behavior under general loading conditions [14]. Therefore, in [15] state variables were used to account for the loading history and build a surrogate for the yield function of foams. An alternative to using state variables is a NNs that takes into account of entire loading history. An Long Short-Term Memory (LSTM) based Recurrent Neural Network (RNN) that considers strain history and microstructure information of elastoplastic composites was described in [16]. Similarly, a RNN surrogate was implemented in [17] for history-dependent homogenization to accelerate nested (micro-macro) finite element (FE\({}^{2}\)) analysis with automatic differentiation of the surrogate model to obtain consistent tangent matrices. Wu, Nguyen, Kilingar, and Noels used a RNN with Gated Recurrent Units (GRUs) to implement a surrogate model for elastoplastic composite microcells in [14]. The methodology was extended to include localization in the surrogate model to extract the micro stress field using dimensionality reduction with PCA [18].

In this paper, we described a surrogate model that can predict the response of arbitrary microstructures with two constituent phases. The model can predict the current (history-dependent, e.g., due to plasticity) stress state average over the Representative Volume Element (RVE) \(\overline{\mathbf{\sigma}}\) at time \(t\) as a result of the strain history \(\overline{\mathbf{c}}_{1:t}\) for a given microstructure defined by parameters \(\mathbf{m}\) and constituent properties \(\mathbf{p}\)[16]:

\[\overline{\mathbf{\sigma}}=f\left(\overline{\mathbf{c}}_{1:t},\mathbf{m},\mathbf{p},t\right). \tag{1}\]

The transformer network architecture was developed by Vaswani _et al._ in [19] was used as an alternative to RNNs. Transformer NNs have been shown to have impressive performance in language modeling tasks, machine translation and a variety of other applications (see Sec. 2). To date, transformers have not yet been used as surrogate models for a history-dependent material response.

Transformers in the context of Natural Language Processing (NLP) applications when trained on diverse data have shown an impressive ability to generalize. A single model has been able to generate text, and write and debug code in different programming languages [20]. We believe that the ability of transformers to generalize can be extended to the field of microstructure homogenization. We hypothesized that a transformer could be trained on a variety of material combinations and microstructures and function as a fast and efficient homogenization tool. Because of the prediction speed, the trained model could be used to perform Uncertainty Quantification (UQ) of the response of composite structures with complex and varying microstructure, where current homogenization techniques still are a computational bottleneck.

In this effort, the feasibility of using a transformer network to predict the history-dependent material response was first trained and validated with homogeneous materials. Following the successful prediction of the history-dependent homogeneous material response, the surrogate model was extended and trained for heterogeneous materials. Furthermore, in NLP, pre-training of transformers, where the network is first trained on large unlabeled datasets followed by task-specific fine-tuning, has shown major performance advantages [21]. Following a similar strategy, we used the weights for homogeneous data as a starting point for training with heterogeneous materials. This strategy reduced the amount of computationally-expensive microstructure homogenization data required for the teach homogenization to the transformer network.

This paper will first introduce the neural network techniques required to solve the homogenization problem and give an overview of similar approaches in the literature. The process to generate the required training data and represent the microstructure information will then be described, followed by a description of the used network architecture and the training process. Finally, the results will be presented.

## 2 Applications and Architecture of the Transformer

Transformer models are recent compared to RNNs to analyze sequential data. First presented in [19], transformers have been successful in NLP applications and more recently have been applied to a wider range of problems. Vaswani _et al._ already showed impressive performance gains in language translation tasks for the original transformer [19]. The recent public release of the transformer-based Large Language Model (LLM) ChatGPT by OpenAI has attracted major media attention [22]. ChatGPT is a version of the Generative Pre-trained Transformer (GPT) models developed by OpenAI [22, 23]. Successively larger models (GPT: 110 million parameters, GPT-2: 1.5 billion parameters, GPT-3: 175 billion parameters [24]) have led to an impressive ability to generate human-like text. The new GPT-4 reportedly has 100 trillion parameters and works with both text and images [25]. The GPT models use a decoder-only architecture with the decoder architecture relatively unchanged from the original transformer model [19, 22, 23]. The GPT models are trained on unlabeled data, allowing the models to be trained on data that require virtually no pre-processing, alleviating the need for manually labeled data.
Aside from NLP, multifaceted successful applications of transformers have been shown recently. Examples include "AlphaFold", a neural network incorporating transformers for predicting a protein's three-dimensional structure based on the amino acid sequence with atomic accuracy [26, 27, 28]. Applications have been shown in computer vision such as for image classification [29, 30] and a wide array of additional tasks such as object detection, image segmentation, or pose estimation [31]. In [32], a transformer was adapted to have as input three-dimensional meshes and perform shape classification and segmentation tasks. Two recent publications use transformers to improve the estimation of counterfactual outcomes in, e.g., healthcare settings [33, 34]. Furthermore, studies have shown that transformers or transformer-based models can be used for time series forecasting based on previously observed data, such as predicting the prevalence of influence [35, 36, 37]. In [38], a transformer-based time series forecasting model named Temporal Fusion Transformer was proposed for multi-horizon forecasting, i.e., prediction of multiple variables of interest. The model incorporates not only previous time series data but a mix of static and time-dependent inputs, and known future inputs.

Following the architecture of most sequence-to-sequence models, the transformer has an encoder-decoder structure [19]. Given an input sequence \(\mathbf{x}=(x_{1},\ldots,x_{n})\), the encoder generates an intermediate sequence \(\mathbf{z}=(z_{1},\ldots,z_{n})\). The decoder uses \(\mathbf{z}\) to element-wisely generate the output \(\mathbf{y}=(y_{1},\ldots,y_{n})\). Previously generated output values are auto-regressively used as supplementary input [19]. Commonly, models to process data with sequential dependencies use recurrent layers to process individual values sequentially [39]. Transformers, in contrast, solely rely on so-called attention mechanisms in the encoder and decoder to represent dependencies between input and output sequences [19, 40]. Attention was initially introduced in [40] to improve the performance of sequence models for long-range dependencies. Vaswani _et al._ used a scaled dot-product attention mechanism given by [19]:

\[\text{Attention}\left(\mathbf{Q},\mathbf{K},\mathbf{V}\right)=\text{softmax}\left( \frac{\mathbf{Q}\mathbf{K}^{T}}{\sqrt{d_{k}}}\right)\mathbf{V} \tag{2}\]

where \(\mathbf{Q}\), \(\mathbf{K}\), and \(\mathbf{V}\) are vectors of dimension \(d_{k}\), \(d_{k}\), and \(d_{v}\), named query, keys, and values, respectively, and softmax is the softmax activation function. In the transformer, instead of using the attention mechanism a single time, Vaswani _et al._ use so-called multi-head attention. The attention function is performed in parallel \(h\) times, where \(h\) is a selected parameter (eight in [19]) and represents the number of attention heads. Additionally, \(\mathbf{Q}\), \(\mathbf{K}\), and \(\mathbf{V}\) are linearly projected to \(d_{v}\), \(d_{k}\), and \(d_{k}\) dimensions, respectively, using projection matrices \(\mathbf{W}\) that are learned during model training [19]:

\[\text{head}_{i}=\text{Attention}\left(\mathbf{Q}\mathbf{W}_{i}^{Q},\mathbf{K}\mathbf{W}_{i}^ {K},\mathbf{V}\mathbf{W}_{i}^{V}\right)\qquad\text{with}\quad i=1,\ldots,h. \tag{3}\]

This operation produces \(d_{v}\)-dimensional output for each attention head. The output of the \(h\) attention heads is then concatenated and linearly projected again using the learned projection matrix \(\mathbf{W}^{O}\) to produce the multi-head attention output [19]:

\[\text{MultiHead}\left(\mathbf{Q},\mathbf{K},\mathbf{V}\right)=\text{Concat}\left(\text{ head}_{1},\ldots,\text{head}_{h}\right)\mathbf{W}^{O}. \tag{4}\]

From this, the dimensions of the projection matrices follow [19]:

\[\begin{split}\mathbf{W}_{i}^{Q}&\in\mathbb{R}^{d_{ model}\times d_{k}}\\ \mathbf{W}_{i}^{K}&\in\mathbb{R}^{d_{model}\times d_{k}} \\ \mathbf{W}_{i}^{V}&\in\mathbb{R}^{d_{model}\times d_{k}} \\ \mathbf{W}^{O}&\in\mathbb{R}^{hd_{v}\times d_{model}}, \end{split} \tag{5}\]

with the selected model size \(d_{model}\) (e.g., \(d_{model}=512\) in [19]).

There are contradicting studies on whether attention can provide reasoning for a model's predictions by laying out the relative importance of inputs [41, 42]. In [41] the authors compared attention weights with feature importance and found a low correlation, arguing that attention provides little transparency. In contrast, in [42] the authors argue that attention can provide a plausible explanation for a model's predictions, however, not _"one true, faithful interpretation of the link...between inputs and outputs"_.

Using the multi-head attention mechanism, the architecture of the transformer is depicted in Fig. 1. The encoder layer consists of a multi-head attention mechanism and a simple fully connected feed-forward layer, both having residual connections [19, 43]. After both attention and feed-forward layers, normalization is performed. The encoder layer is repeated \(N\) times (\(N=6\) in [19]). Keys, values, and queries of the multi-head attention mechanism are the embedded model input (first encoder layer) or the output of the previous encoder layers, making the attention mechanism a self-attention mechanism.

A similar architecture is used for the decoder layers, with an additional multi-head attention mechanism ("encoder-decoder attention", middle sub-layer in Fig. 1) that uses as values and keys the output of the encoder and as queries, the
output of the previous decoder layers [19]. The encoder-decoder multi-head attention mechanism is similarly followed by layer normalization and includes residual connections. Furthermore, the self-attention mechanism attending to the embedded model output (first decoder sub-layer in Fig. 1) is masked and the output sequence is shifted to the right by one position. Masked self-attention effectively prevents the decoder from "looking into the future", i.e., the decoder predictions can only depend on known past output positions. Again, the decoder layer is repeated \(N\) times. Output from the decoder is finally passed through a linear layer followed by softmax activation to obtain the output probabilities.

The inputs and outputs of the model have to be converted to vectors of dimension \(d_{model}\)[19]. In NLP applications, commonly learned embeddings are used for this purpose [44].

Figure 1: Architecture of the transformer model. (adapted from [19])
Unlike RNNs, the transformer does not use any recurrence and, therefore, does not directly have information about the location of a value in a sequence [19]. Instead, positional information is injected by adding a positional encoding to the embedded input and output. The positional encoding consists of sine and cosine functions of different frequencies.

Additionally, transformers can be advantageous in multiple aspects compared to RNNs: Long-range dependencies in time sequences are problematic due to the recurrent nature of RNNs [19, 36]. Transformers seem to be better able to capture long-range dependencies because self-attention mechanisms connect all positions in a sequence with the same number of operations, compared to sequentially in RNNs [19]. Furthermore, RNNs can suffer from exploding or vanishing gradients during training, making training the model more difficult [39].

## 3 Generation of Training data

The objective of this effort was to develop a transformer neural network surrogate to predict the non-linear and history-dependent response of composite microcells. Generation of the training dataset requires analysis of a microcell's response to deformation through scale transition from the micro-scale to the macro-scale. An introduction to Finite Element Method (FEM)-based homogenization to solve the multi-scale Boundary Value Problem (BVP) and generate the stress response of a microcell for finite-strain problems will be given in the following section (Sec. 3.1).

### Solving the Multi-Scale BVP to Calculate the Microcell Stress Response

Given a (homogenized) macro-scale deformation gradient tensor \(\mathbf{F}_{M}\), the homogenized Piola-Kirchhoff stress tensor of the underlying microcell \(\mathbf{P}_{M}\left(\mathbf{X},t\right)\) at (macro) material point \(\mathbf{X}\) and time \(t\) shall be found (see Fig. 2). The homogenized stress tensor \(\mathbf{P}_{M}\) can be found by solving the micro-scale BVP using FEM. In the following equations, \(M\) represents the macro-scale values and \(m\) refers to values on the micro-scale. The given summary of the multi-scale homogenization problem is based on the derivation presented in [14, 18].

The linear momentum equations at the macro-scale in the domain \(\Omega\) and neglecting dynamic effects are given by [14, 18]:

\[\mathbf{P}_{M}\left(\mathbf{X}\right)\cdot\nabla_{0}+\mathbf{b}_{M}=\mathbf{0}\qquad\forall \mathbf{X}\in\Omega, \tag{6}\]

with the gradient operator with regard to the reference configuration \(\nabla_{0}\) and the body force \(\mathbf{b}_{M}\). The linear momentum equations are subject to the Dirichlet boundary conditions on the boundary \(\partial_{D}\Omega\) and Neumann boundary conditions on the boundary \(\partial_{N}\Omega\):

\[\mathbf{u}_{M}\left(\mathbf{X}\right)=\hat{\mathbf{u}}_{M}\qquad\forall\mathbf{X} \in\partial_{D}\Omega \tag{7}\] \[\mathbf{P}_{M}\left(\mathbf{X}\right)\cdot\mathbf{N}_{M}=\hat{\mathbf{T}}_{M}\qquad \forall\mathbf{X}\in\partial_{N}\Omega, \tag{8}\]

with the displacements \(\mathbf{u}_{M}\), the constrained displacements \(\hat{\mathbf{u}}_{M}\) on \(\partial_{D}\Omega\), the unit normal \(\mathbf{N}_{M}\) on \(\partial_{N}\Omega\), and the surface traction \(\hat{\mathbf{T}}_{M}\).

Figure 2: Scale-transition: Given a macro-scale deformation gradient \(\mathbf{F}_{M}\) the corresponding homogenized Piola-Kirchhoff stress \(\mathbf{P}_{M}\) is found by solving the micro-scale BVP. (adapted from [14])
On the micro-scale, the BVP is commonly defined on a parallelepipedic RVE with the boundaries of the RVE being planar faces \(\partial\omega\). The linear momentum equations on the micro-scale are given by:

\[\mathbf{P}_{m}\left(\mathbf{x}\right)\cdot\nabla_{0}=\mathbf{0}\qquad\forall\mathbf{x}\in\omega, \tag{9}\]

subject to boundary conditions:

\[\mathbf{P}_{m}\left(\mathbf{x}\right)\cdot\mathbf{N}_{m}=\hat{\mathbf{T}}_{m}\qquad\forall\mathbf{ x}\in\partial\omega, \tag{10}\]

where the traction \(\mathbf{T}_{m}\) acts on the boundary \(\partial\omega\) with unit normal \(\mathbf{N}_{m}\).

The macro scale deformation gradient tensor \(\mathbf{F}_{M}\) and stress tensor \(\mathbf{P}_{M}\) are given by volume average over \(\omega\) of the micro scale deformation gradient \(\mathbf{F}_{m}\) and micro scale stress tensor \(\mathbf{P}_{m}\), respectively [14]:

\[\mathbf{F}_{M}\left(\mathbf{X},t\right) =\frac{1}{V\left(\omega\right)}\int_{\omega}\mathbf{F}_{m}\left(\mathbf{ x},t\right)d\mathbf{x} \tag{11}\] \[\mathbf{P}_{M}\left(\mathbf{X},t\right) =\frac{1}{V\left(\omega\right)}\int_{\omega}\mathbf{P}_{m}\left(\mathbf{ x},t\right)d\mathbf{x} \tag{12}\]

Energy consistency requires the Hill-Mandel condition to be fulfilled [14, 18, 45]:

\[\mathbf{P}_{M}:\delta\mathbf{F}_{M}=\frac{1}{V\left(\omega\right)}\int_{\omega}\mathbf{ P}_{m}:\delta\mathbf{F}_{m}d\mathbf{x}. \tag{13}\]

On the micro-scale, the displacement field \(\mathbf{u}_{m}\) with an introduced pertubation field \(\mathbf{u}^{\prime}\left(x\right)\), the macro-scale displacement gradient \(\mathbf{F}_{m}\), and unit tensor \(\mathbf{I}\) is defined by [14]:

\[\mathbf{u}_{m}\left(\mathbf{x}\right)=\left(\mathbf{F}_{M}-\mathbf{I}\right)\cdot\left(\mathbf{x }-\mathbf{x}_{0}\right)+\mathbf{u}^{\prime}\left(\mathbf{x}\right). \tag{14}\]

The point \(\mathbf{x}_{0}\) is a reference point in \(\omega\). Using trial and test functions, the weak form of the micro-scale BVP can be formulated:

\[\int_{\omega}\mathbf{P}_{m}\left(\mathbf{u}^{\prime}\right):\left(\partial\mathbf{u}^{ \prime}\otimes\nabla_{0}\right)d\mathbf{x}=0,\qquad\forall\partial\mathbf{u}^{\prime }\in\mathcal{U}\left(\omega\right)\subset\mathcal{U}^{\text{min}}\left(\omega \right), \tag{15}\]

where \(\partial\mathbf{u}^{\prime}\in\mathcal{U}\left(\omega\right)\) is a test function in the admissible kinematic vector field \(\mathcal{U}\left(\omega\right)\), \(\mathbf{u}^{\prime}\in\mathcal{U}\left(\omega\right)\subset\mathcal{U}^{\text{ min}}\left(\omega\right)\) shall be found, and \(\mathcal{U}^{\text{min}}\left(\omega\right)\) is the minimum kinematic vector field.

Given a prescribed macro-scale deformation gradient \(\mathbf{F}_{M}\), the weak form can be solved using FEM to find the homogenized stress tensor \(\mathbf{P}_{M}\). On opposite faces of the RVE \(\mathbf{x}^{+}\in\omega^{+}\) and \(\forall\mathbf{x}^{-}\in\omega^{-}\), Periodic Boundary Conditions (PBCs) are applied:

\[\begin{split}\mathcal{U}^{\text{BC}}\left(\omega\right)=& \left\{\mathbf{u}^{\prime}\in\mathcal{H}\left(\omega\right)\left|\mathbf{u}^{\prime} \left(\mathbf{x}^{+}\right)=\mathbf{u}^{\prime}\left(\mathbf{x}^{-}\right)\right.,\\ &\left.\forall\mathbf{x}^{+}\in\omega^{+}\right.\quad\text{and corresponding}\quad\forall\mathbf{x}^{-}\in\omega^{-}\right\}\subset\mathcal{U}^{\text{min}}\left(\omega \right),\end{split} \tag{16}\]

with the Hilbert space \(\mathcal{H}\). While various kinds of boundary conditions can be applied to fulfill the Hill-Mandel condition [14], PBCs were chosen for this work for ease of implementation.

### Description of the Micro-Scale Composite Volume Elements

To generate the training dataset, a workflow was implemented to:

1. Generate random composite Stochastic Volume Elements (SVEs).
2. Setup and execute FEM simulations with PBCs and given homogenized load to generate the homogenized microcells' response by solving the micro-scale BVPs.
3. Save the microstructure data and homogenized loading and response (as sequences of strain increments and corresponding stress) in a database to use as training data.

The term SVE refers to a volume element of the analyzed composite material that shows varying homogenized properties depending on the SVE realization [46] (e.g., given multiple volume elements with similar Fiber Volume Ratio (FVR), each volume element will yield different results). This is in contrast to an RVE, which assumes statistical representativity.

For simplicity, in this effort, the application was restricted to two-dimensional problems, more specifically to the transversal behavior of continuously reinforced composites assuming plane strain conditions. Therefore, the generated volume elements consisted of circular fibers inside a matrix (see Fig. 3 for three arbitrary SVEs at different volume
fractions). Given plane-strain conditions and with the \(z\)-axis being the out-of-plane direction (normal to an analyzed SVE), the three non-zero strain components \(\varepsilon_{xx}\), \(\varepsilon_{yy}\), and \(\varepsilon_{xy}\) can be prescribed resulting in four stress components \(\sigma_{xx}\), \(\sigma_{yy}\), \(\sigma_{zz}\), and \(\sigma_{xy}\).

Periodic SVEs were generated by randomly placing fibers inside the microcell until a desired FVR was reached. Uniform distributions were used to place the fibers. The uniformly distributed fiber locations inside the generated microstructures may not represent the fiber location distributions found in a composite material [46]. Further analysis of the fiber location distributions could be an interesting objective for future model improvements. After generating the microcell geometry, the SVE was meshed with quadratic (to prevent locking effects) plain strain triangles using the Python API of the open-source software _Gmsh_[47]. Periodic nodes were generated on opposite sides of the SVE during meshing to allow PBCs to be applied. The generated mesh was then used to generate an input file for the FEM solver _CalculiX_[48]. PBCs were applied, mapping the macro-scale deformation gradient to the SVE. An incremental plasticity model with a multiplicative decomposition of the deformation gradient was used to model the matrix behavior. Linear elastic behavior was assumed until reaching the yield strength, followed by perfectly plastic behavior (accumulating plastic strain at constant stress with further loading) [49]. The matrix plasticity introduces a history-dependency of the microcell response. Fibers were modeled linear elastically. Using _CalculiX_, the micro-scale BVP was solved and the strains and stresses at the integration points were recorded. Subsequently, the volume average strain and stress were calculated as the homogenized properties.

Following [14, 18], the generated training data consisted of both cyclic strain sequences and random strain walks. Strain sequences were generated with a maximum length of 100 strain increments. To apply boundary conditions to the micro-scale BVP, the deformation gradient tensor \(\mathbf{F}_{M}\) was required. Therefore, loading paths, again following [14], were derived from a sequence of right stretch tensors \(\mathbf{U}_{M_{0}},\mathbf{U}_{M_{1}},\ldots\mathbf{U}_{M_{n}}\). At each time step, a load increment \(\Delta\mathbf{U}_{M_{n}}=\mathbf{U}_{M_{n}}-\mathbf{U}_{M_{n-1}}\) was applied. The increment of the right stretch tensor can be represented using orthogonal eigenvectors \(\mathbf{n}_{1}\) and \(\mathbf{n}_{2}\) (only the 2D-case is addressed here) that control the loading direction and eigenvalues \(\Delta\lambda_{1}\) and \(\Delta\lambda_{2}\), governing the load increment size [14]:

\[\Delta\mathbf{U}_{M}=\Delta\lambda_{1}\mathbf{n}_{1}\otimes\mathbf{n}_{1}+\Delta\lambda_{ 2}\mathbf{n}_{2}\otimes\mathbf{n}_{2}. \tag{17}\]

To generate the random loading paths, random orthogonal vectors were generated:

\[\mathbf{n}_{1} =\left[\cos\alpha\quad\sin\alpha\right]^{T} \tag{18}\] \[\mathbf{n}_{2} =\left[-\sin\alpha\quad\cos\alpha\right]^{T},\]

with the angle \(\alpha\) being randomly chosen from a uniform distribution with \(\alpha\in[0,\pi)\). Corresponding eigenvalues were generated as:

\[\Delta\lambda_{1} =\sqrt{\mathcal{R}}\cos\left(\theta\right) \tag{19}\] \[\Delta\lambda_{2} =\sqrt{\mathcal{R}}\sin\left(\theta\right),\]

where \(\mathcal{R}\) and \(\theta\) are uniformly distributed variables in the intervals \(\mathcal{R}\in\left[\Delta R_{min}^{2},\Delta R_{max}^{2}\right]\) and \(\theta\in[0,2\pi]\). This guarantees the total step size to be bounded between \(\Delta R_{min}\) and \(\Delta R_{max}\):

\[\Delta R_{min}\leq\sqrt{\Delta\lambda_{1}^{2}+\Delta\lambda_{2}^{2}}\leq \Delta R_{max}. \tag{20}\]

A lower step size of \(\Delta R_{min}=5\times 10^{-5}\) and an upper limit of \(\Delta R_{max}=2.5\times 10^{-3}\) were used for random strain sequences. Furthermore, the maximum total stretch change was limited to 0.05 for each component of the stretch tensor.

Figure 3: Three random composite microcells with varying FCRs.
Additional to random loading paths, cyclic loading paths were generated, randomly choosing between uniaxial loading in \(x\) or \(y\) direction, biaxial loading (\(x\) and \(y\) loading), or pure shear loading. To generate cyclic strain sequences, the lower step size was increased to \(\Delta R_{min}=5\times 10^{-4}\).

To define the boundary conditions using PBCs, the deformation gradient was calculated from the right stretch tensor:

\[\mathbf{F}_{M}=\mathbf{R}_{M}\cdot\mathbf{U}_{M}, \tag{21}\]

where \(\mathbf{R}_{M}\) is a rotation tensor set to \(\mathbf{R}_{M}=\mathbf{I}\) in this effort due to the frame indifference of the micro-scale BVP [14]. The Green-Lagrange strain tensor \(\mathbf{E}_{M}\) as input to the NN surrogate was then calculated:

\[\mathbf{E}_{M}=\frac{1}{2}\left(\mathbf{U}_{M}^{2}-\mathbf{I}\right). \tag{22}\]

The homogenized response of three random microcells with different FVRs to cyclic loading (Fig. 4a) and random loading (Fig. 4b) can be seen in Fig. 4.

### Dimensionality Reduction of Microstructure Representation

Homogenization of heterogeneous materials requires knowledge about the composition and internal structure of the material. To represent the microstructure a set of descriptors is required that can capture the features governing material behavior while being sufficiently low-dimensional to be tractable [50, 51]. Common examples of microstructure descriptors include volume fraction, inclusion size, or inclusion spacing. In [16], FVR, fiber radius, and mean fiber distance were used to train RNNs to predict the homogenized response of continuously reinforced composites. However, selecting microstructure descriptors intuitively might not sufficiently capture information and, when examining different types of microstructures, might not be generally applicable. Instead, an adaptive microstructure encoding that applies to a wide range of microstructures without manual feature selection is preferable.

It has been shown that the relative spatial distribution of local heterogeneity based on two-point statistics can be used to effectively capture microstructure features [51, 52, 53]. Two approaches will be examined in this effort: a microstructure representation based on dimensionality-reduced two-point statistics and a learned microstructure auto-encoding where a CNN reduces the dimensionality and extracts information from an image of the microstructure. An introduction to two-point statistics will be given first.

Two-point statistics or two-point spatial correlations can give information about the distribution of local states and the fractions of local states in a microstructure [52, 53, 54]. Effectively, two-point spatial correlations provide the probability that the start and end points of a vector are on a certain specified local state, respectively. A discretized (e.g., voxelized) microstructure \(j\) is described by a microstructure function \(m_{j}\left[h;s\right]\) that gives the probability distribution for a local state \(h\in H\), where \(H\) are all possible local states, at each position \(s\in S\) with the complete set of all positions \(S\)[52]. In a voxelized two-dimensional microstructure with uniform cell size the position \(s\) can be simply given by indices \(i\) and \(j\) of the voxels. A set of two-point correlations \(f_{j}\left[h,h^{\prime};r\right]\), providing relative spatial information of the local states \(h\) and \(h^{\prime}\), can be calculated from the correlation of the microstructure function:

\[f_{j}\left[h,h^{\prime};r\right]=\frac{1}{\Omega_{j}\left[r\right]}\sum_{s}m _{j}\left[h;s\right]m_{j}\left[h^{\prime};s+r\right], \tag{23}\]

with a discrete vector \(r\) inside of the microstructure domain and a normalization factor \(\Omega_{j}\left[r\right]\) dependent on \(r\). A simple two-phase microstructure (gray and white cells representing the different phases) with the positions \(s\) is shown in Fig. 5. The two-point statistics can be understood as the likelihood of encountering the first and second phase (\(h\) and \(h^{\prime}\), respectively) of the material at the end (tail) and beginning (head) of a randomly placed vector \(r\), respectively [52]. The matrix-phase auto-correlation of the microcells from the training set shown in Fig. 3 are shown in Fig. 6.

In this effort, a computational implementation of two-point statistics in the Python package _Materials Knowledge Systems in Python (PyMKS)_[52] was used.

To reduce the dimensionality of the microstructure representation provided by two-point statistics, PCA was used [55]. Two-point statistics generate a high-dimensional feature space that can contain redundant information, therefore PCA can be applied to efficiently remove redundant information and reduce the dimensionality [52]. Using PCA, the vector \(\mathbf{f}_{j}\left[l\right]\) containing all possible combinations of \(h\), \(h^{\prime}\), and \(r\) in \(f_{j}\left[h,h^{\prime};r\right]\), can be approximated as:

\[\mathbf{f}_{j}\left[l\right]\approx\overline{\mathbf{f}\left[l\right]}+\sum_{n=1}^{k} \mu_{j,n}\mathbf{\phi}_{n}\left[l\right]. \tag{24}\]

In Eq. 24, \(\overline{\mathbf{f}\left[l\right]}\) represent the average values of all microstructures \(j\) from the calibration data for all \(l\), \(\mu_{j,n}\) are the principal component scores (effectively reduced-dimensionality microstructure descriptors), \(\mathbf{\phi}_{n}\left[l\right]\) are orthonormal
vectors (principal components) that retain the maximum variance under projection, and \(k\) is the number of selected components [52, 55, 56]. The principal component scores \(\mu_{j,n}\) are ordered by their contributions to the variances, with \(n=1\) having the highest contribution.

Randomly choosing 1000 microstructures from the training dataset, calculating two-point statistics and performing PCA, the individual and cumulative contribution to the variance of the first ten principal components are plotted in Fig. 7. The first principal component already explains more than 85 % of the dataset variance, while contributions of further components are quickly decreasing. Therefore, with six components explaining more than 90 % of the variance, PCA combined with two-point statistics is a highly efficient method to reduce the dimensionality of the microstructure description (compared to e.g., a geometric representation of the microstructure).

Figure 4: Homogenized stress-strain response of three random SVE realizations with different FURs for cyclic (4a) and random (4b) loading.
In continuously reinforced composites the predominant microstructure property is the ratio between fiber and matrix content [53]. To gain insight into the meaning of the principal component scores in microstructure representation, the scores of the first three principal components for 1000 random microstructures are plotted against the FVR in Fig. 8. A clear correlation between the first principal component score and the FVR is evident, with all microstructures with the same FVR having similar scores. This suggests that the first PCA score is a direct reflection of the FVR, unrelated to other microstructural characteristics. Again, similar values for the second PCA score are observed in microstructures

Figure 5: Schematic depiction of two-point correlations for a simple discretized microstructure. (adapted from [52])

Figure 6: Matrix phase auto-correlation corresponding to the microstructures shown in Fig. 3.
with similar FVRs. A low score is present for intermediate FVR and both high and low FVR exhibit high scores. However, the third PCA score exhibits variability in microstructures with similar FVR, implying that the FVR is not the dominant factor. Further examination is needed to identify the microstructural parameters represented by the principal component scores. Based on diminishing returns of later principal components, the surrogate model was trained using the first 3 principal components to describe the microstructure. These 3 components accounted for 88 % of the dataset variance.

It should be noted that a set of parameters describing the microstructure (such as FVR) might not necessarily capture features that govern strength behavior. As an example, composite stiffness can generally be well predicted purely based on FVR, while strength is generally highly dependent on local material inhomogeneity [57]. The validity of the

Figure 8: Analysis of principal component scores against the FVR for 1000 random microstructures.

Figure 7: Contribution of the principal components to the variance for a random set of 1000 microstructures.
chosen microstructure representation parameters has, therefore, to be evaluated with respect to their expressiveness for both stiffness and strength prediction.

As an alternative to two-point statistics with PCA for dimensionality reduction, a learned microstructure encoding was implemented. An image of the microstructure (as seen in Fig. 3) was used as input to a CNN. The CNN successfully reduced the dimensionality of the microstructure image to a pre-defined number of parameters (e.g., 3 parameters) and created a hidden state that was used as input to the transformer model. Weights of the CNN were trained simultaneously with the transformer network to minimize the model prediction error. The intended purpose of the microstructure encoding learned by the CNN was to enable the automatic selection of microstructural features that control the homogeneous stiffness and strength response. An advantage of a learned microstructure encoding over two-point correlation could be that the CNN better adapts to features governing composite strength. In contrast, dimensionality-reduced two-point statistics provide pre-defined parameters that are not necessarily optimal to capture property governing features.

The used CNN architecture was roughly based on the AlexNet architecture [58]. Five convolutional blocks were used with a progressively increasing number of output channels produced by the convolution. Maximum pooling layers were added in the second and fifth convolutional blocks for dimensionality reduction and batch normalization was used in the first, second, and fifth blocks to make training more stable. The convolutional blocks were followed by a three-layer feed-forward neural network producing the desired number of outputs. Dropout layers were used after the feed-forward layers to prevent overfitting. Throughout the CNN, Rectified Linear Unit (ReLu) activation functions were used. Black and white microstructure images that were input to the CNN were scaled to 32\(\times\)32 pixels. An overview of the CNN network architecture can be seen in Tab. 1. To assess the efficacy of the CNN in encoding microstructure information using dimensionality-reduced two-point statistics, an equal number of parameters were employed to characterize the microstructure, resulting in a CNN with 3 output parameters.

To analyze microstructure features extracted by the CNN, the output values of the trained network are plotted against the FVR in Fig. 9 for 1000 random microstructures. All three output values show a clear correlation to the FVR with larger output values for larger FVRs. The observed clear correlation shows that the CNN was able to extract governing microstructure features. However, because all output values were related to the FVR, likely a single output value would have been sufficient for the used training data. Introducing more variability into the microstructure, such as varying microcell size and fiber diameter, non-constant fiber diameter, or differently shaped inclusions, will require a larger number of parameters to describe the microstructure features.

## 4 Surrogate Model Architecture and Training Procedure

In this effort, a transformer network was chosen as a surrogate to predict the history-dependent response of composite materials. While the original transformer [19] consisted of an encoder-decoder architecture, only the decoder was used to build the surrogate model, similar to the GPT models [23]. The decoder's masked self-attention (see Sec. 2)

\begin{table}
\begin{tabular}{c c c c} \hline \hline Layer nr & Type and nonlinearity & Input size & Output size \\ \hline \multirow{4}{*}{1} & Input & & 32 \(\times\) 32 \(\times\) 1 \\  & Convolution (5\(\times\)5, 16 filters), Batch normalization, ReLu & & 32 \(\times\) 32 \(\times\) 1 & 32 \(\times\) 32 \(\times\) 16 \\  & Convolution (3\(\times\)3, 32 filters), Batch normalization, ReLu, Max Pooling (2\(\times\)2) & 32 \(\times\) 32 \(\times\) 16 & 16 \(\times\) 16 \(\times\) 32 \\  & Convolution (3\(\times\)3, 64 filters), Batch normalization, ReLu & 16 \(\times\) 16 \(\times\) 32 & 16 \(\times\) 16 \(\times\) 64 \\  & Convolution (3\(\times\)3, 128 filters), Batch normalization, ReLu, Max Pooling (2\(\times\)2) & 16 \(\times\) 16 \(\times\) 64 & 16 \(\times\) 16 \(\times\) 64 \\  & Flatten & 8 \(\times\) 8 \(\times\) 128 & 8192 \\  & Fully connected, ReLu, Dropout (0.3) & 8192 & 2048 \\  & Fully connected, ReLu, Dropout (0.3) & 2048 & 1024 \\  & Fully connected, Batch normalization & 1024 & 3 \\ \hline \hline \end{tabular}
\end{table}
Table 1: CNN architecture used to extract microstructure information from images.
prevents the network from seeing "future", i.e., yet unknown, strain values in the input sequence. This guarantees that the predicted stress is only a function of the deformation history up to the current strain increment.

Transformers require an embedding to convert the input and the target sequence to vectors matching the transformer dimension. In the input path, the strain sequence, stress sequence, and the vector containing the microstructure information (either PCA scores or outputs from CNN) were linearly transformed to match the transformer dimension \(d_{model}\). Two stacked Gated Residual Networks (GRNs) were used as encoders and decoder to embed static (time-independent) microstructure parameters. The GRN was adapted from the _Temporal Fusion Transformer_[38]. In addition to the primary input (strain or stress sequence), the GRN has an optional context vector (in our case the microstructure information) as input. By using gating mechanisms, the GRN can filter out insignificant inputs and can skip the GRN entirely if not relevant to the output [38]. A schematic representation of the GRN architecture is shown in Fig. 10. After passing through a linear layer, input and context are added and Exponential Linear Unit (ELU) activation [59] function is applied. This is followed by a second linear layer, dropout, and a Gated Linear Unit (GLU) [60] gating mechanism:

\[\text{GLU}\left(\boldsymbol{\gamma}\right)=\sigma\left(\boldsymbol{W}_{1} \boldsymbol{\gamma}+\boldsymbol{b}_{1}\right)\odot\left(\boldsymbol{W}_{2} \boldsymbol{\gamma}+\boldsymbol{b}_{2}\right), \tag{25}\]

where \(\boldsymbol{\gamma}\) is the input to the GLU, \(\boldsymbol{W}_{1}\), \(\boldsymbol{W}_{2}\), \(\boldsymbol{b}_{1}\), and \(\boldsymbol{b}_{2}\) are the weights and biases, respectively, \(\sigma\left(\cdot\right)\) is the sigmoid activation function, and \(\odot\) the element-wise Hadamard product [38]. A residual connection adds the GRN input and output, allowing the GRN to be skipped entirely. Finally, the output is normalized. It has been shown that residual connections improve network training and performance, and allow for far deeper networks [43]. The input and output size of all GRNs was chosen to be similar to \(d_{model}\) (512) and the dense layers in the GRN were of dimension 2048.

After the GRN encoders, the microstructure information was added again to the input path, followed by normalization and positional encoding. The encoded source (strain) and right-shifted target (stress) sequences were then passed through the transformer decoder. The right-shifted stress sequence required as input to the decoder was generated by inserting a zero vector at the first position of the stress sequence and truncating the sequence at the second to last element. After the transformer decoder, two stacked GRNs with the microstructure as context were used followed by a linear transformation to generate the stress predictions. The neural network architecture is shown in Fig. 11. The model architecture was implemented in _PyTorch_[61] using the _PyTorch Lightning_ interface. Surrogate model parameters are summarized in Tab. 2.

With CNN microstructure encoding, the chosen network had a total of 66.3 million trainable parameters. The total trainable parameters consisted of 19 million parameters for the CNN encoding, 25.2 million parameters for the transformer decoder, and 7.4 million parameters for each GRN encoder (including the linear transformation of the strain and stress sequences, respectively) and 7.3 million parameters for the GRN decoder.

Following [19], the transformer model was trained using the Adam optimizer [62] with learning rate warmup. During training, the learning rate \(lr\) was initially increased linearly, followed by decreasing learning rate proportional to the

Figure 9: Output values of the trained CNN for 1000 random microstructures with varying FVR.
inverse square root of the current step number [19]. The learning rate was calculated by:

\[lr=d_{model}^{-0.5}\min\left(step_{num}^{-0.5},step_{num}\cdot warmup_{steps}^{-1. 5}\right), \tag{26}\]

where \(d_{model}\) is the attention mechanism dimension, \(step_{num}\) is the current training step number and \(warmup_{steps}\) is the number of steps of linearly increasing learning rate. As was chosen in [19], \(warmup_{steps}\) was set to 4000. Using learning rate warmup can prevent divergent training in transformer models [63].

Target and source sequences were standardized to zero mean and unit variance before training to speed up the training process. The training sequences were truncated at random indices to allow the transformer to adapt to different length sequences, resulting in sequence lengths between one and 100. As a loss function, the Mean Squared Error (MSE) between the predicted and true values was used. It should be noted that due to the choice of loss function, the training process was dominated by high-stress data [14]. Therefore, relative prediction performance might decrease for smaller stresses. This effect should be reduced through normalization of the data and the loss function was, therefore, deemed suitable.

### Training the Surrogate Model on Homogeneous Data

To test the feasibility of a transformer network as a surrogate model for history-dependent plasticity, initially, the model was trained on the response of a purely homogeneous elastoplastic material. Using the strain sequence generation described in Sec. 3.2 and a (semi) analytical implementation of incremental J2-plasticity, a large amount of data could be generated at a low computational cost. This way and using the matrix material parameters listed in Tab. 3, 500,000 stress-strain sequences with a maximum sequence length of 100 were generated. Half the generated sequences consisted of random strain walks while the remaining sequences were cyclic. During training, 80 % of the data were used for training and the remaining 20 % were used for model validation. The surrogate model architecture remained

\begin{table}
\begin{tabular}{l c} \hline \hline Model parameter & Size \\ \hline Transformer size (\(d_{model}\)) & 512 \\ Transformer feed-forward size & 2048 \\ Number of decoder layers & 6 \\ Number of attention heads & 8 \\ Dropout probability & 0.1 \\ GRN encoder / decoder size & 2048 \\ Number of stacked GRNs in encoder / decoder & 2 \\ GRN input / output size & 512 \\ \hline \hline \end{tabular}
\end{table}
Table 2: Transformer architecture parameters.

Figure 10: Architecture of the GRN. (adapted from [38])
unchanged except for the GRN context input (microstructure) and the microstructure addition step after the GRN encoders were removed. The training was performed for 800 epochs on four _NVIDIA RTX A6000_ GPUs with a batch size of 500. Total training time was approximately 24 hours.

Figure 11: Architecture of the transformer network to predict the history-dependent response of heterogeneous composites using a CNN to extract features from microstructure images.
After 800 epochs, the lowest validation MSE was \(5.86\times 10^{-4}\). Note that the MSE is for the standardized dataset and not the actual stress prediction. Subsequently, an additional 500 random and 500 cyclic strain sequences were generated to test the model performance on previously unseen data. Using the trained network the average Root-Mean-Square Error (RMSE) of the stress predictions for the ground truth were calculated for the unscaled (non-standardized) data to be 0.804 MPa. Comparisons between the network predictions and the ground truth data can be seen in Fig. 12 for a cyclic strain sequence and in Fig. 13 for a random strain sequence. History-dependence of the material response for both cyclic (Fig. 12) and random strain (Fig. 13) sequences was correctly predicted by the network with only some minor visible differences between the predicted and the true values.
### Adapting the Pre-Trained Surrogate Model for Homogenization

Following the training of the surrogate model for a homogeneous elastoplastic material, a total of 330,000 homogenized stress-strain sequences were generated following the homogenization procedure in Sec. 3.2. Each strain sequence had a length of 100 strain increments. Two-dimensional SVEs modeling the transversal behavior of continuously reinforced composites assuming plane strain conditions were used. The FVR of the generated SVEs was varied between 0.2 and 0.5 (the true values vary somewhat, as the target value cannot necessarily be achieved by randomly placing circular fibers inside of a square microcell) with fixed microcell size and fiber diameter. Around 40 % of the generated sequences were random walks and the remaining sequences had cyclic loading paths. The material parameters were left constant. An overview of the material parameters and SVE properties used to generate the training data are

Figure 13: Random input strain sequence (13a) and comparison between surrogate model prediction and ground truth (13b).
listed in Tab. 3. Simulating 330,000 strain sequences took approximately 4 weeks on two _AMD EPYC 7513 32-Core_ processors.

Using the generated dataset, two surrogate models with different microstructure representations (Sec. 3.3) were then trained. One of the models used a CNN with 3 output parameters, while the other had as input the first 3 PCA scores of the respective microstructure. Both models were initialized using the pre-trained weights from the homogeneous model and subsequently trained for 50 epochs. Because large batch sizes commonly lead to degradation in generalization performance [64], a smaller batch size of 100 was used. Furthermore, the learning rate was reduced by an order of magnitude to avoid "forgetting" all previously learned information. Other training parameters remained unchanged from the ones listed in Sec. 4.1.

After 50 training epochs, the lowest MSE on the validation dataset was \(2.30\times 10^{-3}\) for the standardized data after epoch 39 when using a CNN for microstructure feature extraction, and \(3.04\times 10^{-3}\) after epoch 15 using two-point correlations. Longer required training time to reach minimum MSE when using the CNN can be explained by the randomly initialized CNN weights, while the transformer weights were initialized using the pre-trained weights. Therefore, during training likely mainly the CNN parameters had to be optimized. Larger MSE compared to homogeneous data (minimum homogeneous validation MSE \(8.14\times 10^{-5}\)) showed the difficulty of the network to adapt to heterogeneous data. This was expected to a certain degree, as the homogenization problem represents inherently complex compound elastoplastic behavior, with the fibers storing significant amounts of elastic energy while plastic deformation occurs in the matrix phase [16].

To evaluate the impact of using pre-trained weights, the surrogate model using PCA was also trained using random weight initialization. The minimum validation loss with random initialization of \(3.34\times 10^{-3}\) was reached after 46 epochs, representing considerably longer required training time and a marginally worse minimum loss.

A test dataset containing 3000 strain sequences for random microstructures was generated to validate the performance of the trained surrogate model for previously unseen data. The RMSE between stress predictions and FEM results on the test data was 2.16 MPa using the CNN, 1.51 MPa using PCA with pre-training, and 2.63 MPa using PCA with random weight initialization. Stresses in the test dataset ranged between -387.5 MPa and 360.9 MPa. Both microstructure representation methods were able to successfully relate microstructure information to homogenized properties, with somewhat lower prediction error when using the PCA. Random weight initialization also resulted in worse performance compared to using the pre-trained weights from homogeneous training. The prediction error incrementally rises with increasing sequence length, due to the transformer's approach of predicting one time-step at a time, relying on the preceding data (see Fig. 14).

Fig. 15 demonstrates the comparison between cyclic FEM sequences and the surrogate model performance using both microstructure encoding methods (CNN Fig. 15a, two-point statistics Fig. 15b). Generally, good accuracy between the FEM and predicted curves can be seen. The surrogate models correctly capture changing stiffness with varying FVR and accurately predict the material history-dependence.

Predictions of the surrogate model using two-point correlations for a random strain sequence are depicted in Fig. 16. Generally, a good agreement between surrogate model predictions and FEM results can be seen, including for the out-of-plane stress \(\sigma_{zz}\) component.

However, several strain sequences from the test dataset, particularly random walks, showed relatively large prediction errors. The fact that the majority of training sequences were cyclic loading paths might be a contributing factor to these outliers. Furthermore, as transformers commonly require large amounts of data [33], increasing the number of training sequences could improve prediction performance. In this effort, a relatively large model (around 50 million parameters) was selected. Large models, when trained on small datasets, are more prone to overfitting, leading to poor

\begin{table}
\begin{tabular}{l c} \hline \hline SVE size (\(\mu\)m) & 70 \\ Fiber diameter (\(\mu\)m) & 7 \\ FVR (-) & 0.2 –0.5 \\ Young’s modulus matrix (MPa) & 1000 \\ Poisson’s ratio matrix (-) & 0.4 \\ Yield stress matrix (MPa) & 20 \\ Young’s modulus fiber (GPa) & 200 \\ Poisson’s ratio fiber (-) & 0.18 \\ \hline \hline \end{tabular}
\end{table}
Table 3: Material and microcell parameters used to generate the training data.
generalization performance. Optimization of the network size and hyperparameter tuning could further improve the prediction performance of the surrogate models, reduce the training duration, and lower the required training data.

Three principal components and three output parameters of a CNN were used to encode microstructure information as input to the surrogate model. However, it was shown (see Sec. 3.3) that likely a single parameter, namely the encoded FVR, would have been sufficient for the analyzed microstructures. Increasing the complexity of the microstructures used for training by, for example, varying the RVE size and fiber diameter, or introducing non-circular inclusions will require more parameters to encode RVE features. Similarly, higher fidelity damage models to account for matrix damage will also require more refined structure descriptions to account for fiber clusters and voids initiating and accumulating local damage.

Comparing computational cost, solving a single microstructure subject to a random strain sequence with a length of 100 increments using parallel FEM on two CPUs took approximately 40 s. In contrast, predicting the response of 1000 microcells on a _NVIDIA RTX A6000_ GPU using the surrogate model took less than 25 s. Due to the fast batch processing of the neural network, major time advantages can especially be achieved when predicting a large number of sequences, as would be the case in structural FE\({}^{2}\) analysis.

### Testing Prediction Performance on Unseen Microstructure

To acquire additional insight into the knowledge assimilated by the neural network and test the model's robustness, a completely new microstructure containing a rectangular inclusion (see Fig. 17) was simulated under cyclic loading in \(x\) and \(y\)-directions using FEM and the results were compared with the network predictions. The rectangular inclusion was centered inside the microcell with a FVR of 0.2 and a width (\(x\)) to height (\(y\)) ratio of \(0.25:0.8\). Due to the elongated inclusion shape, anisotropy of the response was introduced (see Fig. 18). Because the surrogate models were only trained on circular inclusions, the prediction performance for new microstructures can give information about the generalizability of the models.

Comparing the simulated data with the neural network, both microstructure encoding methods were unable to capture the microcell's anisotropic response. Instead, the surrogate models predicted almost identical responses for transversal and longitudinal loading. This was expected as the networks were trained on data that was mostly isotropic and the network predictions were predominantly based on the FVR. However, the CNN seemed to be able to extract information about the FVR even for completely unseen microstructures (see Fig. 18a). Using two-point statistics, the transversal response of the microstructure was captured relatively well (see Fig. 18b) with a larger error for longitudinal loading. This shows that the surrogate models were able to learn a FVR-based homogenization for elastoplastic

Figure 14: Average RMSE as a function of number of strain increments using two-point correlation microstructure encoding.
composites. Additional training on a wider variety of microstructures would allow the models to generalize and give improved predictions for a wide range of inclusion geometries.

## 5 Concluding Remarks

A surrogate model to predict the history-dependent homogenized response of elastoplastic composite SVEs based on a NLP transformer neural network has been presented. A decoder-only architecture with GRNs to encode time-dependent (strain, stress) and static data (microstructure information) was chosen. Features were generated from mages

Figure 15: Comparison between surrogate model prediction and FEM for cyclic loading paths at three different FVRs using CNN (15a) and PCA (15b).
of the material microstructure comparing (i) a learned auto-encoding based on a CNN and (ii) two-point correlation functions using PCA for dimensionality reduction.

After training the developed network using random strain walks and cyclic loading paths, the model successfully predicted the history-dependent response of previously unseen composite SVEs. Furthermore, both microstructure representation techniques were able to extract parameters from the microstructure that governed the material response.

Considerable savings in computational cost for microstructure homogenization compared to FEM were demonstrated. Efficiency gains compared to FE\({}^{2}\) make UQ of composite structures with randomness in the underlying constituent distribution attainable. The transformer surrogate model could be used for scale transition to estimate the local ho

Figure 16: Random input strain sequence (16b) and comparison between surrogate model prediction using PCA microstructure representation and ground truth (16b).
mogenized material response based on the micro-scale constituent composition and included in a Monte Carlo (MC) or Quasi Monte Carlo (QMC) loop to sample the response distribution.

The current implementation was restricted to simple two-dimensional composite SVEs of fixed size and fiber diameter and fixed constituent properties. Having shown the feasibility of using a transformer surrogate for homogenization, the model will be extended to account for more complex microstructures, such as randomly sized SVEs and non-circular inclusions, and varying material parameters. Extension to three-dimensional microstructures is possible by using a three-dimensional CNN or three-dimensional two-point correlations to encode the microstructure.

In [2], the authors implement a data-driven material discovery approach using feature importance scores of microstructure parameters to more efficiently cover the training space, thereby escaping the curse of dimensionality. Following their approach, the amount of required training data to extend the transformer model to more complex microstructures could be reduced by carefully designing experiments based on feature importance. Furthermore, physics-informed loss functions could be used to decrease the amount of required training data [65].

## Funding Sources

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Figure 17: Microcell with rectangular inclusion.