Infrared stress measurement is a noncontact imaging technique that enables the visualization of the temperature distribution on an object’s surface by measuring the infrared radiation emitted from it using an infrared sensor. When mechanical stress is applied to a material, a slight temperature change, known as the thermoelastic effect is observed. This phenomenon enables the estimation of variations in principal stress sum in a nondestructive and noncontact manner.

This effect is theoretically described by Kelvin’s equation as shown in Equation 1. Δ T = − k T Δ σ sum , (1) where Δ T denotes the temperature change, T is the absolute temperature, Δ σ sum is the change in the sum of principal stresses, and k is the thermoelastic coefficient, which is given by Equation 2. k = α ρ C p . (2)

In this expression, α represents the coefficient of thermal expansion, ρ is the material density, and C p is the specific heat at a constant pressure.

The infrared stress measurement based on this theoretical framework has been successfully applied to not only metallic materials but also CFRP laminates.

Principal stresses are the eigenvalues obtained by diagonalizing the stress tensor at a given point within a material. They represent the normal stresses acting on mutually orthogonal planes where shear stresses vanish. In a Cartesian coordinate system, the stress tensor σ can be diagonalized such that the diagonal components σ 1 , σ 2 and σ 3 correspond to the principal stresses.

The sum of principal stresses, referred to in this paper as DSPSS, is defined as the trace of the stress tensor, that is, the sum of its diagonal components. This can be expressed in two equivalent forms, namely, by Equation 3 and, σ sum = σ 1 + σ 2 + σ 3 (3) or equivalently, using the stress tensor in Equation 4 and the resulting trace in Equation 5: σ = σ x x τ x y τ x z τ y x σ y y τ y z τ z x τ z y σ z z , (4) σ sum = tr σ = σ x x + σ y y + σ z z . (5)

This scalar value provides a comprehensive measure of the overall mechanical stress intensity on the surface.

GNN (Scarselli et al., 2009) belongs to a class of deep learning models specifically designed for data with graph structures. Unlike conventional neural networks, which are optimized for regular structures such as images and sequences, GNN operates directly on graphs composed of nodes and edges.

GNN updates node and edge features by leveraging the graph's structure, enabling the learning of a holistic graph representation. This makes GNN particularly well suited for utilizing the mesh topology obtained from FEM simulations.

The fundamental mechanism of GNN is message passing (Gilmer et al., 2017), which updates each node's feature vector by aggregating messages from neighboring nodes. The general update process at the l th layer is described by the message function in Equation 6, the aggregation in Equation 7, and the update function in Equation 8. m u → v l = f msg h u l − 1 , h v l − 1 , e u v , (6) where m u → v ( l ) denotes the message from node u to node v , h u ( l − 1 ) and h v ( l − 1 ) are the feature vectors at the ( l − 1 ) th layer, e u v is the edge feature, and f msg is the message function.

Aggregated messages for node v : m v l = ∑ u ∈ N v m u → v l . (7) The updated feature vector for node v is h v l = f upd h v l − 1 , m v l , (8) where f upd denotes the update function and N ( v ) represents the set of neighbors of node v .

In this study, we use a specific GNN architecture called the graph attention network (GAT) (VeliÄkoviÄ et al., 2018), which introduces attention mechanisms to learn the importance of neighboring nodes. Each neighboring node is assigned a learnable weight, allowing the model to focus more on relevant neighbors during feature aggregation.

The basic GAT update for node v at the ( l + 1 ) th layer is given by Equation 9. h l + 1 v = σ ∑ u ∈ N v α v u W l + 1 h u l , (9) where σ is the activation function, W ( l + 1 ) the weight matrix at the ( l + 1 ) th layer, and α v u the attention coefficient between nodes v and u , computed using Equation 10: α v u = e ϕ a ⊤ W l + 1 h v l ‖ W l + 1 h u l ∑ k ∈ N v e ϕ a ⊤ W l + 1 h v l ‖ W l + 1 h k l , (10) ϕ x = x if   x ≥ 0 λ x if   x < 0 ,where   λ ∈ 0,1 . (11) Where ϕ ( x ) denotes the activation function LeakyReLU, which is defined in Equation 11 as a piecewise linear function with a small slope λ for negative inputs. The vector a is a learnable attention weight vector and ‖ denotes the vector concatenation. This mechanism enhances the model's ability to focus on influential neighboring nodes during updates.

The analysis conditions for the CFRP space vehicle structure modeled by FEM are illustrated in Figure 1. The target of the analysis is a scaled-down model representing part of a cylindrical curved interstage structure made of CFRP, similar to those used in the H-IIA rocket. The original structure is a large curved panel with a diameter of approximately 4.0  m , a longitudinal length of about 7.0  m , and an arc length of 12.6  m . This structure is scaled down by a factor of 1/80 with the curvature and geometric characteristics maintained, and the resulting CFRP curved panel is used as the analysis target (Ura et al., 1998).

The dimensions of the curved panel are approximately 2000.0  mm in radius, 1570.0  mm in arc length, 700.0  mm in vertical length, and 22.0  mm in thickness. Square and rectangular holes are introduced at the center to simulate openings typically found in space launch vehicles: 100.0  mm × 100.0  mm for square holes and 200.0  mm × 100.0  mm for rectangular holes. The total thickness of the panel is 22.0  mm , consisting of a 3.0  mm –thick CFRP laminate on the top, a 16.0  mm –thick foam core in the middle, and another 3.0  mm –thick CFRP laminate at the bottom. This sandwich structure design ensures high stiffness while maintaining lightweight structure. The core material used for the CFRP-foam core sandwich structure is Rohacell 110WF, a polymethacrylimide (PMI) rigid foam manufactured by Evonik (Kobayashi, 2023). This material is widely used in aerospace applications owing to its high specific strength and stiffness, and stable mechanical properties even under cryogenic conditions.

The CFRP layers are composed of ten plies of unidirectional prepreg ( 0.3  mm –thick each) laminated on both sides of the foam core (Shimazaki et al., 2015). The fiber orientations of the stacked unidirectional composites are 0 ° , 45 ° , 90 ° , − 45 ° , 0 ° , 0 ° , − 45 ° , 90 ° , 45 ° , and 0 ° . The z -axis in Figure 1 corresponds to the fiber direction of 0 ° .

Defects are inserted away from the red box area, which is expected to be significantly affected by stress concentration around the hole under tensile loading. As shown in the white boxes in Figure 2, each defect measures 100.0  mm × 100.0  mm × 0.3  mm . Defects are implemented in the FEM model by modifying material properties to elements corresponding to the defect regions. Defects are inserted into the 18 internal plies excluding the top and bottom CFRP laminates, i.e., 1st and 20th layers and the foam core. The material properties of the CFRP, foam core, and defect regions are summarized in Table 1. To simulate interlaminar delamination, which is commonly performed in experiments by inserting Teflon sheets between prepregs layers (LÃpezâ et al., 2010), the material properties of Teflon sheets used in a previous study are referenced for the defect region (Kojima et al., 2022).

The mesh size is set to 12.5  mm . As boundary conditions, periodic boundary conditions are applied to both ends of the z -axis ( ∂ X [ + 1 ] and ∂ X [ − 1 ] ). The left edge in the x -axis direction is fully constrained, whereas a uniform displacement boundary condition is applied at the right edge, with an imposed displacement of 10.0  mm in the z -axis direction.

Under these conditions, FEM simulations are conducted to generate paired datasets consisting of the three-dimensional location of defects and the corresponding DSPSS on the curved panel. In total, one dataset without defects and 1,386 datasets with defects are prepared.

Figure 3 illustrates the proposed inverse defect localization framework. As shown in Figure 3a, for GNN training, we use only the DSPSS of the two outermost surface layers. As shown in Figure 3b, this method involves training a GNN using the normalized DSPSS obtained from FEM simulations to predict the three-dimensional location of defects.

To construct the training dataset, defects are inserted into the FEM model, and labels are assigned to the nodes corresponding to their locations, as shown in Figure 4. In this study, the 1st layer (bottommost) and 20th layer (topmost) of the CFRP structure are referred to as the bottom and upper surfaces respectively. Let V be the set of nodes in the mesh graph, and let y i ∈ { 0,1 , … , 18 } denote the class label assigned to node v i ∈ V . The class label corresponds to the defect insertion layer index minus 1, as defined in Equation 12: y i = c * − 1 if node v i lies on the upper surface of a defect inserted in layer c * 0 otherwise . (12)

Using this labeling, we formulate the training objective as a 19-class node classification problem.

The group index for a given class c is defined using the following grouping function, defined in Equation 13: g r o u p c = 0 if   c + 1 ≤ 10 1 if   c + 1 ≥ 11 . (13)

GNN is then trained to solve a 19-class node classification problem on the basis of this input data and outputs the classification results. The training conditions of GNN are summarized in Table 2. Using this approach, we can construct a GNN capable of accurately estimating defect locations from DSPSS.

Algorithm 1 outlines the training pipeline for a GNN model that predicts the three-dimensional location of defects from normalized DSPSS data. It covers data preprocessing, model architecture, distributed training with focal loss, test-time inference, and evaluation. In this study, GAT was constructed using three GATConv layers with hidden dimensions [64, 256, 256] and four attention heads. The attention mechanism is applied at each layer to effectively aggregate the input features. These aggregated features are then passed to a fully connected layer that performs final classification into 19 classes.

Algorithm 1

GNN Training.

Input:Normalized Coordinates { X } ,

To enhance training efficiency, distributed data parallelism was employed, allowing parallel computations across multiple GPUs. To assess the model's generalization capability, stratified five fold cross-validation was carried out. The full dataset was randomly divided into five equal-sized folds. One fold was used as the validation set, whereas the remaining four were used for training, so that every sample was evaluated exactly once. This random splitting procedure guarantees that the model's performance is assessed on diverse, nonoverlapping portions of the data, providing a reliable estimate of its capability to generalize.

During the training process, model evaluation was conducted at each epoch, and early stopping was applied if no performance improvement was observed, thus preventing overfitting. The best-performing model in each fold was saved and evaluated using the test dataset. Evaluation metrics included precision, recall, F1-score, and the confusion matrix, all of which were visualized to interpret performance. Finally, the model that achieved the lowest validation loss among all folds was selected as the final model for performance evaluation.

In this study, we address an imbalanced classification problem in which the number of intact nodes significantly exceeds that of nodes containing defects. To handle this imbalance, focal loss (Lin et al., 2017), rather than conventional cross-entropy loss, is employed. Focal loss increases the loss contribution from misclassified examples whereas it decreases the loss contribution from well-classified ones, thereby encouraging the model to focus more on difficult-to-classify samples: in this case, the nodes contain defects. Since intact nodes dominate the dataset, their contribution to the overall loss is down-weighted accordingly.

The estimated probability p t is defined using p , as shown in Equation 14. p t = p if   y = 1 1 − p otherwise. (14)

Using this definition, we express the focal loss using Equation 15: FL p t = − α t 1 − p t γ ⁡ log p t , (15) where γ is the focusing parameter that controls the degree of down-weighting for well-classified examples and α t is the weighting factor.

In this study, two evaluation metrics are used to quantitatively assess the accuracy of defect prediction: the planar defect location accuracy R and the defect insertion layer prediction accuracy P gauss . These metrics independently measure how accurately the model predicts the planar position and depth of each defect. The total defect prediction score (TDPS), defined as the average of these two metrics, serves as a unified metric for evaluating model performance on each defect case.

In addition, to evaluate the classification performance across the entire test dataset, in this study, we also adopt the macro-averaged F1-score. This metric is used to calculate the F1-score for each class individually and then takes the arithmetic mean across all classes, enabling fair evaluation even when the class distribution is imbalanced.

In the following subsections, we describe the definition of each metric in detail.

As shown in Figure 5a, presents the DSPSS in the physical coordinate system ( x , y , z ) , accurately representing the curvature of the actual structural surface.

Figure 5b projects the same data into a dedicated visualization space. By introducing a virtual coordinate system ( x ′ , y ′ , z ′ ) ,the originally curved surface can be effectively unwrapped and flattened facilitating the inspection of spatial patterns.

Figure 6 presents the surface stress distribution obtained by FEM analysis along with vertical and horizontal stress profiles. The left panel shows the normalized DSPSS when a defect is inserted into the second layer ( 45 ° ) . A significant stress drop is clearly visible around the defect region.

The central graph displays the vertical stress profile along the line passing through the center of the defect. A steep stress decrease is observed in the region corresponding to the defect. The right panel shows the horizontal stress profile across the defect center, indicating stress reduction effects caused by both the defect and the nearby hole.

Figure 7 shows the stress profiles in cases where a single defect is inserted into each of the 2nd–10th layers against the intact scenario. To visually distinguish the effects across layers, layers with the same ply angle are plotted using the same color: 45 ° layers (2nd, 9th) in blue, − 45 ° layers (3rd, 7th) in green, 90 ° layers (4th, 8th) in red, and 0 ° layers (5th, 6th, 10th) in purple. The plots reveal that the impact range and stress reduction patterns vary depending on the ply angle. For instance, in the 90 ° layers, a more abrupt and deeper stress drop is observed than in the other layers, suggesting that the relationship between the fiber orientation and the tensile direction significantly affects stress propagation. Conversely, layers oriented at 0 ° and − 45 ° tend to show more gradual stress gradients.

In the vertical stress profile shown in Figure 7b, each defect-inserted layer exhibits a distinct stress reduction around the defect center, hat clearly deviates from the intact stress distribution. This implies that the presence of defects can be quantitatively identified from surface stress information alone.

The vertical direction in Figure 7b corresponds to the z ′ -axis, which aligns with the fiber direction of 0 ° in Figure 7a. Detailed characteristics according to the ply angle include the following.

• 2nd and 9th layers ( 45 ° ) : Owing to the angled fiber orientation, stress disperses more broadly, resulting in a wider area of stress reduction.

Similarly, in the horizontal stress profiles in Figure 7c, a significant reduction in stress is observed near the defect center, corresponding to the area affected by the defect. Compared with the healthy profile, all layers exhibit consistent stress drops, indicating that the model accurately captures the horizontal spatial positions of defects.

Figure 8 shows the results of predicting the three-dimensional location of defects using the DSPSS obtained from FEM simulations as input. Figures 8a–d correspond to cases where defects were inserted into the 11th layer ( 45 ° ) , 10th layer ( 0 ° ) , 8th layer ( 90 ° ) , and 2nd layer ( 45 ° ) , respectively. Each case includes (i) the input data, (ii) ground truth labels, and (iii) outputs predicted by the model.

Figure 8a represents the casewith the highest TDPS (0.92), which is the average of the planar defect location accuracy R and the defect insertion layer prediction accuracy P gauss . Figure 8b shows the case with second-highest TDPS (0.89), whereas Figures 8c,d correspond to the case with the second-lowest (0.50) and lowest (0.48) TDPS respectively.

Even in cases with low TDPS values, Both the planar position and the insertion layer are generally predicted with reasonable accuracy, as visually confirmed. This suggests that the proposed model successfully learns geometric features of internal defects in multilayered CFRP structures from DSPSS.

For all test data, the minimum value of the planar prediction metric R was 0.55, indicating that the model can generally localize defects in the plane accurately. In many cases, the predicted defect region is slightly overestimated compared to the ground truth. On the other hand, the minimum layer prediction score P gauss was 0.38, indicating that the model is more sensitive to misclassification into neighboring layers or false detections in depth.

Figure 10 shows the distribution of TDPS for each region (left, center, right), corresponding to the planar position where the defect was inserted. Here, TDPS is defined as the average of the planar defect prediction accuracy R and the defect insertion layer prediction accuracy P gauss . By plotting the TDPS for different regions, the figure illustrates how prediction performance varies with spatial location. The horizontal axis indicates each defect insertion layer along with its corresponding ply angle.

As shown in Figure 10a, the data points classified into the center region tend to exhibit higher TDPS than those classified into the left and right regions.

This tendency is primarily attributed to the characteristic of defects in the center region, which vary only in the z ′ –axis direction and are located at a single position along the x ′ –axis in the visualization space ( x ′ , y ′ , z ′ ) . In other words, since the center exhibits the narrowest spatial variation and the most consistent defect pattern, the model can more effectively learn representative features.

On the other hand, defects in the right and left regions are distributed widely in both the x ′ – and z ′ –directions, resulting in greater diversity and making the learning task relatively more difficult. Nevertheless, the performance difference in TDPS among these regions is modest, suggesting that the model can handle asymmetric fields with reasonable accuracy. In summary, prediction accuracy tends to increase in the order from the narrower search space: center, right, and left regions which are indicated by the yellow, red and blue respectively, in Figure 9.

Regarding the planar position prediction metric R shown in Figure 10b, the center region again demonstrates higher accuracy, whereas greater variation is observed in the left and right regions.

Regarding the defect insertion layer prediction metric P gauss shown in Figure 10c, the variation across depth is generally smaller, but the center region still achieves higher accuracy, indicating that the uniqueness of defects in the training data contributes to improved model performance.

Additionally, by examining the layer-wise trend, we observe that the TDPS are highest when defects are inserted into deeper layers (9th–12th), indicating more accurate defect recognition by the model. In contrast, when defects are inserted into the outermost layers (2nd–4th and 17th–19th), TDPS tends to decrease slightly.

Figure 11 illustrates the comparison among the input data, ground truth, and model predictions for various defect locations within the same insertion layer (10th layer). The five subfigures correspond to distinct planar regions as defined in Figure 9. Despite differences in stress field distributions due to proximity to the holes and edges, the model successfully localizes the defect regions with high accuracy. Notably, the prediction performance remains consistent in both symmetric and asymmetric stress regions, demonstrating the robustness and generalization capability of the proposed GNN-based approach.

In all these cases, the TDPS remains high, with the highest being 0.89 and the lowest 0.82. This indicates that even when defects are located in regions with concentrated, peripheral, or nonuniform stress distributions, the proposed method maintains high prediction performance. These results suggest that the proposed approach is not sensitive to particular stress patterns and is robust across various stress distributions. Therefore, the model can stably detect internal defects by appropriately capturing subtle variations in stress fields caused by defects.

Figure 12 shows the confusion matrix for visualizing the prediction results for all nodes across 278 test data instances. The vertical axis represents the ground truth, whereas the horizontal axis represents the predicted classes. Each cell indicates the number of nodes classified into each category. For each ground truth class, precision and recall were calculated to derive the F1-score. The macro-averaged F1-score was obtained by averaging the F1-scores across all classes, resulting in 61%.

In the green-line-enclosed area, the number of nodes misclassified as intact despite actually having defects is extremely low. This indicates that the proposed method rarely fails to detect defects and achieves high detection accuracy. The yellow-line-enclosed area indicates the number of nodes predicted as having defects when there were actually no defects, which correspond to false positives observed around defect edges in Figure 8.

The red-line-enclosed diagonal region represents the correctly classified nodes, and a large number of correct predictions can be observed. However, frequent misclassification into neighboring layers is also noticeable, which contributes to the reduction in macro-averaged F1-score. Misclassifications are most frequent in the deepest layers, specifically the 10th and 11th layers.