Defining the Thermal Features of Sub-Surface Reinforcing Fibres in Non-Polluting Thermo–Acoustic Insulating Panels: A Numerical–Thermographic–Segmentation Approach

Abstract

Natural fibres present ozone-friendly solutions in the field of construction. The attenuation of the sound and heat losses is an important feature in such type of materials above all, when used in non-woven fabrics and fibre-reinforced composites. Hemp fibres show robust insulation performance; this research work should be considered beneficial to the development of a non-destructive thermographic methodology, which can address the thermal barrier (typically applied on multi-layer panel) effects. The intent is to assess the integrity of the sub-surface reinforcing glass fibres; such integrity state will help confer the rigidity and the resistance to mechanical stresses. The testing proposed in this study can be further developed in a laboratory right after the manufacturing process of similar type of components. The testing needs preliminary numerical simulations to help guide the selection of the appropriate pre- and post-processing algorithms combined with or without segmentation operators. A set of numerical and experimental tests were performed through controlled thermal stimulation while recording the thermal responses. The study also highlights the advantages, disadvantages, and future development of the presented technique and methodologies.

1. Introduction

Non-destructive testing (NDT) of composites, including active infrared thermography (IRT), remains an important issue in civil engineering [1]. IRT needs a thermal stimulus to generate thermal imprints [1]. In many construction structures, the composite core is a reinforcing grid made up of basalt [2], carbon [3], or glass fibres [4]. Infrared NDT has proven to be efficient in inspecting glass fibres, but there is still a lack of knowledge in using IRT techniques to assess the composite integrity, when combining glass fibres with hemp fibres scattered in a matrix, to be used as protective coating.

In recent years, a variety of IRT techniques have attracted the attention of many researchers thanks to its capabilities in fatigue assessment, defect detection, and damage repair applications for fibre components or products. The application of IRT for inspecting wind turbine blades, made of GFRP and sandwich panels, was investigated in [5], where surface defects, such as air bubbles, edge bonding, etc., were efficiently detected in blades. In [6], the use of lock-in IRT method for detecting small temperature changes caused by energy dissipation and crack propagation during the tensile test was studied. The evaluation of the IRT effectiveness in identifying different types of defects in masonry structures was performed by Angiuli et al. [7]. Thermographic surveys were also repeated after repair works to identify an acceptance criteria that can be used on areas subjected to total or partial repair. Pitarresi et al. suggested combining techniques of a modified transverse crack tension test coupon and IRT in order to evaluate the static and fatigue behaviour of GFRP composites [8]. In 2020, Xu et al. explored the feasibility of eddy current pulsed thermography to detect sub-surface defects in GFRP composites, a typical non-conductive material. The method introduced a metal part as the heat-generating source to produce the induction heat [9]. The work of Rani and Mulaveesala highlighted advantages of pulse compression in modulated thermal wave imaging, which can be carried out with a moderate peak-power heat source during a certain time period with improved test resolution and sensitivity [10].

Beyond the use of a single IRT method, the combined application of IRT with other NDT techniques, such as acoustic emission (AE), ultrasound imaging, and three-dimensional (3D) computed tomography, demonstrated advantages in numerous scenarios and became an effective solution for NDT practitioners. Karger-Kocsis et al. studied two types of glass fibre composites by applying the AE and IRT [11,12]. Czigány et al. applied AE and IRT on glass fibre-reinforced poly (ethylene terephthalate) composite sheets, concluding that both techniques are suitable for damage detection and assessment of the fibre structure [12]. Ibarra-Castanedo et al. investigated the possibility of detecting defects and assessing the impact severity of delaminations in glass fibre-reinforced polymer (GFRP) by using pulsed IRT and vibro-thermography [13]. Katunin et al. tested various techniques, including piezoelectric ultrasound, IRT, and vibration-based NDT methods, to analyse their applicability for the inspection of aircraft components [14]. Palumbo et al. characterized GFRP adhesive joints by applying ultrasonic imaging and lock-in IRT analysis for the assessment of adhesion quality in composites subjected to static tensile mechanical tests and accelerated aging cycles [15]. Crupi et al. investigated the failure mode and the integrated damage of the sandwich panels by using three-dimensional (3D) computed tomography and IRT. The combined experimental and theoretical studies demonstrated that the use of glass–epoxy reinforcement in aluminium honeycomb sandwiches enhances both energy absorption and load carrying capacities [16]. In 2021, Quattrocchi et al. compared air-coupled ultrasonic testing and active IRT by performing experimental tests on two reference composite panels (GFRP and CFRP) with synthetic defects having different sizes and located at various depths. The pros and cons of the two techniques were highlighted in terms of defect size sensitivity [17].

The continuous development of detection technologies is important to ensure the accommodation of new applications with high precision requirements. Obviously, it is difficult to achieve the desired results via a single detection technique. At the same time, with the development of machine learning algorithms, more image processing methods are used as a post-processing of IRT scans. In order to enhance and expand IRT detection performance and capabilities, Ghorbel et al. investigated the evolution of temperature distributions in glass fibre-reinforced polyamides. In addition to IRT, digital image correlation (DIC) was used to quantify deformed zones and correlate them to heat dissipation sources [18]. Omer and Sultan proposed a new method for post-processing of IRT data by converting temperature read-outs into a pulse train and applying the discrete Fourier transform. Promising results were obtained even for small defects located at the depth of 10 mm below the surface [19]. To increase the detection capability of IRT, a novel IR image processing technique involving the complex Morlet wavelet phase analysis was developed in [20]. It was demonstrated, using the concept of signal-to-noise ratio, that the proposed technique was able to separate the low- and high-frequency noise components from the thermal image signal, thus reliably identifying defects at different depths. Yanjie et al. used the phase Fourier analysis to determine defect location and applied the logarithmic second derivative time to calculate defect depth by using the technique of long pulse thermography [21]. The experimental results were compared with the numerical simulations of a GFRP panel containing fabricated defects [21]. Vijaya Lakshmi et al. proposed a classification- and a regression-based quantitative post-processing technique to characterize sub-surface defects in carbon fibre-reinforced polymers (CFRP) and GFRP specimens by using quadratic frequency modulated thermal wave imaging [22]. In 2019, Sutthaweekul et al. determined the location of flat bottom holes in coated GFRP pipes by using synthetic aperture radar tomography in combination with principal component analysis (PCA), while the defect depth was characterized by analysing the time-of-flight characteristics of the background-subtracted responses [23].

Most of the aforementioned studies are based on using IRT technologies as the primary approach. In fact, numerical simulations and mathematical models are often used in conjunction with experiments to guide and calibrate the test process or assist in developing the pre- and post-processing algorithms. For example, Pracht and Swiderski described the results of numerical simulations in the IRT inspection of multi-layer composite structures, including GFRP [24]. Bernegger et al. presented an approach to quantify delaminations in semi-transparent GFRP. The experimental data was reconstructed by fitting experimental and theoretical profiles using a one-dimensional mathematical model. The delamination depths were estimated with a reasonable accuracy [25]. Hu et al. investigated the characterization and imaging of localized thickness loss (LTL) in GFRP using microwave NDT. A two-dimensional (2D) finite element model with the Ka-band open-ended waveguide and GFRP sample subject to LTL were set up and adopted for this analysis [26].

Herein, we describe specific numerical–experimental tests conducted on a multi-layer specimen. The results were integrated into a segmentation procedure to overcome existing shortcomings in traditional methodologies. Considering the current limited scientific knowledge about non-polluting thermo–acoustic materials [27,28,29] combined with GFRP (that falls into the field of “science and technology for the built environment”), the novelty of presented research will be described in terms of the numerical and experimental techniques and accompanying results. The simulation in Comsol Multiphysics^® software was performed following a purposeful approach aimed at the minimization of computational cost, and the thermal tomography technique developed at Tomsk Polytechnic University (Russia) allowed the detection of detachments at different depths. Surprisingly, the use of local tangent space alignment (LTSA) technique, when applied to the thermographic data and then combined with the Canny segmentation operator, allowed a complete shape (morphology) retrieval of a grid pattern (made by glass fibres), which was located beneath an insulator material in hemp fibres.

The specimen has therefore been analysed by using two numerical models and two different thermal stresses. This helps the reader to understand pros and cons of each of them as the added value of the work. The useful considerations regarding the masking features of noise in respect to the thermal visualization of sub-surface features are presented in the text; also, a comparison between innovative post-processing algorithms applied to the analysis of thermal images (to our knowledge, this is a novel technique applied to composites in the construction sector) is also presented.

The remainder of the article is organized as follows: In Section 2, the description of the multi-layer specimen and main information about numerical–thermographic tests are supplied. Section 3 and Section 5 describe the modalities of numerical simulations performed by Comsol Multiphysics^® and ThermoCalc-3D software; Section 4 and Section 6 describe the different sets of produced thermal images. Section 7 concludes the study by highlighting the outcome of testing and simulation and discussing the results.

2. The Multi-Layer Specimen, the Numerical–Thermographic Tests

The numerical model was developed to study heat conduction peculiarities in a specimen made of materials used in civil construction. Specifically, the specimen was mainly fabricated with a load-bearing base of an expanded Styrofoam (see point 5 both in

Table 1. Dimensional and thermophysical properties of the materials under analysis.

n.	Material	Thickness [m]	Conductivity [W/mK]	Density [kg/m3]	Specific Heat [J/kgK]	Emissivity ε	Porosity %Vol.
1	Scattered Hemp Fibres	0.001	0.038	25	1700	0.9	0.120
2	Cement Mortar	0.0005	1.73	900	0.21	0.54	0.22
3	Fiberglass	0.0002	0.035	21	1030	0.75	0.44
4	Cement Milk	0.0001	1.4	1540	0.87	0.92	0.147
5	Expanded Styrofoam	0.0286	0.03	30	1450	0.6	0.497

and in Figure 1a) and a cement milk layer (see point 4 both in Table 1 and in Figure 1a). A reinforcing grid of fiberglass was embedded in the cement milk layer (see point 3 both in Table 1 and in Figure 1a). This aimed to homogenize possible mechanical stresses that might appear in the base structure and propagate toward the more superficial finishing layers. Moreover, thermal stresses resulting from having different materials with different expansion effects can appear in multi-layer structures as the in-plane deformations of varying magnitude. The studying of such material is important due to its favourable mechanical (i.e., very flexible, etc.) and chemical (i.e., not subject to oxidative phenomena, etc.) nature, yet no established testing regimes are available in the literature.

Above the grid and the cement milk, a layer of cement mortar (see point 2 both in Table 1 and in Figure 1a) was applied while leaving an open space (see Figure 1). In this space, the grid is visually detectable (see point 3 both in Table 1 and in Figure 1a). Such sample design allowed evaluation of the thermal insulation effects conditioned by the cement mortar. Above the cement mortar, the layer of hemp fibres (see point 1 both in Table 1 and in Figure 1a) was present in a part of the inspected sample. The length of this layer was only 8 cm, but it covered the entire width of the specimen. This cover layer highlights the thermal insulation effect of the hemp fibres. For complete details, a numerical description of the stratigraphy in terms of size and thermophysical properties is presented in Table 1.

Figure 1 shows a couple of visual images (a,b), and the numerical models analysed in the first test (c,d).

The specimen was placed on a low-conductive support and illuminated by two projectors equipped with Osram Siccatherm^® 250 (W) lamps and located at −45° and +45°, while the axis of each lamp was oriented towards the centroid of the specimen.

The thermal imaging camera was placed on the bisector of the total angle at a greater distance than that between the specimen and the projectors. In this setup, the camera viewing direction did not intercept the light beams from the projectors during the heating period. The testing scheme was simulated in Comsol Multiphysics^®. In particular, two modelling cases were investigated, of which choice was dictated by the high computational cost required for analysing the model, which will be described below.

In the first modelling, the measurement environment was reproduced by placing the projectors and the specimen at appropriate distances by preserving the real projector-specimen angles.

Figure 2a shows the first test scheme.

In the second test (Figure 2b), the specimen was investigated by applying the classical one-sided thermal non-destructive testing (TNDT) procedure. The sample was heated with a 2.5 kW halogen lamp for 5 s. The temperature distributions were recorded using an Optris PI450 IR imager with an acquisition frequency of 15 Hz.

3. Modelling via Comsol Multiphysics^® Software

The first test was modelled using Catia^® V5 processor reproducing both the heating sources and the specimen without reconstructing the specimen stratigraphy. The external control volume—open boundary—was also defined, and the corresponding governing equations will be explained below.

The modelling defined the electrical power ensured by the lamp filaments. In particular, having selected tungsten as the filament material, the electrical load curve was modelled in Comsol environment to provide Joule heating of the lamps in a time span of 60 s. In order to analyse thermal transients during the cooling phase, the software was developed to evaluate the thermal phenomena when the power supply ceased for 180 s, i.e., simulating the cooling stage. Noting that natural convection was also taken into account.

The model did also accommodate the evaluation of heat distribution and associated fluxes on both the projectors and the sample by taking into account radiation view factors.

All the information obtained from the specimen surface temperature mapping (seen as a temperature footprint) was transferred in the form of input data into the second numerical model. In the framework of the second model, it was possible to analyse the behaviour of the sample in stratigraphic terms without considering the physical environment. Such procedure seemed to be efficient and most appropriate because a hypothetic general model combining the optical features of heating and the entire sample stratigraphy would have resulted in a high computational cost in terms of the accepted numerical mesh (by requiring a RAM of at least 64 Gb [30]).

Figure 3 shows both the image of the mesh (Figure 3a) and the results from the numerical modelling (Figure 3b,c) inherent to the thermal distribution induced by the projectors on the specimen at certain computation times.

Figure 3a shows the tetrahedral mesh, which affects the heating sources with a notable thickening in terms of the number of nodes. Building a fairly dense mesh was necessary to analyse the behaviour of the heating effects in terms of the volumetric analysis of heat flux distribution. Similar mesh specifications were applied to the specimen. In particular, it included a triangular mesh, thus allowing the simulation on the sample sides without calculating sample internal structure. A volumetric tetrahedral mesh was set for the outer cube representing the control volume. For this mesh, a thinning factor of 1.5 was accepted, thus allowing the dense mesh both in the proximity of the agitation source and on the surfaces of the specimen and becoming sparse far from the sample. Such approach has considerably reduced the computational cost without a detailed simulation of the control volume external surfaces.

From Figure 3b, it is possible to evaluate the presence of the control volume represented by the external cube that had the unique function to limit the volume in the analysis of thermal fields. In order to simplify the display, it was decided to use the slice function with the section plane passing through the resistances of the heating sources.

The aforementioned modelling strategy allowed us to study the temperature trends of the tungsten filament. In addition, the isocontour function was used to evaluate the increase in air temperature and its distribution. It is evident that the temperature was in a wide range starting from the maximum value (detected at the tungsten filaments of the heating sources) down to the ambient air temperature of the control volume.

Figure 3c shows the thermal field as generated by the heating sources at 10 s. Both the high temperature from the tungsten filament and the effect of the environment are well demonstrated.

The concept to build two separate models was prompted by the complexity of the specimen shape and its unique internal structure. The apparent regularity of the external dimensions (parallelepiped) and the thickness of the layers constituting the model, however, has prevented the use of a sweep mesh. Therefore, the first model, i.e., the specimen with a swept envelope, was discretized by means of a triangular mesh including only 428 nodes, while the second model was used to meet some problems described below.

The main constraint that complicated the use of a regular mesh along the specimen thickness was the presence of the reinforcing grid. The latter, for a finite areal extension, has no layer of material on it, thus challenging implementation of a swept mesh for all layers of the model. In fact, there was no correspondence between the upper and lower surfaces in terms of the number of nodes.

A tetrahedral approach was used for producing the mesh, which would fit the reinforcing grid. Given its geometrical dimensions, both a nodal thickening equal to 0.0001 m and a resolution of narrow regions equal to 0.25 were necessary. These values were obtained through a series of iterations aimed to evaluate the most efficient thickening with reasonable computational efforts and associated processing cost.

Discretization of the cement milk layer also imposed some difficulties. The addition of the warp and weft of the fiberglass inside the cement milk resulted in the formation of a series of square-based parallelepipeds. The mesh of the cement milk layer required the same degree of accuracy as the fiberglass. Given the regularity of the bottom of the layer (a single parallelepiped), it was possible to reduce the number of nodes through controlling the maximum element grow rate parameter in the direction of the depth. Concerning the rest of the layers, the mesh was built following the automatic generation process. The model simulating the whole specimen had a number of nodal elements equal to 372 346. Because of the mesh complexity involving a high degree of freedom (DOF), it was difficult to ensure a single processing step, which would involve the physics of radiation, including both the projectors and the effects of global convection. In fact, the quantities calculated in the first model were added to the second model, as both the convection and the mutual irradiation of the specimen towards the outside ambient were validated. Furthermore, in the second model, the control volume was no longer present. The latter was not necessary in the first modelling effort, as only the mutual radiation towards infinity was involved. The simplification of the specimen in terms of the first model mesh and the thickening of the specimen DOF in the complete second model prevented the use of the multi-component technique typically provided by Comsol Multiphysics^®. If the specimen remained unchanged, the sets of equations could have been solved separately within the same model at different times by avoiding the physical passage of the variables.

The first model consisted of the heat transfer physics with surface-to-surface radiation set of equations as follows:

$$\rho C_p\frac{\partial T}{\partial t}+\underset{1}{\underbrace{\rho C_p\mathit{u}·\nabla T}}=\nabla ·\left(k\nabla T\right)+Q+Q_r,$$

(1)

where ρ indicates the density in (kg/m³), C_p is the specific heat at constant pressure (J/kg·K), T the represents temperature (K), u represents the velocity field tensor (m/s) (in our case, term 1 of Equation (1) is null because the components are mutually immobile), Q indicates the generic heat source/sink (J), and Q_r is the radiative heat source/sink (J). On the other hand, Equation (2) represents the basis of the radiative phenomena conditioned by the projectors:

$$\underset{1}{\underbrace{S_i·\nabla I_i}}=\underset{2}{\underbrace{k{I}_b\left(T\right)}}-\underset{3}{\underbrace{\beta I_i}}+\underset{4}{\underbrace{\frac{{\sigma }_s}{4\pi }{\displaystyle \sum }_{j=1}^N{\omega }_j{I}_j\mathrm{\varnothing }\left(S_i,S_j\right)}}.$$

(2)

The term 4 of Equation (2) shows the presence of a discretization for the angular space; the term 4 is usually presented in the form of an extended integral from 0 → 4π to analyse all spatial directions. In this case, as there were only two projectors providing the thermal load, the term 4 degenerates into numerical quadratures of discrete directions. Within the summation, σ_s indicates the scattering coefficient divided by 4π, as the function is spatial. The remaining part, reported as a member of the summation, indicates the phase function. The latter evaluated the probability that a generic ray from the direction S_i is projected in the direction S_j. The definition of the phase function depends on the material constituting a single layer, while ω_j indicate the j-directions.

The first model does not include a detailed view of the specimen (i.e., internal structure), yet its main purpose is to compute the radiation effects of the projectors and environment appearing on the surfaces of the specimen. To fully understand this concept, it is necessary to refer to the corresponding problem by following Rosseland [31].

In Equation (2), the summation has N as its upper term, but in the considered case, it is equal to 2. In particular, the term 1 of Equation (2) represents the gradient of the radiative intensity I as the i-th component projected along the S direction. Note that, in the term 2 of Equation (2), k represents the absorbed fraction of the radiative intensity evaluated for the black body I_b(T). In the term 3 of Equation (2), the fraction β indicates the radiant intensity I_i with respect to the generic i-th direction.

Of particular interest are the boundary conditions of the control volume in the first model, which represent open boundaries being defined by Equation (3):

$$\begin{array}{c}T=T_0\hspace{1em}if\hspace{1em}\mathit{n}·\mathit{u}=<0\\ -\nabla T·\mathit{n}=0\hspace{1em}if\hspace{1em}\mathit{n}·\mathit{u}=\ge 0\end{array}$$

(3)

where n is the normal outgoing to the control volume.

Both the lighting bodies and the surfaces of the simplified specimen have boundary conditions of the type:

$$-\mathit{n}·\left(-k\nabla T\right)=\mathit{\epsilon }\left(G-e_b\left(T\right)\right),$$

(4)

where G indicates the incoming radiative heat flux or irradiation (W/m²), e_b(T) is the blackbody hemispherical total emissive power (W/(m²·T³)), and $\mathit{\epsilon }$ provides the surface emissivity, i.e., a dimensionless number in the range 0 ≤ $\mathit{\epsilon }$ ≤ 1.

In particular, as in the first model, the specimen is simplified in terms of geometry. The aim is to set the maximum value of G for each specimen surface; therefore, $\mathit{\epsilon }$ is equal to 1, i.e., being independent of the material. In Equation (4), the k parameter represents the thermal conductivity. Instead, the second model, where the sample is reproduced in all details involves the physics named Heat Transfer with Radiation in Participating Media. In the framework of this physics, the values of G produced by the calculation on the surfaces of the simplified specimen (see the first model) were considered as the input values for the entire model yet to be analysed. However, in the latter model, both the real emissivities of each material and the actual specimen shape (e.g., the geometry of the reinforcement grid) were taken into account to calculate the real thermal load. Then, the governing equation of G on the surfaces becomes:

$$q=G-\mathit{\epsilon } e_b\left(T\right).$$

(5)

Equation (5) allows the calculation of the thermal energy on each surface exposed to the thermal load through considering the actual emissivity and other radiative effects of the media involved. In addition, the set of equations intended for the analysis of the porous materials calculated by means of the physics for porous average in scattered fibres, cement mortar, cement milk, and Styrofoam is given below:

$${\left(\rho C_p\right)}_{eff}\frac{\partial T}{\partial t}+\rho C_p\mathit{u}·\nabla T=\nabla ·\left(k_{eff}\nabla T\right)+Q+Q_{vd}+Q_p,$$

(6)

where,

$$\underset{1}{\underbrace{{\left(\rho C_p\right)}_{eff}}}=\underset{2}{\underbrace{{\mathit{\theta }}_{\mathit{p}}{\rho }_p{C}_{p,p}}}+\underset{3}{\underbrace{\left(1-{\mathit{\theta }}_{\mathit{p}}\right)}} \rho C_p$$

(7)

and k is

$$k_{eff}={\mathit{\theta }}_{\mathit{p}}{k}_p+\left(1-{\mathit{\theta }}_{\mathit{p}}\right) k,$$

(8)

where Q_vd in Equation (6) indicates the energy of the viscous dissipations in the form of heat (J), while Q_p is the point heat source (J). In Equation (7), the term 1 indicates the effective volumetric heat capacity at constant pressure, while the term 2 shows θ_p, i.e., the dimensionless value indicating the fraction of porosity present in the material, ρ_p is the density of the porous material (kg/m³), and C_p,p represents the specific heat capacity (J/(kg·K)) for the porous material. In Equation (7), the term 3 indicates the non-porous percentage of the material expressed as a dimensionless number. Finally, in Equation (8), k_p indicates the thermal conductivity of the porous material.

Figure 4 shows the mesh linked to the second numerical model. This figure is a visual display demonstrating the meshing strategy deployed for the sample.

Numerical Results

The analysis was carried out on the second model, where both the heating and cooling phases can be studied (Figure 5 shows the specimen heating phase).

Figure 5a shows the temperature range calculated in the framework of the second model by using the first model input data. The entire volume of the specimen was subjected to thermal analysis. The slicing technique was not used to take advantage of the volume perspective. It is seen that after only 10 (s) of heating, the sample front surface was thermally loaded in a prevalent way. The calculated thermograms illustrate that the lowest temperature increase appeared on the bottom of the specimen. The insulating effect of hemp fibres was evident and prevailing in relation to the remaining layers. The layers of cement mortar (both in lower and upper positions) were more affected by the temperature increase than the cement milk layer. It is evident that the reinforcement grid, which was installed along the entire specimen, reached a higher local temperature in regard to the cement milk layer.

In Figure 5b, it is possible to notice a uniformity of the temperature field on the front surface, with a reduced thermal gradient among all points up to 1.1 (K). On the one hand, the grid below the cement milk layer became less evident due to the homogenization of the temperature differences between the surface layer and the remaining subsequent layers. On the other hand, it is possible to note the effect of temperature gradient increase along the specimen thickness due to the thermal load, which is mainly caused by the front surface. In Figure 5c, the sub-surface grid below the cement mortar layers (lower and upper position), is less evident. This depends on homogenization of the temperature related to the surface layers.

The diffusion of the thermal front can be also observed along the specimen depth. Figure 5d displays a further homogenization in the temperature map specific to the glass fibre grid; such trend was not evident below the cement milk layer, where the grid contours were clear and well defined. Throughout the remaining part of the specimen, the sub-surface reinforcement grid became less observable.

In Figure 5e, the diffusion process involved both layers of the cement mortar; therefore, it was difficult to define the grid geometry underneath them. On the other hand, the mesh under the cement milk layer was still visible. This can also be attributed to the thickness of this layer, which provided a lining to the glass fibre.

Figure 5f shows the temperature distribution at the end of the heating period. However, the maximum specimen temperature was observed later because of the thermal inertia of the projectors, which continued to illuminate the sample surface. This phenomenon is illustrated by Figure 6a showing the temperature distribution at 10 s after the projectors were switched off.

The temperature distribution in Figure 6b is fairly uniform except the cement milk layer. The hemp fibre area is characterized by a small temperature signal of −0.37 K that can be explained by natural convection. The similar temperature pattern with the temperature gradient on the surface of −0.07 K is shown in Figure 6c; also, the in-depth penetration of the temperature front is clearly seen. Respectively, in Figure 6d, the footprint of the reinforcement grid beneath the cement milk layer becomes less pronounced, and this phenomenon evolves in time (Figure 6e) to reach the temperature gradient between the glass fibres and cement milk of about 0.04 K at 120 s or 180 s after heating started (Figure 6f). In both Figure 6g,h the specimen displays a uniform temperature distribution for all the layers of interest. The reinforcement grid is barely visible below the cement milk layer, and the respective temperature difference is below 0.03 K. To further analyse the in-depth temperature, four different directions were identified. The corresponding segments shown in Figure 7 passed through the intersections of the grid weft and warp in order to calculate influence of glass fibres. The corresponding temperature profiles are presented in Figure 8.

As follows from Figure 8a, there are certain differences between the profiles calculated for the different materials at 10 s. Of particular interest is the comparison between the two layers of cement mortar where the temperature profiles are conditioned by the material thickness. Lower surface temperatures are observed on the surface of cement milk and hemp layers. Respectively, some temperature variations take place at the location of the reinforcing grid.

In general, the penetration depth of the temperature front can be evaluated to be around 8 mm; see the temperature profiles in Figure 8b–f. For example, the temperature increase at the depth of the reinforcing grid is negligible while comparing the two layers of cement mortar. However, between the cement milk layer and the hemp fibres, there is a temperature difference of 0.75 K. A noticeable increase of surface temperature (0.1 K) appears over hemp fibres to compare to the cement milk area (Figure 8c). The abovementioned features of heat conduction in the specimen remain the same for all times (Figure 8d,e).

In Figure 8e, the trends of the temperature field shows a very slight increase. On the one hand, the behaviour of the directrix that affected the hemp fibres can be noticed. On the other hand, no particular differences concerning the abovementioned trends are visible. This did not happen between the cement milk and the hemp fibre layers, as the latter increased its temperature until it exceeded the surface temperature of the former. The same effect was also evident in regard to the temperature of the reinforcement grid.

The profiles in Figure 8f correspond to the end of the heating period, demonstrating a general temperature increase on the sample surface, in particular, in the area over the hemp fibres where the temperature exceeded that of the cement mortar by about 0.5 K.

Figure 9 illustrates changes in the temperature profiles at the cooling phase along the same directions as shown in Figure 7. As mentioned above, for a certain period of cooling, the surface temperature was slightly increasing due to thermal inertia of the projectors (Figure 9a–d) with a slightly higher temperature over the hemp fibres. This can be explained by the higher thermal insulation properties of the hemp fibres. In general, all profiles were similar to those observed within the heating stage, but the temperature penetration depth apparently increased up to about 19 mm at 120 s. In their turn, the profiles in Figure 9e–f show that the residual heating from the projectors ceased at about 180 s; therefore, one may expect the appearance of the temperature inversion phenomenon conditioned by the layers’ different heat capacities. Note also that, at ~180 s, the temperature front reached the depth of 29 mm, which is about the specimen thickness.

4. Processing Experimental Thermograms: The First Session

Experimental IR image sequences (spatial format of 126 × 240), which were provided by active thermal NDT, represented 3D “data cubes” to be further transformed into 2D matrices (spatial format of 30,240 × 240); see Figure 10.

In Figure 10, X is the unfolded 2D matrix. At the next processing step, the data matrix was normalized by subtracting the average of each column (column centralization), and the normalized matrix was still denoted as X. In X, the rows and columns represented samples (thermal images) and variables (pixels), respectively. Then, three manifold learning methods, namely isometric feature mapping (ISOMAP), local linear embedding (LLE), and LTSA, were used to reduce the dimension of X, respectively. It should be noted that it was not necessary to analyse all thermal imagery by applying the aforementioned codes; only the first 60 thermal images were taken as the input for further data processing.

The matrices X were reduced to d-dimension by using the ISOMAP, LLE, and LTSA algorithms. In order to visualize the dimension reduction results, the partial least squares regression (PLSR) algorithm was chosen to approximate the mapping process from X to Y.

4.1. Experimental Results

4.1.1. ISOMAP

While using ISOMAP (Figure 11), the following parameters were accepted: d = 6 and neighbours = 4. After ISOMAP data processing, all six indicator vector (IV) images contained varying degrees of the inhomogeneous background, as well as noise, which obscured information on reinforcing grid. Some grid features can be observed in the IV2 image, but it is difficult to clearly detect the meridians and latitudes of the grid because of the noise. The IV3–IV5 images contained both the severe inhomogeneous background and noise, thus making demarcation at particular specimen sections difficult.

4.1.2. LLE

The parameters of the LLE algorithm were (Figure 12): d = 6 and neighbours = 6. Observing the six IV images of the LLE results, it can be said that the inhomogeneous background (caused by the inhomogeneity of the specimens) is attenuated in the images when compared to the ISOMAP results. In particular, the grid feature is more pronounced in the images of IV4–IV6 as compared to ISOMAP. This may be attributed to the fact that LLE is more sensitive to inhomogeneous information, as it retains the local intrinsic structure of the data unlike the ISOMAP method, which learns the global manifold of the data.

4.1.3. LTSA

The LTSA parameters were set as follows (Figure 13): d = 6 and neighbours = 9. Similarly to the previous two approaches, the 6 IVs of LTSA still had difficulty in identifying the grid lines. Compared with the ISOMAP and LLE results, the LTSA results had more severe inhomogeneous background and associated noise, while the grid features were clearer. In particular, the upper part of the IV2 image contained the weak inhomogeneous background, and the grid lines were more visible in regard to other images.

In summary, none of the results can clearly display the grid shown in area 1 (Figure 1a) because of noise. Therefore, it was necessary to use edge detection algorithms to further analyse the thermographic results.

Here, the authors select indication vector (IV) 2 in LTSA as a case, and use the Canny, Sobel and Laplacian operators to extract signal gradients (Figure 14).

Applying the edge detection algorithm resulted in a better display of the grid outline in the area 1 (Figure 1a). However, the Laplacian operator could only extract a small amount of grid information in the area 2 (Figure 1a), and the Sobel operator also showed a poor denoising ability. In contrast, the Canny operator was able not only to extract the grid information in the areas 1 and 2, but also to effectively remove the noise (Figure 14).

When implementing the Canny operator, some pre-treatments should be conducted to achieve better results. One method is to save the IV2 image (Figure 13) with a resolution of 300 dpi, and another method is to enlarge the IV2 image (Figure 13) directly where each pixel in the original image corresponds to a region of 8×8 in the enlarged image. The respective results are shown in Figure 15.

It is visually seen that the second method produced better results. It was compared with the ground truth image manually obtained by means of the Photoshop^® software through the replicating area 3 above the areas 2 and 1 (Figure 1a). The ground truth image is shown on the left side of Figure 16. However, it is still difficult to give a quantitative comparison between the Canny and ground truth results, since the numbers and positions of grids did not perfectly fit. This means that the ground truth should be obtained in a different way in a future work, perhaps by removing the layers above the grid in the areas 1 and 2 in order to acquire a photograph covering the entire area.

5. Modelling via ThermoCalc-3D^® Software

As shown in Figure 17a, the plane Section 1-1 (Figure 1b) was modelled in order to analyse the possibility to detect the fiberglass grid located beneath scattered hemp fibres and embedded into the layers of cement at the depth of ~1.5 mm from the sample surface. At such depth, air-filled defects can easily be detected due to their high thermal contrast in regard to the host material. The thermal properties of fiberglass were comparable to those of hemp fibres, and the detection of the grid seemed to be possible due its thermal contrast in respect to cement mortar. In the 3D model of Figure 17a (material thermal properties are given in the caption of the figure), the specimen was heated for 5 s using a uniform heat flux, and the purpose of the modelling was to evaluate the dimensionless temperature contrast C = DT/T_nd, where DT = T_d − T_nd, with the subscripts “d” and “nd” specifying temperatures in defect and non-defect areas, as shown in Figure 17a [32].

The images of the contrast C are displayed in Figure 17b, and its time evolution, in Figure 17c. The contrast reaches the maximum within the heating period; then, it decays and finally experiences one more maxima at the end of heating that is a typical behaviour of the contrast in active TNDT (note that the differential signal DT is characterized by a single maxima). At longer times, the contrast becomes negative because of the higher heat capacity of fiberglass with respect to that of the cement mortar, which is the abovementioned signal inversion, but its amplitude is low; the reader may compare the first and last images in Figure 17b.

The maximum contrast is about C = 0.012, or 1.2%. This is a fairly low value to ensure reliable detection of the grid; for example, if the sample is heated up to 10 °C above the ambient, the temperature signal DT produced by the grid will be only 0.12 °C. In the ideal case of no noise, the grid pattern can be underlined, for example, by the thermal tomogram shown in Figure 17d; it is worth noting that the grid pattern in the tomogram is distorted because of heat diffusion. In general, any material can be characterized by the surface noise contrast $C_n={\sigma }_n/{\displaystyle \stackrel{=}{T_{nd}}}$, which includes both the additive and multiplicative noise components; here, ${\sigma }_n$ is the standard deviation of noise within a chosen reference area, and $\stackrel{=}{T_{nd}}$ is the mean temperature in the reference area [33]. Noise contrast varies in time and depends on an area chosen for calculating ${\sigma }_n$. For hemp fibers, the value of $C_n$ is presented in Figure 17c with the dotted line. The noise is higher within the heating stage because of reflected radiation with a fixed value of about 1.2% within the cooling stage. By comparing the contrasts produced by the grid and the surface clutter (Figure 17c), one may conclude that the fiberglass grid cannot be detected on the noisy background (in some cases, a regular structure of a defect pattern may be helpful for identifying the corresponding structures).

6. Processing Experimental Thermograms: The Second Session

By applying the Principal Component Analysis (PCA) to the experimental results, one can clearly observe the non-uniform structure of the mortar layers (Figure 18), which appears due to varying thickness of the mortar, as well as some inter-layer detachments occurring on the boundaries of the fibre glass grid.

The analysis of the area 1 (Figure 1a and Figure 18c) shows a vague pattern of the grid structure underneath the hemp fibres, which act as a thermal insulator. As mentioned above, in one-sided TNDT, some depth-resolved results can be obtained by using the technique of dynamic thermal tomography [34]. The basic concept of this technique involves the analysis of time delays of differential temperature signals between a defect and sound areas. Each class of structural defects can be described by a particular calibration curve allowing the evaluation of a defect depth by the corresponding time delays. Material diffusivity and thickness also affect calibration results. Some thermal tomograms of the tested specimen are shown in Figure 19 to reveal no grid pattern, as predicted by the noise analysis. However, thermal tomography illustrates variations of mortar thickness and some peculiarities of inter-layer boundaries.

7. Discussion and Conclusions

New building enclosure materials, its assemblies, and systems used for minimizing and/or regulating space heating and cooling modes are today at the focus of several research projects and gaining more interest from many sectors.

In this work, a specific multi-layer composite (useful to limit heat losses from buildings) was analysed numerically and experimentally for integrity. In particular, Comsol Multiphysics^® and ThermoCalc-3D^® software and a couple of active IRT experiments were implemented to detect a grid in glass fibres positioned under an external layer made of hemp fibres. The numerical part revealed some matches, such as the evaluation of the thermal/signal inversion and exposed challenges in detecting the grid on the noisy background. In general, the ThermoCalc-3D software^® can be considered more user-friendly in regard to the modelling performed on the Comsol platform, even if the latter approach was optimized to reduce computational cost.

In order to detect the reinforcing grid, some new data processing approaches were proposed and demonstrated. Namely, the IV2 and LTSA techniques were complemented by the use of the Canny, Sobel, and Laplacian operators. The best results were provided by the Canny algorithm for the presented case.

The one-sided TNDT (usually used in building and material inspections) [35,36] was also applied in this study, specifically in the second experimental session. Data processing was performed using the techniques of PCA and thermal tomography. No grid pattern, as predicted by the noise analysis, was visible by thermal tomography, although the varying thickness of the mortar and some peculiarities of inter-layer boundaries were revealed. The second component of PCA slightly revealed the presence of the sub-surface grid.

No quantitative analysis was possible at this stage of research as the extraction of an accurate ground truth image inherent to the grid pattern still requires an objective strategy to allow quantitative assessment.

In conclusion, we have demonstrated that: (a) the use of advanced software, such as Comsol Multiphysics, is useful to precisely simulate a multi-layered specimen within an accepted heat conduction model; (b) a thin external layer of hemp fibres ensures a significant shift of thermal phase in construction materials; (c) the use of ThermoCalc-3D^® software allowed us to link the maximum thermal contrast with the noisy background and understand why a sub-surface fiberglass grid cannot be detected by using only optical heating; (d) among some advanced data processing approaches used, such as ISOMAP, LLE, and LTSA, the LTSA algorithm demonstrated better results; (e) the thermal tomography technique was able to reveal detachments between specimen layers better than other processing techniques; and (f) the LTSA algorithm provided the best input image for the subsequent segmentation procedure by means of the Canny operator.