Mixed formulation of physics-informed neural networks for
  thermo-mechanically coupled systems and heterogeneous domains
Physics-informed neural networks (PINNs) are a new tool for solving boundary
value problems by defining loss functions of neural networks based on governing
equations, boundary conditions, and initial conditions. Recent investigations
have shown that when designing loss functions for many engineering problems,
using first-order derivatives and combining equations from both strong and weak
forms can lead to much better accuracy, especially when there are heterogeneity
and variable jumps in the domain. This new approach is called the mixed
formulation for PINNs, which takes ideas from the mixed finite element method.
In this method, the PDE is reformulated as a system of equations where the
primary unknowns are the fluxes or gradients of the solution, and the secondary
unknowns are the solution itself. In this work, we propose applying the mixed
formulation to solve multi-physical problems, specifically a stationary
thermo-mechanically coupled system of equations. Additionally, we discuss both
sequential and fully coupled unsupervised training and compare their accuracy
and computational cost. To improve the accuracy of the network, we incorporate
hard boundary constraints to ensure valid predictions. We then investigate how
different optimizers and architectures affect accuracy and efficiency. Finally,
we introduce a simple approach for parametric learning that is similar to
transfer learning. This approach combines data and physics to address the
limitations of PINNs regarding computational cost and improves the network's
ability to predict the response of the system for unseen cases. The outcomes of
this work will be useful for many other engineering applications where deep
learning is employed on multiple coupled systems of equations for fast and
reliable computations.
Ali Harandi, Ahmad Moeineddin, Michael Kaliske, Stefanie Reese, Shahed Rezaei
Mixed formulation of physics-informed neural networks for thermo-mechanically coupled systems and heterogeneous domains

###### Abstract

Deep learning methods find a solution to a boundary value problem by defining loss functions of neural networks based on governing equations, boundary conditions, and initial conditions. Furthermore, the authors show that when it comes to many engineering problems, designing the loss functions based on first-order derivatives results in much better accuracy, especially when there is heterogeneity and variable jumps in the domain [1]. The so-called mixed formulation for PINN is applied to basic engineering problems such as the balance of linear momentum and diffusion problems. In this work, the proposed mixed formulation is further extended to solve multi-physical problems. In particular, we focus on a stationary thermo-mechanically coupled system of equations that can be utilized in designing the microstructure of advanced materials. First, sequential unsupervised training, and second, fully coupled unsupervised learning are discussed. The results of each approach are compared in terms of accuracy and corresponding computational cost. Finally, the idea of transfer learning is employed by combining data and physics to address the capability of the network to predict the response of the system for unseen cases. The outcome of this work will be useful for many other engineering applications where DL is employed on multiple coupled systems of equations.

keywords: Physically informed neural networks, Thermo-mechanically coupled problems, Heterogeneous solids +
Footnote †: journal:

## 1 Introduction

Deep learning (DL) methods are a branch of machine learning (ML) that have remarkable flexibility in approximating arbitrary functions and operators [2; 3]. They are able to perform the latter task by finding the hidden patterns for the training data set (also known as supervised learning). Moreover, DL methods are capable of discovering the solution to a boundary value problem (BVP) by solely considering physical constraints (unsupervised learning). Therefore, DL models find their way into many engineering applications from material design to structural analysis [4]. One advantage of ML algorithms compared to the classical solver is the huge speedup that one achieves after successful training. The efficiency of these models makes them extremely relevant for multi-scaling analysis where the bottleneck is in the process of information
transformation between the upper scale and the lower scale (see [5; 6; 7] and references therein). An important issue is that the network has to be trained for enough observations considering all the relevant scenarios. Otherwise, one cannot completely rely on the network outcome. Adding well-known physical constraints to the data training process improves the network's prediction capabilities (see [1] and references therein). The latter point needs to be investigated properly for different test cases and physical problems to gain the most advantage from applying the physical laws to the classical data-driven approaches.

Here, we intend to examine the potential of the DL method for the field of computational engineering, where we deal with a complex heterogeneous domain as well as a coupled system of governing equations. In what follows, we shall review some of the major contributions and address the new contributions of the current study.

Training ML algorithms based on the available physical data obtained from our numerical solvers can speed up the process of predicting the solution. Park et al. [8] utilized a multi-scale kernel neural network to enhance the prediction of strain fields for unseen cases. In traditional ML approaches, where NNs are only trained based on data, the underlying physics of the problem is embedded in the data itself, which is the result of other numerical solvers. Yang et al. [9] employed the DL method to predict complex stress and strain fields in composites. Khorrami et al. [10] employed a convolutional neural network to predict the stress fields for history-dependent material behavior. See also [11], where authors employed a Fourier neural operator to predict a stress field at a high resolution by performing the training based on data of low resolution. NN surrogate models are powerful models to solve complex systems of equations in coupled thermal, chemical, and hydrological engineering problems [12]. Data-driven models are used to estimate the state of charge based on the driving conditions [13] or as digital twins of complex systems like proton exchange membrane fuel cells [14]. See also [15; 16] for a comparison between the performance of different neural operators which are suitable for industrial applications.

Interestingly enough, following the idea of physics-informed neural networks (PINNs) [17], and by employing the physical laws in the final loss functions, one can reasonably rely on the network outcome. This is due to the fact that the network outcomes automatically respect and satisfy important physical constraints. However, a proper and in-depth theoretical understanding of the convergence and generalization properties of the PINNs method is still under investigation. Therefore, the idea behind PINNs is also further extended and explored by several authors to improve the network performance even for unseen scenarios. For some examples see [18; 19; 20]. In case the underlying physics of the problem is completely known and complete, the NN can be trained without any initial data and solely based on the given BVP [1]. The latter point is also referred to as unsupervised learning since no initial solution to the problem is needed beforehand. Hu et al. [21] initiate discussions on how extended PINNs that work based on domain decomposition method show effectiveness in modeling multi-scale and multi-physical problems. The authors introduce a tradeoff for generalization. On the one hand, the decomposition of a complex PDE solution into several simple parts decreases the complexity and on the other hand, decomposition might lead to overfitting and the obtained solution may become less generalizable, see also investigations by Jagtap et al. [19], D. Jagtap and Em Karniadakis [22], Henkes et al. [23] and Wang et al. [24] related to the idea of domain decomposition. Wang et al. [25] investigated how PINNs are biased towards learning functions along the dominant eigen-directions of the neural tangent kernel and proposed architectures that can lead to robust and accurate PINN models. Different loss terms in the PINNs formulation may be treated differently on their importance. McClenny and Braga-Neto [26] proposed to utilize adaptive
weights for each training point, so the neural network learns which regions of the solution are difficult and is forced to focus on them. The basic idea here is to make the weights increase as the corresponding losses increase. Xiang et al. [27] also addressed that the PINNs' performance is affected by the weighting of losses and established Gaussian probabilistic models to define the self-adaptive loss functions through the adaptive weights for each loss term. See also [28; 29]. Another promising extension is the idea of reducing the differential orders for the construction of loss functions. For this purpose, researchers combined different ideas from the variational formulation of partial differential equations and apply them to various problems in engineering applications [1; 23; 30; 31].

One important aspect of engineering applications is to consider the multi-physical characteristics. In other words, in many realistic problems, various coupled and highly nonlinear fields have to be considered simultaneously. As an example in solid mechanics, besides the mechanical deformation, often one needs to take into account the influence of the thermal field [32], chemical concentrations [33], electrical and/or magnetic fields [34], damage field [35; 36], and others.

The idea of employing PINNs for solving underlying equations within the multi-physical context is addressed by several authors. In such problems, we are dealing with different loss terms which represent different physics and might be completely different from each other in terms of the magnitude of numerical values. Therefore, a balance between different loss terms is often necessary. Therefore, classical PINN models need to be further improved. Amini et al. [37] investigated the application of PINNs to find the solution to problems involving thermo-hydro-mechanical processes. The authors proposed to use a dimensionless version of the governing equations as well as a sequential training strategy to obtain better accuracy in the multiobjective optimization problem. Raj et al. [38] studied the thermo-mechanical problem using PINN for functionally graded material and the ability of the network to approximate the displacements and showed that the network leads to an order higher error when it comes to approximating the stress fields. Laubscher [39] presented a single and segregated network for PINNs dry air humidification fluid properties heat diffusion to predict momentum, species, and temperature distributions of a dry air humidification problem. It is reported that the segregated network has lower losses when compared to the single network. Mattey and Ghosh [40] also reported that the PINN's accuracy suffers for strongly non-linear and higher-order PDEs such as Allen Cahn and Cahn Hilliard equations. The authors proposed to solve the PDE sequentially over successive time segments using a single neural network. See also studied by Chen et al. [41], Nguyen et al. [42] and Bischof and Kraus [28]. In modeling crack propagation in solids via a smeared approach, one usually assigns a field variable to damage. In this case, equations for the evolution of the damage field are usually coupled to those from the mechanical equilibrium [36]. Goswami et al. [35] proposed a physics-informed variational formulation of DeepONet for brittle fracture analysis, see also investigations by [43; 44].

The computational efficiency and cost of the training can be improved by utilizing a transfer learning algorithm [45]. The key idea of transfer learning is to capture the common features and behavior of the source domain and transfer it to a new target domain. To do so, one can transfer the pre-trained hyperparameters of the NN and initialize them as hyperparameters of the NN for the new domain. The latter decreases the computational cost of training [1; 46]. Tang et al. [47] utilized transfer learning-PINNs to solve vortex-induced vibration problems. Goswami et al. [48] employed the idea of transfer learning to partially retrain the network's hyperparameters to identify the crack patterns in a domain. Pellegrin et al. [49] exploited the transfer learning idea to solve a system of equations by having various configurations (initial conditions) and showed
that utilizing the trained multi-head network leads to a significant speed-up for branched flow simulations without losing accuracy. Desai et al. [50] utilized transfer learning to solve linear systems of PDEs such as Poisson and Schrodinger equations by PINNs. Gao et al. [51] proposed a singular-valued-based transfer learning algorithm for PINNs and showed its performance for Allen Cahn equations.

According to the reviewed articles, it becomes clear that having fast and reliable results by employing the PINNs method for engineering applications has to be done with care. The naive version of PINN may not converge to the correct solution, and even so, the training process may take significantly more time compared to other numerical solvers based on the finite element method. This work discusses ideas on how to overcome the previously mentioned issues. Here, we intend to focus on the idea of first-order PINNs formulation for a stationary thermo-mechanically coupled system in heterogeneous solids and compare the performance of the coupled and sequential training (see also Fig. 1). In Section 2, we provide the derivation of the governing equations and their variational form. In Section 3, the architecture of the new mixed formulation of PINNs is described. The results and conclusions are then provided in Sections 4 and 5, respectively.

## 2 Formulation of the problem

### Stationary thermo-elasticity in heterogeneous solids

The formulation behind the thermo-mechanical problem is first summarized in what follows. In this problem, we are dealing with the displacement field denoted by the vector \(\mathbf{u}\) and the temperature field denoted by the scalar parameter \(T\). Through different coupling terms,

Figure 1: Sequential training and coupled training algorithms for solving a multi-physical boundary value problem as well as the mixed PINN architecture. \(\mathcal{L}_{\mathrm{Phy}}\) stands for the relevant physics law. The variable \(\mathbf{x}\) denotes the location of collocation points. The variable \(\mathbf{u}_{i}\) shows the \(i\)-th primary variable (solution) and \(\nabla\mathbf{u}_{i}^{o}\) is their spatial gradient. Separated neural networks are then utilized to predict primary variables and their spatial gradients.
these two fields influence each other. Starting with kinematics, the total strain \(\mathbf{\varepsilon}\) is additively decomposed into the elastic part and the thermal part as

\[\mathbf{\varepsilon}\,=\,\mathbf{\varepsilon}_{e}\,+\,\mathbf{\varepsilon}_{t}=\,\text{sym} \left(\text{grad}(\mathbf{u})\right)\,=\,\frac{1}{2}\left(\nabla\mathbf{u}\,+\,\nabla \mathbf{u}^{T}\right). \tag{1}\]

In Eq. (1), \(\mathbf{\varepsilon}_{e}\) represents the elastic strain tensor, and \(\mathbf{\varepsilon}_{t}\) presents the thermal strain caused by the temperature field. Considering isotropic material behavior, the thermal strain tensor follows a linear expansion law and reads

\[\mathbf{\varepsilon}_{t}\,=\,\alpha(x,y)\left(T-T_{0}\right)\,\mathbf{I}. \tag{2}\]

Here, \(\alpha\) is the thermal expansion coefficient and it can vary spatially as we are dealing with a heterogeneous microstructure. Moreover, \(T_{0}\) is the initial temperature and \(\mathbf{I}\) is the second-order identity tensor. Utilizing the fourth-order elasticity tensor \(\mathbb{C}\), one can write the stress tensor as

\[\mathbf{\sigma}\,=\mathbb{C}(x,y)\,\left(\mathbf{\varepsilon}-\mathbf{\varepsilon}_{t} \right). \tag{3}\]

Note, that at this point, we assume the elastic properties to be temperature-independent. Considering temperature-dependent material properties shall be investigated in future work and the investigations in the current work can easily be extended to such problems. Writing Eq. (3) in Voigt notation reads

\[\hat{\mathbf{\sigma}}\,=\hat{C}(x,y)\,\,\left(\hat{\mathbf{\varepsilon}}\,-\,\hat{\bm {\varepsilon}}_{t}\right). \tag{4}\]

In Eq. (4), \(\hat{\bullet}\) presents the tensor (\(\bullet\)) in the Voigt notation. For the two-dimensional plane strain setup, the position-dependent elasticity tensor \(\hat{\mathcal{C}}\) is written as

\[\hat{C}(x,y)\,=\,\frac{E(x,y)}{(1-2\nu(x,y))(1+\nu(x,y))}\,\left[\begin{array} []{ccc}1-\nu(x,y)&\nu(x,y)&0\\ \nu(x,y)&1-\nu(x,y)&0\\ 0&0&\frac{1-2\nu(x,y)}{2}\end{array}\right]. \tag{5}\]

Here, Young's modulus \(E\) and Poisson's ratio \(\nu\) represent the elastic constants of the material. The isotropic assumption holds for every point of material. It is worth mentioning that having a heterogeneous domain causes the dependency of material parameters on the coordinates of collocation points. For the mechanical field, the balance of linear momentum by having no body force vector, and the Dirichlet and the Neumann boundary conditions read

\[\text{div}(\mathbf{\sigma}) =\mathbf{0}\quad\text{ in}\quad\Omega, \tag{6}\] \[\mathbf{u} =\bar{\mathbf{u}}\quad\text{on}\,\,\,\,\Gamma_{D},\] (7) \[\mathbf{\sigma}\cdot\mathbf{n}=\mathbf{t} =\bar{\mathbf{t}}\quad\text{ on}\,\,\,\Gamma_{N}. \tag{8}\]

Next, we write the equation of energy balance in the absence of a heat source, and for a steady-state problem ([32])

\[\text{div}(\mathbf{q}) =0\quad\text{ in }\Omega, \tag{9}\] \[T =\bar{T}\quad\text{on }\bar{\Gamma}_{D},\] (10) \[\mathbf{q}\cdot\mathbf{n} =q_{n} =\bar{q}\quad\text{ on }\bar{\Gamma}_{N}. \tag{11}\]
In the above equations, \(\Omega\), and \(\Gamma\) are the material points in the body and on the boundary of the mechanical field, respectively. The temperature boundary condition is also satisfied on \(\bar{\Gamma}\). In Eqs. (9)  and (11), \(\mathbf{q}\) denotes the heat flux and is defined based on Fourier's law as

\[\mathbf{q}=-k(x,y)\,\nabla T. \tag{12}\]

Here, \(k(x,y)\) shows the position-dependent heat conductivity coefficient.

To obtain the Galerkin weak form (WF) for Eq. (6), the multiplication with a test function \(\delta\mathbf{u}\) is done, which also satisfies the Dirichlet boundary conditions of the mechanical field. After integration by parts and employing the Neumann boundary conditions, one finally gets

\[\int_{\Omega}\delta\mathbf{\hat{\varepsilon}}_{e}^{T}\,\hat{C}(x,y)\, \mathbf{\hat{\varepsilon}}_{e}\,\,\,dV\,-\,\int_{\Gamma_{N}}\delta\mathbf{u}^{T}\, \bar{\mathbf{t}}\,\,dA=0. \tag{13}\]

Considering the balance between the mechanical internal energy \(E_{\rm int}^{M}\) and the mechanical external energy \(E_{\rm ext}^{M}\), we write the weak form in Eq. (13) as a minimization of the total energy, see [30], as

\[\text{Minimize:}\,\,\underbrace{\int_{\Omega}\frac{1}{2}\mathbf{\hat {\varepsilon}}_{e}^{T}\,\hat{C}(x,y)\,\mathbf{\hat{\varepsilon}}_{e}\,\,\,dV}_{E _{\rm int}^{M}}-\underbrace{\int_{\Gamma_{N}}\mathbf{u}^{T}\,\bar{\mathbf{t}}\,\,dA}_ {E_{\rm ext}^{M}}, \tag{14}\]

which is subjected to BCs in Eq. (7) and Eq. (8). By multiplying the Eq. (9) with a test function \(\delta T\) and following a similar procedure for the mechanical field and applying the boundary condition in Eqs. (10)  and 11, the Galerkin weak form for the thermal field leads to

\[\int_{\Omega}\,k(x,y)\,\nabla^{T}T\,\delta(\nabla T)\,\,dV\,+\, \int_{\bar{\Gamma}_{N}}\bar{q}\,\,\delta T\,\,dA\,=\,0. \tag{16}\]

One can also reformulate the weak form in Eq. (16) as a variation of the energy. The latter initiates the balance of thermal internal and external energies \(E_{\rm int}^{T}\) and \(E_{\rm ext}^{T}\), respectively, which leads to

\[\text{Minimize:}\,\,\underbrace{\int_{\Omega}\frac{1}{2}\,k(x,y)\, \nabla^{T}T\,\nabla T\,\,\,dV}_{E_{\rm int}^{T}}+\underbrace{\int_{\bar{ \Gamma}_{N}}\bar{q}\,T\,\,dA}_{E_{\rm ext}^{T}}. \tag{17}\]

The above equation is minimized subjected to BCs in Eq. (10) and Eq. (11).

It is important to note that in Eq. (14), the minimization is performed for both the displacement field \(\mathbf{u}\) and temperature field \(T\) in the coupled training process. However, by using sequential training and fixing one field variable, the referring loss function is minimized with respect to the displacement field \(\mathbf{u}\). For Eq. (16), in the absence of the other fields in the formulation, the minimization is done with respect to the temperature field \(T\) in both training procedures.

Similarly to the previous works [1, 30], expressions in Eqs. (13)  and (16) are added to the neural network's loss function for the primary variables \(\mathbf{u}\) and \(T\), respectively.

Two different sets of geometries according to Fig. 2 are considered. The first geometry is selected to represent one of the most common cases in engineering problems, the fiber matrix
setup. The second geometry stands for an advanced engineering alloy microstructure that contains rectangular inclusions. Properties of phase 1 (matrix) are denoted by the sub-index "mat" while the material properties for phase 2 (inclusion) are denoted by the sub-index "inc". Next, we define the ratio \(E_{\text{mat}}/E_{\text{inc}}=k_{\text{mat}}/k_{\text{inc}}=R\). The value \(R\) is equal to \(0.3\) and \(2\) for the left and right geometries, respectively. The stationary thermal elasticity problem is solved for these two setups. Finally, the boundary conditions are summarized in Fig. 3. For the mechanical field, the left and right parts are restricted to move in the \(x\)-direction and the upper and lower edges are fixed along \(y\)-direction. For the thermal field, constant temperatures are applied on the left and the right edges while on the upper and lower edges, there is no heat flux normal to the surface.

Figure 3: Boundary conditions for the thermo-mechanically coupled problem. Left: mechanical field boundary conditions. Right: Thermal field boundary conditions.

Figure 2: Selected geometries for studies on the thermal elasticity problem. The color variation shows the different values of Young’s modulus and thermal conductivity.
## 3 PINNs' architecture

### Thermo-elasticity problem

Based on the mixed formulation proposed by Rezaei et al. [1], the input for the network is the location of the collocation points (i.e. the coordinates \(x\) and \(y\) in a 2D setting). The outputs are then divided into the field variables from the mechanical and thermal sub-problem. For the mechanical field, outputs are the components of the displacement vector as well as the stress tensor (i.e. \(u_{x}\), \(u_{y}\), \(\sigma_{xx}\), \(\sigma_{xy}\), and \(\sigma_{yy}\)). For the thermal field, the output layer includes the scalar temperature value as well as the components of the heat flux (i.e. \(T\), \(q_{x}\), and \(q_{y}\)). It is important to mention that we intend to use the separate fully connected feed-forward neural network (FFNN) for each individual output variable (see also Fig. 4 and explanations in [1]). The structure of each neural network takes the standard form where it can be split into a single input layer, several possible hidden layers, and the output layer. Each layer is connected to the next layer for transferring information [52]. In every single layer, there is no connection among its neurons. Therefore, we represent the information bypass from the \(l-1\) layer to \(l\) via the vector \(\mathbf{z}^{l}\). Every component of vector \(\mathbf{z}^{l}\) is computed by

\[z_{m}^{l}=a(\sum_{n=1}^{N_{l}}w_{mn}^{l}z_{n}^{l-1}+b_{m}^{l}), \quad l=1,\ldots,L. \tag{18}\]

In Eq. (18), \(z_{n}^{l-1}\), is the \(n\)-th neuron within the \(l-1\)-th layer. The component \(w_{mn}\) shows the connection weight between the \(n\)-th neuron of the layer \(l-1\) and the \(m\)-th neuron of the layer \(l\). Every individual neuron in the \(l\)-th hidden layer owns a bias variable \(b_{m}^{l}\). The number \(N_{I}\) corresponds to the number of neurons in the \(l\)-th hidden layer. The total number of hidden layers is \(L\). The letter \(a\) stands for the activation function in each neuron. The activation function \(a(\cdot)\) is usually a non-linear function. In this work, the hyperbolic-tangent function is utilized which is defined as

\[\tanh(x)=\frac{\exp^{x}-\exp^{-x}}{\exp^{x}+\exp^{-x}}. \tag{19}\]

The proper choice of the activation function is problem dependent and shall be obtained based on hyperparameter studies on the final performance of the network. See also discussions in [1], [53] and [54].

The input layer of the PINNs for the described thermo-mechanical problem consists of the position vector \(x\) and \(y\) coordinates. For the output layer, we have \(u_{x}\), \(u_{y}\), \(\sigma_{xx}\), \(\sigma_{xy}\), and \(\sigma_{yy}\) in a mechanical field and \(T\), \(q_{x}\), and \(q_{y}\) in the thermal field. In the current work, the focus is on a 2-D setting, and extension to a 3-D setting is trivial. The rest of the network architecture will be discussed in what follows. The trainable set of parameters of the network is represented by \(\mathbf{\theta}=\{\mathbf{W},\mathbf{b}\}\), which are the weights and biases of a neural network, respectively. Considering each neural network structure as \(\mathcal{N}\), the final outcome of the network for the mechanical and thermal field reads

\[u_{x}=\mathcal{N}_{u_{x}}(\mathbf{X};\mathbf{\theta}),\ u_{y}=\mathcal{N}_{u_{y}}(\bm {X};\mathbf{\theta}),\ \sigma_{xx}^{o}=\mathcal{N}_{\sigma_{xx}}(\mathbf{X};\mathbf{\theta}),\ \sigma_{xy}^{o}=\mathcal{N}_{\sigma_{xy}}(\mathbf{X};\mathbf{\theta}),\ \sigma_{yy}^{o}=\mathcal{N}_{\sigma_{yy}}(\mathbf{X};\mathbf{\theta}), \tag{20}\]

\[T=\mathcal{N}_{T}(\mathbf{X};\mathbf{\theta}),\ q_{x}^{o}=\mathcal{N}_{q_{x}}(\mathbf{X}; \mathbf{\theta}),\ q_{y}^{o}=\mathcal{N}_{q_{x}}(\mathbf{X};\mathbf{\theta}). \tag{21}\]
The neural network outputs are functions of the trainable parameters and the training is done by minimizing different physical loss functions. Next, we build the BVP via the defined input and output layers. To do so, the partial differential equation and its corresponding initial and/or boundary conditions will be defined in terms of loss functions for the neural networks. Therefore, we require to obtain the derivatives of the output layer with respect to the input layer (for instance, the derivative of stress with respect to \(x\)- or \(y\)-direction). Current deep learning packages, i.e, Pytorch [55] and Tensorflow [56] are capable of computing the derivatives based on the automatic differentiation algorithms [57]. The algorithms developed in the current work are implemented using the SciANN package [58] but the methodology can easily be transferred to other platforms.

Denoting the summation of mechanical-related loss terms by \(\mathcal{L}_{M}\), it is defined based on a set of equations, Eqs. (2) - (8), as

\[\mathcal{L}_{M}=\underbrace{\mathcal{L}_{M}^{EF}+\mathcal{L}_{M}^{DBC}}_{\text {based on }\boldsymbol{u}}+\mathcal{L}_{M}^{cnc}+\underbrace{\mathcal{L}_{M}^{SF}+ \mathcal{L}_{M}^{NBC}}_{\text{based on }\boldsymbol{\sigma}^{o}}. \tag{22}\]

In Eq. (22), the energy form of the problem (\(\mathcal{L}_{M}^{EF}\)) and the prescribed Dirichlet boundary conditions losses (\(\mathcal{L}_{M}^{DBC}\)) are minimized by means of the displacement vector \(\boldsymbol{u}\) output. For Neumann boundary conditions (\(\mathcal{L}_{M}^{NBC}\)) and the strong form (\(\mathcal{L}_{M}^{SF}\)) of the problem requiring the first-order displacement derivatives, relevant loss terms are applied to the \(\boldsymbol{\sigma}^{o}\),

\[\mathcal{L}_{M}^{EF} =\text{MAE}_{M}^{EF}\left(\int_{\Omega}\frac{1}{2}\hat{\boldsymbol {e}}_{e}^{T}\,\hat{C}(x,y)\,\hat{\boldsymbol{e}}_{e}\,\ dV-\int_{\varGamma_{N}} \boldsymbol{u}^{T}\,\bar{\boldsymbol{t}}\;dA\right), \tag{23}\] \[\mathcal{L}_{M}^{DBC} =\text{MSE}_{M}^{DBC}\left(\boldsymbol{u}-\overline{\boldsymbol{ u}}\right),\] (24) \[\mathcal{L}_{M}^{cnc} =\text{MSE}_{M}^{cnc}\left(\boldsymbol{\sigma}^{o}-\boldsymbol{ \sigma}\right),\] (25) \[\mathcal{L}_{M}^{SF} =\text{MSE}_{M}^{SF}\left(\text{div}(\boldsymbol{\sigma}^{o}) \right),\] (26) \[\mathcal{L}_{M}^{NBC} =\text{MSE}_{M}^{NBC}\left(\boldsymbol{\sigma}^{o}\cdot \boldsymbol{n}-\overline{\boldsymbol{t}}\right). \tag{27}\]

The loss function related to the thermal field is also computed based on Eqs. (9)- (17).

Analogously to the mechanical losses, the thermal loss \(\mathcal{L}_{T}\) consists of \(\mathcal{L}_{T}^{EF}\) (energy form of the thermal diffusion problem), \(\mathcal{L}_{T}^{DBC}\) (Dirichlet boundary conditions), \(\mathcal{L}_{T}^{NBC}\) (Neumann boundary conditions), \(\mathcal{L}_{T}^{SF}\) (strong form of diffusion problem), and \(\mathcal{L}_{T}^{cnc}\) (connection term) and it reads

\[\mathcal{L}_{T} =\underbrace{\mathcal{L}_{T}^{EF}+\mathcal{L}_{T}^{DBC}}_{\text {based on }T}+\mathcal{L}_{T}^{cnc}+\underbrace{\mathcal{L}_{T}^{SF}+\mathcal{L}_{T}^{ NBC}}_{\text{based on }\boldsymbol{q}^{O}}, \tag{28}\] \[\mathcal{L}_{T}^{EF} =\text{MAE}_{T}^{EF}\left(\int_{\Omega}\frac{1}{2}\boldsymbol{q} ^{T}\cdot\nabla T\,dV+\int_{\bar{\varGamma}}(\boldsymbol{q}\cdot\boldsymbol{ n})\ T\,dA\right),\] (29) \[\mathcal{L}_{T}^{DBC} =\text{MSE}_{T}^{DBC}\left(T-\overline{T}\right),\] (30) \[\mathcal{L}_{T}^{cnc} =\text{MSE}_{T}^{cnc}\left(\boldsymbol{q}^{o}-\boldsymbol{q} \right),\] (31) \[\mathcal{L}_{T}^{SF} =\text{MSE}_{T}^{SF}\left(\text{div}(\boldsymbol{q}^{o})\right),\] (32) \[\mathcal{L}_{T}^{NBC} =\text{MSE}_{T}^{NBC}\left(\boldsymbol{q}^{o}\cdot\boldsymbol{n}- \overline{\boldsymbol{q}}\right). \tag{33}\]

**Remark 1** The energy form Eqs. (23) and (29) are evaluated at the global level and a single loss term is available for all of the collocation points. The other loss terms are minimized
at every single collocation point.

For completeness, the mean squared and mean absolute error are defined as

\[\text{MSE}(\bullet)_{type}=\frac{1}{k_{type}}\sum_{i=1}^{k_{type}}( \bullet)^{2},\quad\text{MAE}(\bullet)_{type}=\frac{1}{k_{type}}\sum_{i=1}^{k_{ type}}|\bullet|. \tag{34}\]

The mathematical optimization problem is written as

\[\mathbf{\theta}^{*}=\arg\min_{\mathbf{\theta}}\mathcal{L}(\mathbf{X};\mathbf{ \theta}), \tag{35}\]

where \(\mathbf{\theta}^{*}\) are the optimal trainable parameters (weights and biases) of the network which is the result of the optimization problem.

We investigate two distinct training algorithms to solve the coupled system of equations. In one approach, we utilize a so-called sequential training, and in the other approach, we intend to use a coupled training as depicted in Fig. 1. The sequential training algorithm is done by minimizing losses of one individual physics while the other physical parameters are kept frozen. Readers are also encouraged to see Haghighat et al. [59] for more details. Sequential training has a similar procedure as staggered numerical schemes [32; 36; 60]. For sequential training, we define the following loss terms

\[\mathcal{L}_{Seq1}=\mathcal{L}_{\text{Phy1}}=\mathcal{L}_{T}, \tag{36}\] \[\mathcal{L}_{Seq2}=\mathcal{L}_{\text{Phy2}}=\mathcal{L}_{M}. \tag{37}\]

During the minimization of the loss function in Eq. (36), the output related to the mechanical field (i.e. \(u_{x},u_{y},\sigma_{xx}^{o},\sigma_{xy}^{o},\text{ and }\sigma_{yy}^{o}\)) is kept frozen. While the loss function in Eq. (37) is minimized, the output related to the thermal field (i.e. \(T,q_{x}^{o},\text{ and }q_{y}^{o}\)) is kept frozen. On the other hand, we have coupled training where we minimize all the loss terms, simultaneously.

For the coupled training, the total loss function is the summation of all the loss terms and it is written as

\[\mathcal{L}_{Cpl}=\mathcal{L}_{\text{Phy1}}\,+\mathcal{L}_{\text{ Phy2}}=\mathcal{L}_{M}+\mathcal{L}_{T}. \tag{38}\]

Fig. 4 depicts the mixed PINN formulation for solving the thermoelasticity problem. In addition to the position vector (\(x\), \(y\)), the material parameters \(C(x,y)\) and \(k(x,y)\) are added to the input layer. The latter avails the capability of a network to learn the solution for other material properties.

In the proposed network structure, the primary variables and their spatial derivatives are directly evaluated. For example, when it comes to the mechanical problem, the displacement vector \(\mathbf{u}\) and the stress tensor \(\mathbf{\sigma^{o}}\) are defined as network output. Moreover, the strain tensor \(\mathbf{\varepsilon}\) is calculated by means of the predicted values for deformation (\(u_{x}\) and \(u_{y}\)), and the predicted temperature \(T\) according to Eqs. (1) and (2). By inserting the computed strain tensor (\(\mathbf{\varepsilon}\)) into Eq. (4), one computes the stress tensor \(\mathbf{\sigma}\) for the coupled system of equations. Consequently, the resulting stress tensor from the evaluated displacement \(\mathbf{\sigma}\) is linked to the network's direct evaluation \(\mathbf{\sigma}^{o}\) through the connection loss term \(\mathcal{L}_{M}^{enc}\) in Eq. (25). The latter is a soft constraint to connect \(\mathbf{\sigma}^{o}\) to \(\mathbf{\sigma}\) values.

The analogous strategy is employed for the thermal field in which the heat flux vector \(\mathbf{q}\) is derived on the basis of Eq. (12) where we compute the thermal field gradients by means of the automatic differentiation. The \(\mathbf{q}\) is connected to \(\mathbf{q}^{o}\) by the \(\mathcal{L}_{T}^{enc}\) described in Eq. (31).
According to Fig. 4, the optimization problem in a sequential training procedure starts with the first desired field loss function to enforce the minimization of the first physical field loss functions.

Figure 4: Network architecture and loss functions for the multi-physical PINNs, sequential and coupled training flowcharts for the thermoelasticity problem.
In the thermoelasticity problem, the training starts by finding the optimum network parameters to satisfy the thermal loss function (\(\mathcal{L}_{T}\)). This choice is made based on the fact that the thermal-induced strain tensor will affect the results of the mechanical field. In our proposed algorithm, the training is stopped for the thermal field loss function (first physics' loss function) after reaching a certain number of epochs \(n_{T}\). Subsequently, the training is done for the next physics' loss function (mechanical field loss function \(\mathcal{L}_{m}\)) for the same or a different number of epochs \(n_{M}\). The whole training will be stopped if the total number of epochs is reached the desired value \(n_{A}\). While the desired value for the total number of epochs has not been reached, the same procedure will be continued to reach that particular value (minimizing \(\mathcal{L}_{T}\) for \(n_{T}\) epochs and minimizing \(\mathcal{L}_{M}\) for \(n_{M}\) epochs). For the coupled training procedure, the total loss function is computed by summing the different physics loss functions (\(\mathcal{L}_{M}\) and \(\mathcal{L}_{T}\)) into a single one \(\mathcal{L}_{Cpl}\), see Eq. (38) and Fig. 4.

**Remark 2** It is shown that employing separated networks and utilizing only first-order derivatives in the network's loss function leads to a better performance [1]. See also studies by Chiu et al. [61] on this matter. Reducing the order of derivatives also has other benefits. For example, one has more flexibility in choosing the activation functions (i.e. we avoid problems related to vanishing gradients due to higher-order derivatives). Moreover, the computational cost is decreased compared to the cases where we have second-order derivation.

**Remark 3** Another possibility to stop each training procedure in the sequential training is to reach a certain desired value of the loss for each field. This is achieved by having an infinite number for the number of epochs while optimizing each field's loss function. The training procedure will stop whenever a certain value for the desired loss function is accomplished.

**Remark 4** The current mixed PINN approach can also be interpreted as a multi-objective optimization problem. Therefore, having a proper balance between different loss terms is essential. Often in multi-physical problems, it is important to normalize the quantities before the minimization starts [23]. Also, the loss terms \(L_{T}^{EF}\) and \(L_{M}^{EF}\) are defined based on the means of absolute error to be in the same order as other loss terms.

In this work, the Adam optimizer is employed, and the relevant network parameters are summarized in Table 1. Please note that in the case of a simple PINNs formulation, where the energy form of the problem is absent, one can use multiple batches for the training of the network.

\begin{table}
\begin{tabular}{l l} \hline \hline Parameter & Value \\ \hline Input, Output & \(\mathbf{X}\), \((\mathbf{u}\), \(\mathbf{\sigma}^{o}\), \(T\), \(\mathbf{q}^{o})\) \\ activation function & tanh \\ number of layers and neurons per layer (\(L\), \(N_{l}\)) & (5, 40) \\ batch size & full batch for the mixed PINNs formulation \\ (learning rate \(\alpha\), number of epochs) & (\(10^{-3},10^{5}\)) \\ \hline coupled training & \\ total number of epochs (\(n_{A}\)) & 100k epochs \\ \hline sequential training & \\ total number of iterations (\(n_{A}\)) & 200k epochs \\ iterations for the thermal field before switch (\(n_{T}\)) & 20k epochs \\ iterations for the mechanical field before switch (\(n_{M}\)) & 20k epochs \\ \hline \hline \end{tabular}
\end{table}
Table 1: Summary of the network parameters for coupled and sequential training procedures.
Based on the boundary conditions discussed in Fig. 3, we write all the loss terms more specifically for the given BVP in what follows. The components of the total strain tensor and temperature gradient vector are calculated by taking derivatives of displacements and temperature fields, respectively. Considering Eqs. (20) and (21), we define the following networks

\[\mathcal{N}_{\varepsilon_{x}}(\mathbf{X};\mathbf{\theta}) =\partial_{x}\mathcal{N}_{u_{x}}, \tag{39}\] \[\mathcal{N}_{\varepsilon_{y}}(\mathbf{X};\mathbf{\theta}) =\partial_{y}\mathcal{N}_{u_{y}},\] (40) \[\mathcal{N}_{\varepsilon_{xy}}(\mathbf{X};\mathbf{\theta}) =\frac{1}{2}\left(\partial_{y}\mathcal{N}_{u_{x}}+\partial_{x} \mathcal{N}_{u_{y}}\right),\] (41) \[\mathcal{N}_{\nabla_{x}T}(\mathbf{X};\mathbf{\theta}) =\partial_{x}\mathcal{N}_{T},\] (42) \[\mathcal{N}_{\nabla_{y}T}(\mathbf{X};\mathbf{\theta}) =\partial_{y}\mathcal{N}_{T}. \tag{43}\]

Components of the (isotropic) thermal strain tensor are then derived by

\[\mathcal{N}_{\varepsilon_{tx}}(\mathbf{X};\mathbf{\theta})=\mathcal{N}_{\varepsilon_{ ty}}=\alpha\left(\mathbf{X}\right)\left(\mathcal{N}_{T}-T_{0}\right). \tag{44}\]

Considering the total and thermally induced strain tensors, the elastic part of the strain tensor reads

\[\mathcal{N}_{\varepsilon_{ex}}(\mathbf{X};\mathbf{\theta}) =\mathcal{N}_{\varepsilon_{x}}-\mathcal{N}_{\varepsilon_{tx}}, \tag{45}\] \[\mathcal{N}_{\varepsilon_{ey}}(\mathbf{X};\mathbf{\theta}) =\mathcal{N}_{\varepsilon_{y}}-\mathcal{N}_{\varepsilon_{ty}},\] (46) \[\mathcal{N}_{\varepsilon_{exy}}(\mathbf{X};\mathbf{\theta}) =\mathcal{N}_{\varepsilon_{xy}}. \tag{47}\]

Utilizing the constitutive relations (Hooke's law) as well as Fourier's law, one computes the components of the stress tensor and heat flux vector according to the following relations

\[\mathcal{N}_{\sigma_{x}}(\mathbf{X};\mathbf{\theta}) =\frac{E(\mathbf{X})}{(1+\nu(\mathbf{X}))(1-2\nu(\mathbf{X}))}\left((1-\nu( \mathbf{X}))\mathcal{N}_{\varepsilon_{ex}}+\nu(\mathbf{X})\mathcal{N}_{\varepsilon_{ xy}}\right), \tag{48}\] \[\mathcal{N}_{\sigma_{y}}(\mathbf{X};\mathbf{\theta}) =\frac{E(\mathbf{X})}{(1+\nu(\mathbf{X}))(1-2\nu(\mathbf{X}))}\left(\nu(\mathbf{ X})\mathcal{N}_{\varepsilon_{x}}+(1-\nu(\mathbf{X}))\mathcal{N}_{\varepsilon_{y}} \right),\] (49) \[\mathcal{N}_{\sigma_{xy}}(\mathbf{X};\mathbf{\theta}) =\frac{E(\mathbf{X})}{2(1+\nu(\mathbf{X}))}\mathcal{N}_{\varepsilon_{exy}},\] (50) \[\mathcal{N}_{q_{x}}(\mathbf{X};\mathbf{\theta}) =-k(\mathbf{X})\mathcal{N}_{\nabla_{x}T},\] (51) \[\mathcal{N}_{q_{y}}(\mathbf{X};\mathbf{\theta}) =-k(\mathbf{X})\mathcal{N}_{\nabla_{y}T}. \tag{52}\]

To simplify equations for formulating loss functions, in what follows, the dependency of neural network output \(\mathcal{N}\left(\mathbf{X};\mathbf{\theta}\right)\) on the collocation points' coordinates \(\mathbf{X}\) and trainable parameters \(\mathbf{\theta}\) are not shown.

Finally, the mathematical expressions for calculating loss functions for the mechanical and thermal field are expanded in what follows. See also Eqs. (23) to (33) as well as Fig. 4,
\[\mathcal{L}_{M}^{EF} =\left|-\frac{1}{2N_{b}}\sum_{\mathbf{X}\in\Omega}\left(\mathcal{N}_{ \sigma_{x}}\mathcal{N}_{\varepsilon_{e_{x}}}+\mathcal{N}_{\sigma_{y}}\mathcal{N }_{\varepsilon_{e_{y}}}+\mathcal{N}_{\sigma_{xy}}\mathcal{N}_{\varepsilon_{e_{ xy}}}\right)\right|, \tag{53}\] \[\mathcal{L}_{M}^{DBC} =\frac{1}{N_{DBC}^{L}}\sum_{\mathbf{X}\in\Gamma_{D}^{L}}\left( \mathcal{N}_{u_{x}}-0\right)^{2}+\frac{1}{N_{DBC}^{R}}\sum_{\mathbf{X}\in\Gamma_{ B}^{B}}\left(\mathcal{N}_{u_{x}}-0\right)^{2},+\frac{1}{N_{DBC}^{T}}\sum_{\mathbf{X} \in\Gamma_{D}^{T}}\left(\mathcal{N}_{u_{y}}-0\right)^{2}\] \[+\frac{1}{N_{DBC}^{B}}\sum_{\mathbf{X}\in\Gamma_{B}^{B}}\left( \mathcal{N}_{u_{y}}-0\right)^{2},\] (54) \[\mathcal{L}_{M}^{cnc} =\frac{1}{N_{b}}\sum_{\mathbf{X}\in\Omega}\left(\left(\mathcal{N}_{ \sigma_{x}^{\circ}}-\mathcal{N}_{\sigma_{x}}\right)^{2}+\left(\mathcal{N}_{ \sigma_{y}^{\circ}}-\mathcal{N}_{\sigma_{y}}\right)^{2}+\left(\mathcal{N}_{ \sigma_{y}^{\circ}}-\mathcal{N}_{\sigma_{xy}}\right)^{2}\right),\] (55) \[\mathcal{L}_{M}^{SF} =\frac{1}{N_{b}}\sum_{\mathbf{X}\in\Omega}\left(\left(\partial_{x} \mathcal{N}_{\sigma_{x}^{\circ}}\,+\,\partial_{y}\mathcal{N}_{\sigma_{xy}^{ \circ}}\right)^{2}+\left(\partial_{y}\mathcal{N}_{\sigma_{y}^{\circ}}\,+\, \partial_{x}\mathcal{N}_{\sigma_{xy}^{\circ}}\right)^{2}\right),\] (56) \[\mathcal{L}_{M}^{NBC} =0. \tag{57}\]

In the above equations, \(\Gamma^{L}\), \(\Gamma^{R}\), \(\Gamma^{T}\), and \(\Gamma^{B}\) denote the location of points at left, right, top, and bottom boundaries, respectively. _DBC_ and _NBC_ represent for Dirichlet and Neumann boundary conditions.

Note, that the loss term related to the Neumann BCs (Eq. (57)) is zero according to the defined BVP. We do not apply any external traction to the microstructure of the material. For the computation of the integral in the energy form, we use the summation over the collocation points with equal weights as the collocation points are uniformly distributed in the domain. Therefore, \(N_{b}\) denotes the number of collocation points within the domain (body). Moreover, \(N^{L}\), \(N^{R}\), \(N^{T}\), and \(N^{B}\) represent the number of collocation points at left, right, top, and bottom boundaries, respectively. See also [31] for other possible methodologies.

The loss functions related to the thermal field read as

\[\mathcal{L}_{T}^{EF} =\left|\frac{1}{2N_{b}}\sum_{\mathbf{X}\in\Omega}\left(\mathcal{N}_{ q_{x}}\mathcal{N}_{\nabla_{x}T}+\mathcal{N}_{q_{y}}\mathcal{N}_{\nabla_{y}T} \right)+\frac{1}{N_{NBC}^{L}}\sum_{\mathbf{X}\in\Gamma_{N}^{L}}\left(\mathcal{N}_ {q_{x}}\mathcal{N}_{T}\right)\right|, \tag{58}\] \[\mathcal{L}_{T}^{DBC} =\frac{1}{N_{DBC}^{L}}\sum_{\mathbf{X}\in\Gamma_{D}^{L}}\left( \mathcal{N}_{T}-1\right)^{2}+\frac{1}{N_{DBC}^{R}}\sum_{\mathbf{X}\in\Gamma_{D}^{ B}}\left(\mathcal{N}_{T}-0\right)^{2},\] (59) \[\mathcal{L}_{T}^{enc} =\frac{1}{N_{b}}\sum_{\mathbf{X}\in\Omega}\left(\left(\mathcal{N}_{ q_{x}^{\circ}}-\mathcal{N}_{q_{x}}\right)^{2}+\left(\mathcal{N}_{q_{y}^{\circ}}- \mathcal{N}_{q_{y}}\right)^{2}\right),\] (60) \[\mathcal{L}_{T}^{SF} =\frac{1}{N_{b}}\sum_{\mathbf{X}\in\Omega}\left(\partial_{x}\mathcal{ N}_{q_{x}^{\circ}}\,+\,\partial_{y}\mathcal{N}_{q_{y}^{\circ}}\right)^{2},\] (61) \[\mathcal{L}_{T}^{NBC} =\sum_{\mathbf{X}\in\Gamma_{N}^{T}}\frac{1}{N_{NBC}^{T}}\left(\mathcal{ N}_{q_{y}^{\circ}}-0\right)^{2}\] \[+\sum_{\mathbf{X}\in\Gamma_{N}^{B}}\frac{1}{N_{NBC}^{B}}\left( \mathcal{N}_{q_{y}^{\circ}}-0\right)^{2}. \tag{62}\]
## 4 Results

### A circular heterogeneity in a matrix (geometry 1)

We start with the first geometry which is described in Fig. 2. The collocation points are set as depicted on the left-hand side of Fig. 5. As explained in [1], to enhance the accuracy of the results, one needs to increase the collocation points in the heterogeneity, domain boundaries, and matrix. Approximately, 5000 collocation points are considered in a rectangular domain from 0 to 1 [mm] in both \(x\)- and \(y\)-directions. Minimization of the loss functions (described in Sec. 3.1) is done at the collocations points. The material parameters for the current boundary value problem are mentioned in Tab. 2.

To evaluate the network response, the thermoelasticity problem is also solved for the same domain by the finite element program FEAP [62]. The discretization of the domain for FE analysis is done in a way that the elements coincide with the regular collocation points for evaluating the response of the trained network. Moreover, in Fig. 5, the variation of color shows the different Young's moduli and thermal conductivity values in the PINNs collocation points. The latter comes from the fact that some points are living at the interface of the matrix and the inclusion. The same procedure exists during the assembly procedure in the FE method.

The minimization of the loss functions is done by two different procedures, coupled and sequential training (see also Fig. 4). The training is performed on the Skylake CPU, Platinum 8160 model, and Intel HNS2600BPB Hardware Node Type for both cases. For the coupled training, every loss term is minimized by performing 100 k epochs (iterations). To have a fair

\begin{table}
\begin{tabular}{l l} \hline  & Value/Unit \\ \hline \hline
**Thermal elasticity problem of geometry 1** & \\ Matrix Young’s modulus and Poisson’s ratio (\(E_{\text{mat}}\), \(\nu_{\text{mat}}\)) & (0.3 GPa, 0.3) \\ Inclusion Young’s modulus and Poisson’s ratio (\(E_{\text{inc}}\), \(\nu_{\text{inc}}\)) & (1.0 GPa, 0.3) \\ Matrix heat conductivity (\(k_{\text{mat}}\)) & 0.3 W/mK \\ Inclusion heat conductivity (\(k_{\text{inc}}\)) & 1.0 W/mK \\ \hline \end{tabular}
\end{table}
Table 2: Geometry 1 input parameters for the linear thermal elasticity problem

Figure 5: Left: collocation points for minimizing the loss functions in PINNs. Middle: points for evaluating the trained network response. Right: discretized domain using FEM to solve the same boundary value problem in thermoelasticity. The color variation shows the different values of Young’s modulus.
comparison for the sequential training, we perform a total of 200 k epochs. In the sequential training, we minimize each loss term for 100 k epochs. The number of epochs is chosen based on the final loss value and accuracy of the results. This number may easily change based on the optimizer type and the preference for the desired accuracy.

#### 4.1.1 Coupled training for the geometry 1

The decay of prominent loss terms during the training procedure is shown in Fig. 6.

In the coupled training minimization, we try to find the optimum network parameters to satisfy all of the loss terms, simultaneously. Note that for the sake of clarity and a better

Figure 6: Prominent loss terms for the thermal and mechanical field decay during the training procedure utilizing the coupled training algorithm.
comparison with the sequential training, the loss terms related to the mechanical and thermal sub-problem are plotted individually in Fig. 6.

The deformed configurations of PINNs and FEM are plotted in Fig. 7. The deformed configurations are almost identical. Note that all the deformations are induced due to the introduced temperature gradient in the system. To have a more detailed view, the comparison is made in Fig. 8 for the displacements and the stress components in PINNs and FEM results. The PINNs results are presented in the first column while the FEM results and the difference between PINNs and FEM are shown in the second and third columns, respectively.

According to the results in Fig. 8, the mechanical field results from PINNs for the thermoelasticity problem are in good agreement with the outcomes from FE analysis. The averaged relative difference lies at \(0.59\,\%\) for displacements in \(x\)-direction and is about \(0.87\,\%\) for stresses in \(x\)-direction. The maximum relative difference for the displacements and stress components are about \(4.1\,\%\), and \(13.5\,\%\), for \(u_{x}\) and \(\sigma_{x}\), respectively. The previously measured errors only exist at a few points in the domain around the inclusion boundaries.

The maximum amount of difference is observed at the interface of the matrix and heterogeneity. This point is the result of having a rather sharp transition in the material properties at the interface. One possible remedy to overcome the latter problem is to add more collocation points near the interface or by applying domain decomposition methods, see [23, 63, 64]. Increasing the number of epochs for training is also helpful.

The agreement of the thermal field results between PINN and FE analysis in the thermoelasticity problem for the coupled training is also addressed in Fig. 9. The relative averaged error for \(T\) is about \(0.3\,\%\) and for \(q_{x}\) is around \(0.76\,\%\). The maximum relative difference between PINN and FE analysis for the thermal field is about \(1.0\,\%\) and \(24.93\,\%\) for the temperature field and the heat flux in the \(x\)-direction, respectively. The measured difference again occurs mostly at a few points on the interface. In the next section, a more detailed comparison is made by applying sections in \(x\)-direction and having a closer look at the response of PINN and FE analysis.

Figure 7: Left: undeformed configuration of evaluation points in PINNs. Middle: deformed configuration of evaluation points in PINNs employing coupled training. Right: deformed configuration of the discretized domain in FE analysis.
Figure 8: The mechanical field results of the thermoelasticity problem employing the coupled training procedure. Left: PINNs, Middle: FEM, Right: FEM with PINNs difference.
#### 4.1.2 Sequential training

The identical collocation points from Fig. 5 are utilized to minimize the thermal and mechanical loss functions (\(\mathcal{L}_{T}\) and \(\mathcal{L}_{M}\)) by employing the sequential algorithm. As it is discussed beforehand, the minimization is performed in total for \(200\,\mathrm{k}\) epochs which results in \(100\,\mathrm{k}\) for each field loss function. The decay of main losses for each field with respect to the number is plotted separately for the mechanical and thermal fields in Fig. 10.

In Fig. 10, one can see the effect of chosen values for \(n_{T}\) and \(n_{M}\). These two numbers determine when we switch between the minimization of thermal and mechanical loss functions. The value for these two parameters is set to \(20\,\mathrm{k}\). Take into consideration that the loss terms for the thermal field are not evaluated during the mechanical part training, and vice versa. Therefore, to show the whole loss functions the losses are saved for each individual part during its training and later on appended together. One interprets that by minimization of the one field loss functions, other fields' losses may increase. The network tries to find the optimum parameters for the active field that may change the previously found optimum parameters for the nonactive fields. The latter point is hardly seen for each loss function right after upright lines in Fig. 10, specifically for the thermal field losses. Finally, one needs to state that having some iterations between the fields during training is essential to find the global minimum of the system. For this example, \(n_{M}\), \(n_{T}\), and \(n_{A}^{Seq}\) are chosen in a way to make five iterations possible.

Figure 9: The thermal field results of the thermoelasticity problem employing the coupled training procedure. Left: PINNs, Middle: FEM, Right: FEM with PINNs difference.
Further studies are required to find the optimum switching times between the fields and the number of epochs while minimizing one field loss function.

Utilizing a mixed PINN by employing sequential training leads to \(0.81\,\mathrm{\char 37}\) and \(0.8\,\mathrm{\char 37}\) averaged relative difference for the displacement and stress components in \(x\)-direction with FE analysis, respectively, while the maximum relative errors are \(5.9\,\mathrm{\char 37}\), and \(14.89\,\mathrm{\char 37}\). Again, one should emphasize that the error is localized in a few points around the inclusion. In the thermal field, the maximum relative temperature error stands at around \(1.3\,\mathrm{\char 37}\) and at almost \(29.16\,\mathrm{\char 37}\) for the heat flux in \(y\)-direction, \(q_{y}\), see Fig. 15. The averaged relative errors lie at \(0.26\,\mathrm{\char 37}\), and \(0.66\,\mathrm{\char 37}\)

Figure 10: Top: prominent loss terms for the thermal field decay during the training procedure utilizing the sequential training algorithm. Bottom: prominent loss terms for the mechanical field decay during the training procedure utilizing the sequential training algorithm.
for the mentioned values.

To emphasize the capability of the mixed PINNs formulation in comparison with the standard PINNs structure [20; 37], the standard PINN is used to solve the thermoelasticity boundary value problem with an identical number of epochs for training. Note that the standard PINNs structure utilizes a separate neural network to predict output which are only \(u\) and \(T\). Therefore, second-order derivatives appear in the formulation. The network has five hidden layers and forty neurons in each of them (\(5\times 40\)) with having the _tanh_ as an activation function. The training is done by both coupled and sequential training algorithms for each of the normal PINNs and mixed PINNs formulations.

The \(x\)-components of displacement, stress, and heat flux as well as the temperature field are depicted in Fig. 11 for the section which is made at \(x=0.5\,\mathrm{mm}\). The FE analysis results are also added as a reference solution to the figures.

The standard PINNs formulation is not able to show the effect of heterogeneity for the displacement component \(u_{x}\), and for \(T\), the temperature field. It also shows poor performance for

Figure 11: Comparison of main output of the networks through the section at \(x=0.5\mathrm{mm}\). Mix. PINNs: denotes the formulation of this paper by having the \(\mathbf{\sigma}\) and \(\mathbf{u}\) as output and utilizing the energy form loss in addition to conventional losses. PINNs: stands for the standard same formulation which has only \(\mathbf{u}\) as output in the absence of the energy form of the problem in the loss function. Coupled: minimizing loss function utilizing coupled training algorithm (all fields loss function are minimized simultaneously). Sequential: minimization is done with the proposed sequential training method.
\(q_{x}\), and \(\sigma_{x}\) prediction, regardless of having coupled or sequential training procedures. However, employing sequential training leads to slightly better performance for stress and flux components. The mixed PINNs formulation shows much better performance and is able to capture the effect of heterogeneity for both coupled and sequential training algorithms. By having a closer look at the results, one sees the modestly better performance of coupled training algorithm compared to the sequential one. Due to having a pixel-based mesh in FEM discretization, some small jumps are exited for FE analysis results in Fig. 11 which does not exist for the mixed PINNs formulation. One can state that the mixed PINNs perfectly solve the current multi-physical problem by utilizing sequential or coupled training procedures and lead to the almost same results as the FE analysis.

**Remark 5** To address the computational efficiency of each training procedure, the average value of one epoch for coupled training is about 0.67 seconds while for the sequential training, this value decreased to 0.16 seconds. However, since one needs to perform at least two iterations to achieve the solution in the case of sequential training, the computational gain is about 0.35 seconds which reduces the computational time more than two times compared to the coupled training method.

### Complex microstructure containing rectangular inclusions (geometry 2)

To further evaluate the performance of the proposed approach for solving multi-physical in the heterogeneous domains, the second geometry in Fig. 12 is now under consideration. The pattern of heterogeneity is defined to be more challenging with sharp corners. Such a pattern also corresponds to some specific material microstructures such as nickel alloy systems, see [65].

For the second geometry, we introduce 2601 collocation points to minimize the loss functions (see the left-hand side of Fig. 12). The corresponding mesh, which discretizes the domain in order to solve the thermoelasticity problem by means of FEM, is depicted on the right-hand side of Fig. 12. To evaluate the response of the system, 10201 points are employed. In this case, approximately 75 % of points have different positions concerning the collocation points, see Fig. 12 middle part. The latter addresses the applicability of the network to interpolate the solution. The material parameters utilized in the second boundary value problem are stated in Tab. 3.

Figure 12: Left: collocation points for minimizing the loss functions for the second geometry. Middle: points for evaluating the trained network response. Right: discretized domain using FEM to solve the same boundary value problem in thermoelasticity. The color variation shows the different values of Young’s modulus.
#### 4.2.1 Coupled training of the second geometry

The loss functions according to the thermoelasticity problem are minimized by the coupled training procedure for the second geometry. Fig. 13 illustrates the deformed (magnification is one time) and undeformed configurations for the microstructure utilizing PINNs and FEM. The PINNs mixed formulation is capable of predicting the response of this multi-physical problem, despite the shape of heterogeneity and also the large deformation of inclusions. However, on the left-hand side of the deformed configuration in Fig. 13, one can observe that the left boundary condition is not satisfied. The latter is the main source of error for the mechanical field responses of the network. See the first row in Fig. 14, where the left Dirichlet boundary condition is not fully satisfied in certain regions.

A remedy for tackling this problem is using hard constraint boundary conditions or adding extra collocation points at boundaries. Similar to the results of the first geometry, the difference is accumulated mainly at interface regions. The maximum differences between the mixed PINNs approach and the FEM are \(9.7\,\mathrm{\char 37}\) and \(13.4\,\mathrm{\char 37}\) for the \(u_{x}\) and \(\sigma_{x}\), respectively. The averaged relative errors for the mentioned fields are \(1.6\,\mathrm{\char 37}\), and \(3.3\,\mathrm{\char 37}\) which shows again the localization of the error around inclusions. The maximum relative difference between the trained network and the FEM lies at around \(3.4\,\mathrm{\char 37}\) for the temperature field and at \(21.5\,\mathrm{\char 37}\) for \(q_{y}\) which is the critical component of the flux vector. For the temperature prediction, the error lives near the left boundary while for the fluxes right boundary is also critical, see Fig. 15. For the mentioned values, the averaged relative errors are \(1.3\,\mathrm{\char 37}\), \(1.8\,\mathrm{\char 37}\).

Figure 13: Undeformed and deformed configurations result from the mixed PINNs formulation employing coupled training algorithm and FEM respecting the second geometry. No scaling is applied.

\begin{table}
\begin{tabular}{l l} \hline \hline  & Value/Unit \\ \hline \hline
**Thermal elasticity problem of geometry 2** & \\ Matrix Young’s modulus and Poisson’s ratio (\(E_{\mathrm{mat}}\), \(\nu_{\mathrm{mat}}\)) & (1.0 GPa, 0.3) \\ Inclusion Young’s modulus and Poisson’s ratio (\(E_{\mathrm{inc}}\), \(\nu_{\mathrm{inc}}\)) & (0.5 GPa, 0.3) \\ Matrix heat conductivity (\(\mathbf{k}_{\mathrm{mat}}\)) & 1.0 W/mK \\ Inclusion heat conductivity (\(\mathbf{k}_{\mathrm{inc}}\)) & 0.5 W/mK \\ \hline \hline \end{tabular}
\end{table}
Table 3: Geometry 2 input parameters for the linear thermal elasticity problem
Figure 14: The mechanical field results of the thermoelasticity for geometry 2 using coupled training. Left: PINNs, Middle: FEM, Right: FEM with PINNs difference.
#### 4.2.2 Sequential training of the second geometry

The same number of collocation points as in Sec. 4.2.1 are utilized to minimize the loss function in the sequential training approach. The corresponding results of the network are partially plotted in Fig. 16.

We also observed that the same issue occurs also for this training strategy and the Dirichlet boundary condition of the left-hand side is not fully satisfied. In general, the results are in good agreement with the results from FEM, the averaged relative difference is about \(0.49\,\%\), and \(2.12\,\%\) which are for \(u_{x}\) and \(\sigma_{x}\), respectively. The maximum calculated relative errors for the mentioned fields are \(2.8\,\%\), and \(10.19\,\%\). For the thermal field, relative maximum and relative averaged errors lie at \(2.5\,\%\), and \(0.7\,\%\). These values for the heat flux in \(y\)-direction are \(21.5\,\%\), and \(1.82\,\%\).

For carrying out a closer inspection of the results employing coupled and sequential methods, the components of the mechanical field outputs and heat flux in \(x\)-direction as well as the temperature field are compared by drawing a section at \(x=0.5\)mm. The mentioned sections are made in such a way that at least they go through one inclusion part. The results of PINNs and coupled training are compared by having FE analysis results as the reference solution in Figs. 16. The results of sequential training slightly outperform the coupled training results with respect to the FEM. More specifically, near the inclusion, the results of the sequential

Figure 15: The thermal field results of the thermoelasticity for the second geometry using coupled training. Left: PINNs, Middle: FEM, Right: FEM with PINNs difference.
training are closer to the reference solution. The related relative difference values mimic the outperformance of the sequential training strategy.

### Transfer learning by combining physics and data

In this section, we aim to obtain full-field solutions for unseen cases by combining data and physics. The trained network reaches the solution of the given BVP promptly. Therefore, it makes the neural network enormously faster than a conventional solver like FE.

For the given BVP, one can change the topology of the microstructure, boundary conditions as well as material properties. For the sake of simplicity, we focus only on different material properties for different involved phases. Therefore, geometry 2 from the previous section is taken with the set of different material parameters. To this end, the ratio of Young's modulus between inclusion and the matrix as well as their thermal conductivity ratio is changing in the range of 1 to 10. For each case, the thermoelasticity problem is solved by means of FE analysis, and the corresponding results, later on, be used to train the neural network. The FE data is utilized to facilitate the process of training and reduce computational costs.

For training the network, one needs to insert Young's modulus and thermal conductivity of each point alongside the collocation points' coordinates as inputs. The designed network's architecture is depicted in Fig. 17. In addition to the loss based on data, the losses based on

Figure 16: Comparison of the network’s outputs through the section at \(x=0.5\)mm of the geometry 2.
the mixed PINNs formulation are taken into account, \(\mathcal{L}_{T}\) and \(\mathcal{L}_{M}\) see Eqs. (36) and (37). The total loss function in the presence of data from FE analysis is written as

\[\mathcal{L}_{total}=\underbrace{\mathcal{L}_{data}}_{\text{FE simulations}}+\,w\left(\underbrace{\mathcal{L}_{T}+\mathcal{L}_{M}}_{\text{Physics}}\right). \tag{63}\]

In Eq. (63), \(w\) denotes the weight of the physics in the total loss function. The network's predictions are compared for four different cases in which the \(w\) parameter is varying from \(0.0\) to \(1.0\). In the case of \(w=0.0\), the data-driven case is investigated. After studying a set of different values for the weight value, the value of \(w=0.1\) is taken that reads the best results. The training is done by employing the transfer learning procedure where the \(\mathcal{L}_{total}\) in Eq. (63) is minimized considering a specific ratio for Young's modulus \(E(x,y)\) and thermal conductivity \(k(x,y)\) of different phases. Next, the network's optimum parameters are found for the first set of values after \(10\,\mathrm{k}\) epochs. By having the previously tuned parameters, the minimization of the losses is done for the new value of the ratio for \(10\,\mathrm{k}\) iterations.
To have a better comparison, a section is made at \(x=0.35\) mm. The results from the combination of data and physics as well as the case with the pure data-driven case are reported. The interest outputs in these sections are \(u_{x}\), \(\sigma_{x}\) for the mechanical field. From the thermal field, \(T\), and \(q_{x}\) are chosen for comparison.

In general, corresponding results from the mentioned approaches are in good agreement with the reference solution which is computed by the FE analysis, see Fig. 20. The boundary values are approximated better by the mixed PINNs formulation combined with data. Both methods are able to observe the effects of inclusions and the use of physical constraints enhances the accuracy of predictions in areas where heterogeneity exists.

Figure 18: Mechanical and thermal fields results of the thermoelasticity for the case of \(E_{r}=K_{r}=15\). Left: PINNs + Data (\(w=0.1\)), Middle: FEM, Right: FEM with PINNs difference.
**Remark 6** The computational cost of the data-driven network employing transfer learning is almost \(82\,\%\) less than the network which is trained by combining data and physics. For constructing physical constraints, one needs to differentiate the output parameters with respect to inputs which is the main reason for the higher computational cost.

Figure 19: Mechanical and thermal fields results of the thermoelasticity for the case of \(E_{r}=K_{r}=15\). Left: Data-driven (\(w=0.0\)), Middle: FEM, Right: FEM with a data-driven difference.
## 5 Conclusion and outlooks

This work shows the capability of physics-informed neural networks in the absence of any ground truth data (unsupervised learning) to solve a multi-physical problem (coupling PDEs) in heterogeneous domains. The latter is achieved by utilizing the mixed PINNs formulation proposed by Rezaei et al. [1], where the authors used only first-order derivatives to formulate the loss functions which leads to a better approximation of the solution of a boundary value problem compared to standard field approaches where second-order derivatives are existing in the loss functions. It is shown that employing sequential training (minimization of the loss functions is done step-wise: first one field is minimized by freezing the other fields and in the same way losses of other fields are minimized) is beneficial in terms of computational cost and also the accuracy of the network. The proposed network is employed to solve a quasi-static thermoelasticity problem. The corresponding results are compared to the results of the coupled training where all losses from different fields are minimized simultaneously.

The prescribed network architecture is employed in combination with data (in order to reduce the computational cost) to predict unseen cases. To do so and by having Young's modulus and thermal conductivity ratios as additional input parameters, the training is done for several values

Figure 20: Comparison of main outputs of the networks through the section at \(x=0.5\)mm of the second geometry using mixed PINNs formulation combined with data utilizing transfer learning training procedure or pure data-driven scheme.
of ratios, and the resulting network is used to predict (extrapolate) an unseen value of the ratio. The trained network can be used to predict the solution of a given boundary value problem with almost no computational cost. The latter can be utilized in digital twins, where the solution of the system is required in a fraction of a second. Further studies are needed to predict the response of arbitrary microstructures under the same boundary and initial conditions.

In future works, one could aim to solve other coupled partial differential equations such as the phase-field damage model [36, 48]. Furthermore, one can combine more complex physics for predicting cracking in multi-physical environments [32, 66]. As an example, the trained networks can be utilized for structural optimizations and optimal design of materials created through additive manufacturing. The methodology discussed in this paper for transfer learning can be employed to generalize the network for various problem configurations. It is also rewarding to combine and compare these ideas with those from operator learning [3, 67].

**Acknowledgements**:

The authors acknowledge the computing time granted via the high-performance computers of RWTH Aachen University. Ali Harandi and Stefanie Reese would like to acknowledge the financial support of the German Science Foundation (DFG) for project A6 "Multi-scale modeling of damage and fracture behavior of nano-structured layers" as project number 138690629.

**Author Statement**:

Ali Harandi: Methodology, Software, Writing - Review & Editing. Ahmad Moidendin: Methodology, Software, Writing - Review & Editing. Michael Kaliske: Supervision, Review & Editing. Stefanie Reese: Funding acquisition, Review & Editing. Shahed Rezaei: Conceptualization, Methodology, Supervision, Writing - Review & Editing.
