Second-Order Parabolic Equation to Model, Analyze, and Forecast Thermal-Stress Distribution in Aircraft Plate Attack Wing–Fuselage

: During a flight, the steel plate attack wing–fuselage of an aircraft is subjected to cyclical thermal stress caused by flight altitude variation that could compromise the functionality of the plate. Thus, it is compulsory after a sequence of flights to evaluate the state of plate health. In this work, we propose a new dynamic model on the basis of the physical transmission of heat by conduction governed by a second-order parabolic partial differential equation with suitable initial and boundary conditions to analyze and forecast thermal stresses in the plate of a P64 OSCAR B airplane. Developing this model in the COMSOL Multi-Physics R (cid:13) environment, a finite-element technique was applied to achieve the thermal-stress map on the plate. The achieved results, equivalent to those obtained by a campaign of infrared thermographic experiment measurements (not yet used in the aeronautical industry), highlight the evolution of the thermal load of the steel plate attack wing–fuselage, adding evidence of possible incoming fatigue phenomena to identify in advance if the steel plate must be replaced.


Problem Introduction
Design in the aeronautical industry is currently increasingly based on the use of accurate calculation codes and on the precise knowledge of the physical characteristics of the materials constituting the elements of an aircraft [1][2][3][4]. Moreover, to reduce both management and production costs, it is imperative to require extreme performance of the materials through in-depth knowledge of their mechanical resistance [1,5,6], especially when fatigue phenomena occur due to cyclic loads such as lift on wings, gusts, take-off/landing maneuvers, and fuselage pressurization [7]. Even if cyclic loads are smaller than breaking ones, they contribute to the formation of microcracks, which inevitably lead to the collapse of the structure [8]. In the history of civil aviation, air disasters caused by collapse due to material fatigue have been numerous [9]. For this reason, the aeronautical industry pays attention to studying the problem of fatigue due to cyclic loads [7,8]. In an aircraft, plate attack wing-fuselages, characterized by reduced geometries, are the structural elements most exposed to the risk of failure because they are subjected to high stresses [8]. In fact, these plates, from take-off to landing, are loaded from bottom to top by the aerodynamic lift and from top to bottom by mass forces due to the weight and loads distributed along the wingspan (tanks, engines, and carriages) [10]. Moreover, highly dangerous thermal stresses due to high temperature gradients arising from the flight altitude variations during the flight could occur. In fact, there is a sudden drop in temperature during take-off until the cruising altitude is reached, while, during the descent maneuver, temperature increases drastically until landing, generating a load (which becomes cyclical due to the high number of flights operated by an aircraft), whose effects must be carefully monitored [10]. If cyclic mechanical loads (wind, maneuvers, fuselage pressurization) do not arouse particular concerns thanks to a highly sophisticated and certified design, analysis, and maintenance procedures, thermal stresses are still a source of concern because they are evolving phenomena in delayed time regimes compared to heat sources. In other words, the effects (mechanical stresses) are hysteretic concerning the causes (high temperature gradients). Because the simultaneous use of specific analysis techniques [11] is required to limit management/production costs, it is preferable to exploit Nondestructive Testing and Evaluation (NDT and E) for analyzing structural elements ensuring their functionalities and, if necessary, to guarantee their reuse [12]. Usually, extraordinary maintenance interventions on an aircraft require NDT and E to evaluate each element through certified quality-control procedures after a certain number of flight hours to provide indications, with sufficient safety, on the number of load cycles beyond which the element must be replaced. Technologies offer a wide possibility of NDT procedures that assess the state of health of the structural elements in an aircraft. NDT techniques based on penetrating liquids are traditionally used in the aeronautical industry because they are simple to perform and offer fast responses with a good degree of accuracy in defining the "superficial" state of health of the element under investigation, while, on a punctual level, and especially in the deeper layers, they hardly provide comprehensive information [13]. In fact, this technique exploits the "wetting" properties of particular liquids and is used to detect surface defects. Furthermore, the impossibility to carry out remote checks with penetrating liquids, together with the almost impossibility to automatize the procedure for real-time applications, obliges using other techniques. Ultrasonic Techniques (UT) do not lend themselves to this type of flanking because the Doppler effect on the echo signal invalidates the measurements [12]. Furthermore, the presence of roughness in many metallic structural elements severely limits the use of Eddy Currents (ECs) because, besides presenting the skin-effect problem [14,15], they do not detect defects arranged parallel to the magnetic field [16,17]. Then, to satisfy the need, on the one hand, to also assess the integrity of a structural element locally and, on the other hand, to assess dynamically the evolution of thermal stresses, it is necessary to use, together with penetrating liquids (or alternatively), NDT techniques able to analyze in detail the local variations of temperature as a function of time. In this regard, Infrared thermography (IR), not yet used in the aeronautical industry, appears to be an excellent candidate because it quickly and accurately returns in real time surface temperature maps of the investigated elements (thermograms). Furthermore, by analyzing the surveyed areas in real time and repeating the measurements over time, it is possible to check changes in surface temperature on several points simultaneously, even when the environmental conditions change (video storage) [18]. Moreover, the possibility of data storage and processing by specific software makes IR particularly attractive for our purposes. In this work, a new approach to evaluate the conditions of structural elements of civil aircraft subjected to cyclic thermal stress, on the basis of implementation in the COMSOL Multiphysics R environment, is proposed. In particular, attention was paid to the implementation and analysis of the steel plate attack wing-fuselage of a P64 OSCAR B subject only to cyclical thermal stresses caused by different temperatures due to flight altitude variation in both static and dynamic conditions. The proposed model, based on the physical transmission of heat in solids by conduction, led to a dynamic model governed by a parabolic Partial Differential Equation (PDE) [19] that, by means of appropriate initial and boundary conditions, was made by Finite Element Method (FEM) techniques [20,21]. The obtained results, superimposable on those obtained by a campaign of IR experiment measurements under the same conditions, assess the state of punctual solicitation of the plate and, specifically, its time evolution, highlighting possible incoming fatigue phenomena, a priori identifying the plate intended to be replaced. The paper is structured as follows. After a brief overview of the technical characteristics of a P64 OSCAR B aircraft (Section 2), the proposed model is detailed in Section 3. Specifically, Sections 3.1 and 3.2 discuss in detail the mathematical physics background and initial and boundary conditions, respectively. Then, germane and useful lemmas are presented in Section 4. The description of the realization of the model in the COMSOL Multiphysics R environment is described (Section 5) detailing the physics and geometry (Section 5.1), the determination of the heat-transfer coefficients, and both mesh creations and the choice of significant temperatures, respectively. Relevant results are presented in Section 6 depending on whether they are in steady-state or dynamical conditions. The experiment investigations by thermal imaging and interesting comparisons with the proposed model are reported in Section 7. Specifically, after having shown through an investigation with penetrating liquids that the plate has no detectable defect (Section 7.1), a brief overview of IR thermography is described in Section 7.2. In Sections 7.3-7.5, we show the thermographic investigation at different temperatures. Finally, some conclusions and future perspectives end this work.

Technical Characteristics of the P64 OSCAR B Aircraft
At the Laboratory of Micromechanics and Materials for Aeronautical and Aerospace Technology of the Department of Civil, Energy, Environmental, and Material Engineering (DICEAM), Mediterranea University of Reggio Calabria (Italy) is located a P64 OSCAR aircraft under renovation (Figure 1), whose technical characteristics are detailed in Table 1. Each wing of this aircraft is made up of two half-wings, each of which is connected to the fuselage by two pin-type connections and an upright connection. Each half-wing consists of a central body connected at the front to polyester resin and a glass fiber shell reproducing the front part of the wing profile. The used alloy is aluminum-copper 2024 characterized by good workability characteristics for machine tools, by high resistance and toughness, good resistance to atmospheric corrosion, and high resistance to high temperatures. The central body is made up of a front (main) side element and a "false" element (secondary) that make up the vertical walls. A series of ribs and cladding panels complete the structure. Flaps and ailerons are both made by a structure consisting of a single spar, ribs, and cladding panels, also in alloy 2024. Each half-wing houses a fuel tank; at the end of the main spar and of the false element are the plate attack wing-fuselages (see Figure 2a,c).

Plate-Attack Wing-Fuselage: Geometrical and Physical Characteristics
The plate is 14 cm long, has a thickness of 0.5 (cm), and weighs 146.3 (g). It also has four holes whose diameter is 0.7 (cm) to accommodate the bolts for fixing to the main longitudinal element and a hole of 1.5 (cm) in diameter to insert the pin for fixing to the fuselage (see Figure 2b,d). The used material was steel 3630N containing about 0.3% C, about 1% Cr, and about 0.2% Mo, which give the material high mechanical properties suitable for aerospace applications. In it, there are also ions Mn and Si. Impurities due to S and P that are present with a lower percentage do not compromise the mechanical properties (the thermal and mechanical parameters of steel 3630N are listed in Table 2); unlike other alloys, it does not undergo alterations due to aging.

Heat Equation
The heat equation (or linear diffusion equation) for a function T(v, t) (in our case, the absolute temperature), with v ∈ R n the spatial variable and t > 0 the temporal variable, is a PDE taking the well known form [22][23][24]: where D is a positive constant (diffusion coefficient), ∆ is the Laplacian operator, and f (v, t) is the external heat source. Equation (1) describes the evolution of T(v, t) in a homogeneous and isotropic medium, with constant density ρ, receiving heat from a source. To achieve Equation (1), we indicate by r the rate of heat per unit of mass supplied to the body from the outside, and we consider a volume V. The energy-balance law requires that the rate of change in internal energy in V equal the flow of heat through ∂V (frontier of V), due to conduction, plus the contribution due to the external source. If e is the internal energy per unit of mass (internal energy density), the total amount of internal energy in V is given by: Let us denote by q the heat-flux vector. Then, given an infinitesimal surface dσ of the edge of V centered in v ∈ V with external unit vector n, q · ndσ expresses the amount of energy flowing through dσ. The scalar product is consistent with the observation that there is no flow of energy if vector q is normal to n. The incoming heat flow through ∂V is given by: where we used the Gauss theorem. Finally, the contribution due to the external heat source is given by: The energy balance, therefore, requires the following relationship:

Constitutive Laws
Let us now assume two constitutive relations. The first is Fourier's law for heat conduction, according to which q is proportional to T(v, t) as follows: where thermal conductivity k is a constant positive related to the properties of the material [25]. Obviously, in some applications, k cannot be considered constant, but depends on T. The minus sign in Equation (6) is due to the fact that the heat flows from high temperature towards low temperature zones. The second constitutive equation links internal energy e and absolute temperature T(v, t): where c p is the specific heat at constant pressure of the material. Combining Equations (5)- (7) and taking into account that the Reynolds theorem is applicable, we obtain [19]: from which, since Equation (8) is true for each V, then: We thus obtain Equation (1) in which D = k/c p ρ and f = r/c p .

Initial and Boundary Conditions
To complete the formulation of the problem under study, we need to specify both initial and boundary conditions for Equation (1). In particular, regarding the initial conditions, we impose T(v, 0) ∀v ∈ V in the plate at t = 0. That is: Regarding the boundary conditions, it is possible to set the Dirichlet condition, specifying ∀v ∈ ∂V: in which T q is the external temperature. Alternatively, the Neumann boundary condition specifies flux of heat ∀v ∈ ∂V as follows: Finally, Newton's law of cooling, can be exploited, in which h represents the heat-transfer coefficient. Condition (13) treats the heat flux through the surface, which depends on the temperature difference between the surface and environment. If radiant heat transfer between the body and the surrounding medium (or vice versa) becomes considerable, the following nonlinear condition: can be exploited. In this relation, σ is the well known Stefan-Boltzmann constant, A is the surface, and is surface emissivity, respectively [19].

Remark 1 (Parabolic frontier).
No conditions were assigned at t = t * , with t * the final time. The condition was only assigned on the so-called parabolic frontier ∂ P Q T of cylinder Q T given by the union of base V × {t = 0} and lateral part S T = ∂V × (0, t * ]:

Steady-State Case
If T is independent of t (that is, T(v)), Equation (1) becomes Poisson's equation: otherwise, it takes the form of the Laplace equation: Remark 2 (General formulation of the steady-state case). Equation (16) can be considered as a particular case of the following problem (in Dirichlet form): where D, σ 1 , f , and b are assigned functions (or constant). Generally, In these cases, the solution can give limit layers, that is regions next to ∂V where the solution is characterized by strong gradients of T. Obviously, if both b · ∇T(v) (convective term) and σ 1 T(v) (reactive term) are negligible, Problem (18) becomes Equation (16).

Preliminary Lemmas
As described above, both the Laplace Equation (17) and heat diffusion Equation (1) need suitable boundary conditions and initial-boundary conditions, respectively. However, for applicative reasons, it is imperative to ensure that T(v) and T(v, t) distributions, in both cases (steady-state and dynamical conditions), are uniquely determined to know the physical-mathematical conditions ensuring their existence and uniqueness. The following lemmas provide exhaustive indications about the existence and uniqueness of the solutions in both cases.
Lemma 1 (Laplace equation: maximum principle). If T ∈ C(V) has the mean property in V (i.e., the value at the center of each sphere B R ⊂⊂ V equals the average of the values at the boundary ∂B R ) and v * ∈ V is a global maximum or minimum point for T, then T is constant. In particular, if V is bounded and T ∈ C(V) is not constant, then ∀v ∈ V: Lemma 2 (Laplace equation: uniqueness, comparison, and stability for the Dirichlet problem). Let us suppose that g ∈ C(∂V). Then, the problem has, at most, a solution belonging to C 2 (V) ∪ C(V). In addition, if g 1 ≥ g 2 on ∂V and g 1 = g 2 in at least one point on ∂V, then T g 1 > T g 2 in V. Finally: Heat always flows towards regions where the temperature is lower. This has as the consequence that the solution of Equation (1) takes its global maximum and minimum values on ∂ p Q T (maximum/minimum principle). Moreover, T(v, t) is independent of any changes in the data after t. Mathematically, we have the following lemma, which applies to continuous functions up to the edge of Q T with the temporal first derivative and spatial second derivatives continuous inside Q T . Here, these functions are indicated by C 2,1 (Q T ) ∩ C(Q T ) [19].
Then, the maximum (resp. minimum) of T is on ∂ P Q T : In particular, if T(v, t) is negative (positive) on ∂ P Q T , then it is negative (positive) in Q T .
As a consequence of Lemma 3, we have that if: the maximum and the minimum of T are located on ∂ P Q T . In particular, Lemma 3 is a version of the weak maximum principle because it does not exclude the possibility that the maximum (minimum) can also be located on a point not belonging to ∂ P Q T . Lemma 4 (Heat equation: comparison, uniqueness, and stability for the Cauchy-Dirichlet problem). Let us consider T 1 (v, t) and T 2 (v, t) as solutions in C 2,1 (Q T ) ∩ C(Q T ) of: and: respectively, with f 1 (v, t) and f 2 (v, t) bounded in Q T . Then, 2. the following estimation of stability holds: In particular, the Cauchy-Dirichlet problem admits a unique solution.

Model Physics
By the heat-transfer module of COMSOL-Multiphysics R , combined with the structural-mechanics module, we analyzed the thermal conduction on the plate attack wing-fuselage of a P64 OSCAR B aircraft. This is because these interfaces combine mechanical problems with problems due to thermal exchanges. This combination occurs when T(v, t) acts as a thermal load for the mechanical interface, producing thermal expansion. As a first step, a steady-state condition was selected. After that, we studied the dynamical step.

Model Geometry
To perform the simulations, we first designed the steel plate in a 2D environment space of the COMSOL CAD interface, reproducing the geometry as shown in Figure 3a. Then, we implemented an extrusion, obtaining the final 3D model as depicted in Figure 3b.

Parameter Setting
This step required setting the mechanical and thermal parameters that are characteristic of the material constituting the steel plate as listed in Table 2.

Fixed-Constraint Setting
For our purposes, the fixed constraints of the plate on the main longitudinal component were selected through five holes, as shown in Figure 3b.

Plate Heat Transmission
As is known, heat transmission occurs by conduction, convection, and radiation [26]. However, heat transmission does not exclusively occur in one of the three modes. Heat conduction is governed by Equation (1); Equation (13) governs convection mode; and finally, the Maxwell-Hertz-Plank equation describes the radiation mode [23,26]. However, for our purposes, we just considered heat transmission by conduction because the heat through convection and radiation, in our case, can be considered negligible [26].

Mesh Creation
First, the partition of the 3D domain was made by extending the selected faces of the geometry by extending a face on the adjacent face of a domain to partition it. Partitioning was used for the creation of tetrahedral meshes, subdividing the complex form domain (highly nonconvex) into a simpler mesh domain [23]. Furthermore, this could increase the level of parallelization (shared memory) during mesh generation. If the original geometry was meshed, unstructured tetrahedrons had to be used, which may not be the optimal choice. In order to achieve a more robust tetrahedral mesh and to avoid errors concerning failure to respect boundary-edge elements on the geometry face or internal errors in boundary respecting, we used a new algorithm involved in the "free tetrahedral" function in COMSOL MultiPhysics R . In particular, the tessellation method was automatically set. First, the standard method was applied. If it failed, the new algorithm was used. However, in both cases, the Delaunay procedure was utilized [23]. In particular, the generation of the mesh was unstructured, so the triangulation was carried out in such a way that the sphere circumscribed to each tetrahedron did not contain vertices within it. This kind of triangulation was unique and maximized the minimum angle of the grid tetrahedra. However, to make the Delaunay procedure applicable to non-convex domains (such as of the plate at hand), we used the constrained Delaunay triangulation, which allowed fixing a set of sides of the mesh to be generated: the resulting grid necessarily associated these sides with some tetrahedron. In particular, we imposed the sides that defined the mesh boundary. It is worth pointing out that the constrained Delaunay triangulation is not a Delaunay triangulation because some of its tetrahedra may contain vertices belonging to the initial set. However, the vertices were only and exclusively the original ones specified in the set, and no others were added [23]. Finally, we obtained a mesh consisting of 3623 volume elements, 2208 surface elements, and 407 edge elements (see Figure 4).

Choice of Significant Temperatures
Once the geometry of the plate was implemented and the problem was set from a physical-mathematical point of view, it was necessary to choose significant external temperatures to take into account the most important phases of a flight. This link derived from the fact that T q (external temperature), starting from sea level, almost linearly decreases by 3.25 ( • C) for each 500 (m) of altitude, up to 11.000 (m) (36,089 (ft)), beyond which and up to 22,000 (m) (82,021 (ft)), it remains constant at the value of −56.5 ( • C). We therefore chose to simulate the behavior of the plate at sea level and 2500 and 5000 (m) corresponding to 0, 15, and 35 ( • C), respectively, taking into account that a P64 OSCAR B in safe conditions can reach a maximum altitude of 5000 (m) (for details, see Table 3). Table 3. Significant temperatures depending on flight altitude.

Boundary-Condition Setting
In this work, we used temperature T(v) = T q on ∂V (Dirichlet's condition (11)) as a boundary condition. However, COMSOL-Multiphysics R requires, by default, the computation of heat-transfer coefficient values h related to significant temperatures depending on flight altitudes as specified in the following subsection.

Determination of Heat-Transfer Coefficients h
To determine h, we combined Equations (12) and (13), obtaining [26,27]: in which q represents the exchanged heat flow (Watt), A is the area of the heat-exchange surface (m 2 ), and T q − T s is the difference of temperature ( • K). Moreover, combining Equations (12) and (14), we obtain: so that Equation (28), considering Equation (29), becomes: where = 0.07 is the emissivity of polished steel, σ is the Stefan-Boltzmann constant (W/m 2 • K 4 ), and T q and T s are the room and the plate temperatures, respectively. Setting = 0.07, T s = 311.15 ( • K) (polished steel), σ = 5.67 × 10 −8 ( W/m 2 • K), we obtained the required values of h related to significant temperatures depending on flight altitudes (see Table 4).

Simulations in Steady-State Conditions
As above specified, heat transfer in the plate in this work was assumed to be by conduction so that Equation (1) could be selected in COMSOL MultiPhysics R in its steady-state version (∂T(v)/∂t = 0).
Moreover, since the external heat source was null, in which T q represents the temperature to which the plate is subjected during each simulation. As is known, Problem (31) can be considered as a particular case of the general Problem (18) (in Dirichlet form) in which the reactive term is neglected and b · ∇T(v) = 0 because b represents a velocity that, in steady-state conditions, vanishes. Then, Problem (18) becomes: If D is considered as a constant, then Problems (31) and (32) coincide. Then, a set of simulations in steady-state conditions was performed. Figure 5a-c depicts the 3D simulations subjecting the plate to different T q values. In particular, in Figure 5a, T q = 0 ( • C) at sea level, T q = −16.25 ( • C) if the flight altitude is 2500 (m), and T q = −32.50 ( • C) when the flight altitude is 5000 (m) (see Table 3). Analogously, in Figure 5b, T q = 15 ( • C) at sea level, T q = −1.25 ( • C) when the flight altitude is 2500 (m), and T q = −17.50 ( • C) when the flight altitude is 5000 (m). Finally, in Figure 5c, T q = 30 ( • C) at sea level, T q = 13.75 ( • C) when the flight altitude is 2500 (m), and T q = −2.50 ( • C) when the flight altitude is 5000 (m). The simulations of the thermal maps were obtained with an uneven color over the entire surface, but they are clearly distinguishable from the more and less intense areas. Red highlights the critical zones due to both the plate geometry and the type of material used. In contrast, blue shows plate areas that are poorly affected by thermal stress.

Remark 3 (On the existence and uniqueness of thermal maps in steady-state conditions).
It is worth noting that the fundamental solution of Problem (31) is [19]: where n is the dimension of the space and α(n) is the volume of the ball of unit radius, so that, for n = 3, it makes sense to write: Then, ∀v ∈ V, the solution T(v) exists. In fact, in the steady-state simulations, v = 0, so that T(v) is physically bounded. In addition, fixing a values for T q , Lemma 2 guarantees the uniqueness for T(v), ensuring that the thermal maps in each simulation are unique when T q is fixed.
Remark 4 (On functional space for T(v)). By this set of simulations, it can be seen that the most critical areas concerned the portion of ∂V next to the holes where thermal stresses were strongly present. This was congruent with the fact that solution T provided limit layers next to ∂V characterized by a strong gradient of T (see Remark 2). It was also congruent with Lemmas 3 and 5. Experimentally, T(v) ∈ C 2,1 (Q T ) ∩ C(Q T ) because both the thermal and mechanical properties of the steel used constituting the plate admit continuous-time variation of the temperature and very reduced but continuous changes of its second spatial derivatives in Q T .

Remark 5 (On the location of the minimum and maximum values of thermal stress on the plate).
Consider that ∂T/∂t = 0, f (v) = 0, and ∆T = 0, as highlighted in Equation (21). Then, maximum T values took place on ∂V (zones around the holes). On the other hand, far from the holes, thermal stresses were very low. In particular, in some portions of ∂V far from the holes, the thermal stress assumed minimum values (see Equation (21)). This was again congruent with Lemmas 3 and 5 because the maximum and minimum values were located on ∂V. To confirm this experimentally, we chose two generic points of the plate far from the holes, whose vectors' positions were v 1 and v 2 , respectively. It follows that v 1 − v 2 was sufficiently high (because they were internal to the connected portions of V), so that T(v 1 − v 2 ) resulted in being greatly attenuated. Equation (34) results: then in the area of the plate far from the holes, the temperature was not high. On the other hand, near the holes of the plate, there was an increase of temperature, as can be also noticed by the results reported in Figure 5 (to reproduce as much as possible the edge of the holes, it was necessary to consider small values of |v 1 − v 2 |).
The same considerations described in Remarks 4 and 5 are still valid for the depicted simulations in Figure 5b,c. Although these simulations were in steady-state conditions, they could be considered as a good model giving interesting information about more thermally stressed holes (and consequently, from those mechanically stressed). To give further information about the evolution of the thermal-stress distribution required by both designers and technicians, we needed to develop dynamical simulations when thermal shock was applied. These simulations could be used as a prediction model of T(v) to also provide designers with useful and indispensable support for the optimization of the plate geometry/material or to help technicians with deciding if the plate can be reused after a sequence of flights.

Simulations in Dynamic Conditions
Our problem in dynamic conditions can be written as follows: in which the initial condition is T(v, 0) = g(v) = T q and the boundary condition is T(v, t)| v∈∂V = T q . Equation (36) and related initial and boundary conditions were implemented in COMSOL MultiPhysics R in order to obtain a set of thermal maps in dynamical conditions. Remark 6 (On the existence and uniqueness of thermal maps in dynamical conditions). As is known, the fundamental solution of Problem (36) can be written as follows: which represents a curve very close to a Gaussian one. It can be observed that whatever v in the plate, the solution exists because Γ D is bounded. Furthermore, Lemma 4 guarantees the uniqueness for T(v, t) ensuring that the thermal maps in each dynamical simulation are unique when both initial and boundary conditions are fixed. T(v, t)). Even in dynamic simulations, the same considerations apply as in Remark 4.

Remark 7 (On the functional space for
Remark 8 (On the location of the minimum and maximum values of thermal stress on the plate). As in the previous steady-state simulations and far from the holes of the plates, the thermal stress assumed low values also in dynamical simulations. In fact, considering v 1 and v 2 two different position vectors and for a fixed t = t * , Equation (37) can be written as follows: If v 1 and v 2 lie in the areas far from the holes of the plates, Γ D (v 1 − v 2 , t) assumes low values. In other words, in dynamical conditions and far from the holes, the thermal stresses were low. The contrary happened next to the holes: here, the temperature increased considerably, as well as the thermal stress (to reproduce as much as possible the contour of each hole, it was imperative to consider small values of |v 1 − v 2 |).
To understand how thermal stresses evolve in time t in the plate when thermal shock T q is initially applied, we needed to perform a set of simulations that focused on thermal-stress distributions in the plate for t ∈ [0, +∞) (s). In particular, we focused on the evolution of the thermal-stress distribution for t ∈ {0, 1, 2, 3} (s), with ∆t = 1 s, because it represented a more dangerous time interval for the plate after the application of a thermal shock. In fact, in the plate after 3 s, T started to decrease, so that the more dangerous time interval was [0, 3] s [26].
Remark 9 (On the choice of the dangerous time interval for the plate after the application of a thermal shock). After having fixed a vector position v (in other words, after having chosen a particular point on the surface of the plate), the fundamental solution (37) can be written as Γ D (t). For each t, this was a bell shaped definite positive function such that Γ D (t) → 0 as the time t → 0 and Γ D (t) → 0 as the time t → +∞. Then, setting v and varying the time t, the temperature increased, reaching the stationary point (at the time t * ), and finally decreased, as highlighted in the simulations. For each simulation, t = t * = 3 s, so that [0 3] s represented the dangerous time interval for the plate after the application of the thermal shock.
For dynamic simulations, we also used the physics and geometry of the model, the setting of the parameters and continuum mechanics, and also the mesh creation related to the steady-state simulation as detailed in Sections 5.1-5.6, respectively. Moreover, for analysis uniformity, we chose the same significant temperatures here as those described in Section 5.7. The simulations produced Figures 6a-c, 7a-c and 8a-c. In particular, Figure 6a-c shows the performance in dynamic conditions with initial thermal shock T q = 0 • C (sea level), −16.25 • C (2500 (m)), and −32.5 • C (5000 (m)), respectively. It is easy to argue that after 3 (s) of the application of T q , the plate presented the critical areas strongly highlighted in red. Moreover, in blue, the areas less thermally loaded were evidenced. As for the steady-state simulations, the same consideration discussed in Remarks 4 and 5 were still valid. The obtained results notably did not depend on the use of other parameters such as atmospheric pressure and wing weight, which weighed even more on the areas highlighted in red. The discussion can be easily extended to the other simulations whose performance is depicted in Figures 7a-c and 8a-c. To validate the performance of the results achieved by the implementation of the parabolic model in COMSOL MultiPhysics R , we compared the obtained results with those achieved by an experimental investigation by IR thermal imaging.

Pretreatment with Penetrating Liquids
The plate was pretreated with penetrating liquids to verify the absence of defects (the presence of at least one defect obliges that the plate not be reused). Particularly, its surface was cleaned with a degreasing agent and heated to 35 ( • C), then coated with a fluorescent penetrating liquid. The excedent fluid was removed with a hydrophilic emulsifier. After washing, the penetrating liquid left inside the defects rose to the surface and left a signal with a much larger size than the defect in the detector. The surface was analyzed with a UV radiation lamp that guaranteed high contrast with high detection sensitivity [13]. Analysis did not show any type of defect in the plate, so it made sense to perform a survey for the detection of thermal stresses.

IR Thermography: Brief Overview
This is an NDT measurement of the IR radiation emitted by a body to determine its surface temperature. IR radiation is located between the visible spectrum and radio frequencies/microwaves with a wavelength between 1 (mm) and 0.7 µ(m) (frequencies between 3 × 10 11 and 4 × 10 14 (Hz)). Even if the human eye cannot see it, its presence is perceived by the thermal effect it produces. IR thermography visualizes and quantifies thermal flows through the generation of maps (thermograms) that associate a color with a detected temperature to evaluate the condition of materials. It is necessary to know (or to hypothesize) the distribution of surface temperatures in the absence of defects to compare them with the actual distribution measured during a test. An anomaly in temperature distribution indicates a possible defect, visible by a color spot. In this work, we used a FLIR SC660 camera with a resolution of 640 × 480 and a wavelength between 0.7 and 1 (D m). The achieved images (thermograms) were processed by FLIR RESEARCH IR software. After appropriate calibration of the camera, three types of simulations were conducted: at below-zero temperature (−18 ( • C)), at room temperature (15 ( • C)), and at about 30 ( • C). These temperatures were slightly different from the significant temperatures of −16.25 ( • C), −15 ( • C), and 35 ( • C) as shown in Table 3. This was due to technical laboratory requirements.

Thermographic Investigation at
The plate was placed in the freezer for 15 min to reach −18 ( • C). The achieved thermogram (Figure 9a) shows three critical areas identified by A, B, and C. Using the marker Li1 (shown in Figure 9a in the area identified by A), we obtained the trend of the temperature (see Figure 9b). Figure 9c depicts the trend of temperature identified by the marker Li2 (see Figure 9a) in the area identified by B) in which a sudden change in temperature was detected, but in a wider space with respect to the area A. Finally, Figure 9d shows a clear change in temperature near the point C.
Comparing this result with the performance achieved by COMSOL MultiPhysics R , it is evident that the critical areas identified in the simulation corresponded exactly to those highlighted with the thermographic analysis.

Thermographic Investigation at 15 ( • C)
As in Section 7.3, the procedure was performed at 15 ( • C), obtaining the thermogram shown in Figure 10a. Through analysis, we deduced that the temperature along its upper area (marker Li1 in Figure 10a) showed an erratic trend (Figure 10b) with evident thermal jumps. Analysis along the lower side (Figure 10c) showed the existence of a series of sudden temperature changes confined in limited areas. Even near the holes (Figure 10d), we noted critical areas. Comparing the experiment performance with the simulated performance at 15 ( • C), one could see critical areas, common to the two models, concentrated mainly near the holes.

Thermographic Investigation at 30 ( • C)
With a heat gun, the plate was heated to up to 30 ( • C), obtaining the thermogram shown in Figure 11a. Analysis of the thermogram highlighted the existence of two sudden temperature changes (see Figure 11b), identifying two critical areas. Figure 11c,d highlights the presence of areas that need to be more carefully analyzed. In this case, the overlapping of the critical areas obtained in the simulation phase and in the experiment tests was also evident.

Conclusions and Perspectives
In this work, the authors presented a second-order parabolic model implemented in COMSOL Multiphysics R to analyze thermal stresses (distribution of temperature T) in a plate attack wing-fuselage of a P64 OSCAR B aircraft. Starting from the details of the technical characteristics of the aircraft, the proposed model was mathematically detailed in terms of equations and boundary and initial conditions. In addition, important well known results of the existence/uniqueness and localization of minimum/maximum T values for both steady-state and dynamical conditions, respectively, were presented. Then, the authors detailed the model realization in COMSOL MultiPhysics R , modeling both the physics and geometry of the problem and then setting the parameter and constraints. The 3D mesh was created by an innovative technique exploiting the constrained Delaunay triangulation suitable for the plate under study (non-convex domain).
Once the authors chose significant temperatures of aeronautical interest, they performed a set of simulations in steady-state conditions at different flight altitudes. The simulated results put in evidence critical areas of the plate where the thermal stress appeared significant. This fact was also demonstrated exploiting the fundamental solution for the heat equation in steady-state conditions (Laplace equation), highlighting the minimum and maximum values of the temperature on the plate.
The results also highlighted the model's capability to evaluate the presence and exact position of any critical points (in the steady-state case), also showing the T distribution evolution in time on the plate (in dynamical conditions), putting into evidence and quantifying a dangerous time interval for the plate after the application of a thermal shock. Moreover, in the dynamic conditions, the location of the minimum and maximum values of the thermal stress on the plate were proven by exploiting the fundamental solution of the heat equation. Finally, the validity of the model was confirmed via comparison with the obtained results by IR thermographic analysis. In particular, the results achieved by the proposed model in terms of critical areas of thermal stress on the plate were totally comparable (both in position and intensity) with the results achieved by IR thermographic analysis. The proposed model presented here was a reliable prediction system of critical plate areas, helping technicians to establish if the plate can be reused. Moreover, the dynamic approach fixed the number of flights beyond which the plate could not be reused. The competitiveness of this model lied in the quality and quantity of information that could be extracted from the simulations, making it possible to identify the geometric characteristics of critical areas both in terms of shape and size. However, despite the good agreement between numerical simulations and thermographic images, the proposed model resulted in not being well suitable for real-time applications due to the related heavy computational burden. Anyway, taking advantage of the fact that the transition between low thermal loads and high thermal loads is not crisp, but gradual, it would be worthwhile to model the problem using soft computing techniques. In particular, fuzzy computation techniques could be exploited to create a system handled by a bank of IF-THEN fuzzy rules easily readable by non-experts and, in any case, updatable by the knowledge of an expert. Finally, the reduced computational load of these systems should allow obtaining low cost hardware implementation. This will be the scope of a future work.

Abbreviations
The following abbreviations are used in this manuscript: NDT Heat-transfer coefficient Q T Cylinder ∂ P Q T Parabolic frontier S T Lateral part of ∂ P Q T −∇ · (D∇T(v, t)) Diffusive term b · ∇T(v) Convective term Bounded functions in Q T T s Plate attack wing-fuselage temperature Li1, Li2, Li3 Markers