Dimensionless numerical analysis of cracked materials by means of infrared lock-in thermography

MASTER’S DEGREE IN Space Science and Technology MASTER’S THESIS DIMENSIONLESS NUMERICAL ANALYSIS OF CRACKED MATERIALS BY MEANS OF INFRARED LOCK-IN THERMOGRAPHY Student Sagarduy Marcos, David Supervisor Rodríguez Aseguinolaza, Javier Department Applied physics Academic year 2021-2022 Bilbao, June 24th of 2022 Acknowledgments I would like to thank to the Photothermal Technique Laboratory research group for all the help they have given me, with special mention to my supervisor Javier for the dedication and effort he has put into this work. The many hours we have spent trying to decipher the bash scripts and their entrails are now remembered as the titanic effort that has been made. In these lines I cannot forget my parents. I would like to thank them for supporting me from the very beginning and giving me solutions when I didn’t see them (even when I thought there were none!). Although only my name appears in the authorship of this work, it is as much mine as theirs. Last but not least, I would like to thank my grandparents for the constant encouragement I have received from them over the years. Without it this would have been a much harder road. i Abstract In this master’s thesis a dimensionless numerical model of the lock-in infrared thermography experiment for the detection of open-surface cracks is developed. Starting with the constitutive equation and the boundary conditions associated with the laser and the cracks, their dimensions are removed by introducing length, time and temperature scales related to the physical problem. As a result, a set of dimensionless parameters is obtained, that allows to give a broader vision to the problem, while maintaining the mathematical simplicity of the dimensional model. Once the dimensionless equations are obtained, they are implemented in a finite element method multiphysics software called OpenFoam. With this objective, in this work the full numerical calculation process has been developed: the pre-processing (or meshing), processing and post-processing stages. After implementing the equations in the finite element method software a parametric analysis has been performed by means of different simulations in order to analyze the effect of each of the dimensionless parameter in the resulting amplitude thermogram. Key words: Dimensionless equations, lock-in infrared thermography, finite element methods, OpenFoam, parametric analysis. ii Resumen En este trabajo fin de m´ aster se desarrolla un modelo num´ erico adimensional del experimento de termograf´ ıa infrarroja modulada para la detecci´ on de grietas superficiales. Partiendo de la ecuaci´ on constitutiva y de las condiciones de contorno asociadas al l´ aser y a la grieta, se lleva a cabo el proceso de adimensionalizaci´ on introduciendo escalas de longitud, tiempo y temperatura, asociadas al problema f´ ısico. Como resultado se obtiene un conjunto de par´ ametros adimensionales que permite dar una visi´ on m´ as amplia al problema, manteniendo a su vez la sencillez matem´ atica del planteamiento dimensional. Una vez adimensionalizadas, estas ecuaciones se implementan en un software multif´ ısico de m´ etodos de elementos finitos llamado OpenFoam. Para ello en este trabajo se contemplan las tres etapas de c´ alculo que se han seguido: el pre-procesado (o mallado), procesado y post-procesado. Tras implementar las ecuaciones en el software del m´ etodo de elementos finitos se realiza un an´ alisis param´ etrico mediante diferentes simulaciones con el fin de analizar el efecto de cada uno de los par´ ametros adimensionales en el termograma de amplitud resultante. Palabras clave: Ecuaciones adimensionales, termograf´ ıa infrarroja modulada, m´ etodos de elementos finitos, OpenFoam, an´ alisis param´ etrico. iii Laburpena Master amaierako lan honetan, gainazalera irekitako pitzadurak detektatzeko termografia infragorri modulatuko esperimentuaren zenbakizko eredu adimentsional bat garatzen da. Eratze-ekuaziotik eta laserraren eta pitzadurak sartutako mugalde baldintzetatik abiatuta, adimentsionalizazio prozesua gauzatzen da, arazo fisikoari lotutako luzera, denbora eta tenperatura eskalak sartuz. Ondorioz, parametro adimentsionalen multzo bat lortzen da, arazoari ikuspegi zabalagoa emateko aukera ematen duena, eta, aldi berean, planeamendu dimentsionalaren sinpletasun matematikoari eusten diona. Dimentsioak ezabatu ondoren, ekuazio horiek OpenFoam izeneko elementu finituen metodoen software multifisikoan inplementatzen dira. Horretarako, lan honetan jarraitu diren hiru kalkuluetapak jasotzen dira: aurre-prozesatzea (edo saretzea), prozesatzea eta prozesatu ostekoa. Elementu finituen metodoaren softwarean ekuazioak inplementatu ondoren, analisi parametriko bat egiten da hainbat simulazioren bidez, parametro adimentsional bakoitzak anplitude-termograman duen eragina aztertzeko. Gako hitzak: Ekuazio adimentsionalak, termografia infragorria modulatua, elementu finituen metodoak, OpenFoam, analisi parametrikoa. iv Contents 1 Introduction 1 1.1 Context ......................................... 1 1.2 Objectives........................................ 3 2 Experimental background 4 3 Theoretical model 6 3.1 Dimensionlessequations ................................ 6 3.2 Dimensionlessparameters ............................... 9 4 Computational numerical simulations 12 4.1 OpenFoam........................................ 12 4.2 Numericalschemes................................... 12 4.3 Solution and algorithm control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 4.4 Themesh ........................................ 14 4.4.1 Infinitecracks.................................. 15 4.4.2 Semi-infinitecracks .............................. 16 4.4.3 Finitecracks .................................. 17 4.5 Meshresolutioncontrol................................. 18 4.5.1 Mesh element distribution gradient . . . . . . . . . . . . . . . . . . . . . . 18 4.5.2 Meshrefinement ................................ 21 4.6 Processing........................................ 22 4.7 Post-Processing..................................... 22 5 Results 23 5.1 Validationofthemodel................................. 23 5.2 Parametricanalysis ................................... 25 6 Discussion 31 6.1 Validationofthemodel................................. 31 6.2 Parametricanalysis ................................... 32 6.2.1 Π1........................................ 32 6.2.2 Π2........................................ 32 6.2.3 Π3........................................ 32 6.2.4 ¯y0........................................ 33 6.2.5 Πw,Πl,Πdand θ................................ 33 6.2.6 Non-dimensional parameter combinations . . . . . . . . . . . . . . . . . . . 34 7 Conclusions 37 References 38 List of Figures 1 Typical laboratory lock-in IR thermography experimental setup. . . . . . . . . . . . 4 2 Scheme of the evolution of the surface temperature of the sample with harmonic laser heating. ......................................... 5 3 Scheme of the thermal waves in different points of the sample. . . . . . . . . . . . . 5 4 (left) Hexahedral and (right) tetrahedral meshes. . . . . . . . . . . . . . . . . . . . . 14 5 (left) Tetrahedral and (right) hexahedral mesh elements. . . . . . . . . . . . . . . . . 14 6 Schemeofaninfinitecrack............................... 16 7 (a) The two blocks that make up the infinite crack and (b) the final result. . . . . . . 16 8 (a) The four blocks that make up the semi-infinite crack and (b) the final result after joiningthem. ...................................... 17 9 (a) Lateral view in the +¯xdirection showing the corresponding four blocks and (b) zenital view in the +¯zdirection showing the corresponding three blocks. . . . . . . . 17 10 Final result after joining the corresponding twelve blocks. . . . . . . . . . . . . . . . 18 11 Diagram of the spacing between mesh edges. Left: regular distribution. Right: downwardspatialgrading. .................................. 18 12 (a) Three block mesh without gradient and (b) same mesh with γ= 1 calculated by (40)............................................ 20 13 (a) Single box, (b) triple cylinder refinements and (c) both of them combined. . . . . 21 14 Temperature evaluation of an arbitrary point over cycling until the stationary stability criterionissatisfied.................................... 22 15 Comparison between the non dimensional numerical results and analytical solution for infinite vertical cracks. (a) AISI 304 and (b) Cu. . . . . . . . . . . . . . . . . . . 23 16 Comparison between the non dimensional results and the DG FEM model in the semiinfinite crack case in AISI 304 for two different crack inclinations. (a) 90º and (b) 45º. 24 17 Comparison between the non dimensional results and the DG FEM model in the finite crack case in AISI 304 for two different crack inclinations. (a) 90º and (b) 45º. . . . . 24 18 Natural logarithm plot of the thermal amplitude (a) non-normalized and (b) normalized on the transverse sample profile for Π1= 4 ×103,1×105,2×105,3×105. ............................................. 26 19 Natural logarithm plot of the thermal amplitude (a) non-normalized, showing the zoomed amplitude jump associated to the crack in the insert and (b) normalized on the transverse sample profile for Π2= 1 ×10−7,2×10−7,3×10−7,1×10−8. . . . 26 20 Natural logarithm plot of the thermal amplitude (a) non-normalized and (b) normalized on the transverse sample profile for Π3= 50,500,1000,4000........... 27 21 Natural logarithm plot of the thermal amplitude (a) non-normalized and (b) normalized on the transverse sample profile for ¯y0= 0.24,0.32,0.40,0.48. ......... 27 22 Natural logarithm plot of the thermal amplitude (a) normalized and (b) crack region zoomed on the transverse sample profile for Πw= 0.001,0.002,0.005,0.01. . . . . . 28 23 Natural logarithm plot of the thermal amplitude (a) normalized and (b) crack region zoomed on the transverse sample profile for Πl= 0.2,0.4,0.5,0.75.......... 28 24 Natural logarithm plot of the thermal amplitude (a) normalized and (b) crack region zoomed on the transverse sample profile for Πd= 0.1,0.2,0.3,0.5. ......... 29 25 Natural logarithm plot of the thermal amplitude (a) normalized and (b) crack region zoomed on the transverse sample profile for θ= 50º, 75º, 90º, 100º. ......... 29 26 Natural logarithm plot of the thermal amplitude on the transverse sample profile for two different (Π1,Π3,Πw) triples which satisfy (41). . . . . . . . . . . . . . . . . . 34 27 Natural logarithm plot of the thermal amplitude on the transverse sample profile for two different (Π3,Πw) tuples which satisfy (49). . . . . . . . . . . . . . . . . . . . . 36 2 Experimental background through the domain, characterized by their amplitudes and phases (see figure 3). In this project, the parametric analysis is going to be done with the amplitude. Figure 2: Scheme of the evolution of the surface temperature of the sample with harmonic laser heating. Figure 3: Scheme of the thermal waves in different points of the sample. 5 3 Theoretical model 3 Theoretical model 3.1 Dimensionless equations In order to develop a dimensionless computational numerical model of the crack detection problem by means of lock-in thermography, the first aspect that must be taken into account is the manner in which the dimensions of the set of equations that describe problem are removed. Due to the nature of the experiment, the constitutive equation is the heat equation: ∇2T=1 α ∂T ∂t (4) where αis the thermal diffusivity. Aiming to remove the dimensions of this equation the following variables are introduced: ¯x=x Lx ,¯y=y Ly ,¯z=z Lz ,¯ T=T T0 ,¯ t=t tc (5) being Lx,LyyLzcharacteristic lengths of the problem in the three spatial directions, T0an arbitrary temperature (e.g., room temperature) and tca characteristic time value of the problem. Explicitly writing the heat equation and applying the chain rule: ∇2T=∂2(¯ TT0) ∂¯x2 ∂2¯x ∂x2+∂2(¯ TT0) ∂¯y2 ∂2¯y ∂y2+∂2(¯ TT0) ∂¯z2 ∂2¯z ∂z2=1 α ∂(¯ TT0) ∂¯ t ∂¯ t ∂t (6) Using the previously introduced parameters, this equation reduces to: 1 L2 x ∂2¯ T ∂¯x2+1 L2 y ∂2¯ T ∂¯y2+1 L2 z ∂2¯ T ∂¯z2=1 tcα ∂¯ T ∂¯ t(7) A simplification that can be carried out at this point, and which is particularly convenient for computational calculations, is to choose the same length scale for all three spatial directions. Thus, imposing Lx=Ly=Lz=Lthe dimensionless heat equation reduces to: ¯ ∇2¯ T=∂2¯ T ∂¯x2+∂2¯ T ∂¯y2+∂2¯ T ∂¯z2=L2 tcα ∂¯ T ∂¯ t(8) where ¯ ∇2is the dimensionless Laplacian operator. In the process of removing dimensions of partial differential equations there are multiple criteria to choose the characteristic length and time scales. In particular, for the heat equation there are two main different criteria depending on the characteristics of the problem [14]: 1. tc=L2/α, if heat diffusion occurs significantly throughout the domain. Here Lis the length of the sample. 6 3 Theoretical model 3.1 Dimensionless equations 2. tc= 1/ω, if the change in temperature is significant only up to a certain limited distance, l. Here, ωis an angular frequency associated with the problem. In lock-in IR thermography experiments, heat diffusion is a process that occurs up to distances from the heat source close to µ, where the heat is attenuated by a factor of 1/e. Consequently, in this particular case, the second choice of scales is the most natural one, being l=µ1. Choosing the characteristic length as µthe equation (8) reduces to: ¯ ∇2¯ T= 2∂¯ T ∂¯ t(9) Here it has been used that ω= 2πf and the definition of thermal diffusion length: µ=pα/(πf). For purely aesthetic reasons, the computational model has implemented the dimensionless heat equation with tc= 2/ω, so that the previous equation becomes: ¯ ∇2¯ T=∂¯ T ∂¯ t(10) The second equation from which the dimensions must be eliminated is the boundary condition associated with the laser. In this case the laser, of power P, amplitude modulated at a frequency f and centered at (x0,y0,0), is assumed to have a Gaussian profile and to be focused to a radius rg(at 1/e2of the maximum intensity) [6]. Therefore, the boundary condition in the stationary stage is: −κ∂T ∂z   z=0 =ηP πr2 g e −2hx−x0 rg2 +y−y0 rg2icos(2πft)(11) being κthe thermal conductivity of the sample material and ηthe power fraction absorbed by the sample. In this work, the negative sign is introduced as a phase in the modulation term so it will no longer appear again. On the other hand, ηis taken equal to one and a factor 2is introduced because when the right hand side of the equation is integrated over xand y, without having that factor, the result is P/2. Taking this three details into account, this equation can be rewritten as follows: κ∂T ∂z   z=0 =2P πr2 g e −2hx−x0 rg2 +y−y0 rg2icos(2πft)(12) This boundary condition becomes dimensionless using the previously introduced variables: (x, y, z)→(µ¯x, µ¯y, µ¯z), t →2 ω¯ t, T →¯ TT0(13) Therefore, applying the chain rule, the equation (12) becomes: 1In addition, this choice greatly simplifies the subsequent analysis of the problem based on the dimensionless parameters. 7 3 Theoretical model 3.1 Dimensionless equations ∂¯ T ∂¯z  ¯z=0 =2Pµ πr2 gT0κe −2h¯x−¯x0 rg/µ 2 +¯y−¯y0 rg/µ 2icos(2¯ t)(14) It must be noticed that the oscillation frequency of the laser heat source, in this dimensionless formulation, is no longer depending in any material or experimental parameters, being a constant value in any case. Once the dimensions of heat equation and laser boundary condition are removed, the two remaining equations are those associated with the crack. The constraints that must be satisfied are: • The continuity of heat flow over the crack. This condition is written as follows: [[ ˙ Q]] = 0 (15) where the [[ ]] operator stands for the change on the flux over the crack [6]. Being equal to zero, removing the dimensions of this equations is straight forward: [[ ˙ ¯ Q]] = 0 (16) • A temperature discontinuity at crack position: The crack is modeled as a thermal contact resistance Rth [6], related to the width of the crack wthrough: Rth =w κair (17) being κair the thermal conductivity of the air, which is assumed to fill the crack. Thus, the temperature discontinuity in the crack is given by: ∆T=κRth∇T=κw κair ∇T(18) Making use of the previously introduced dimensionless variables, this condition becomes: T0∆¯ T=κw κair T0 µ¯ ∇¯ T(19) Here the dimensionless gradient operator ¯ ∇≡∇/µ has been introduced. Rearranging the terms, the dimensionless temperature jump condition due to the presence of a crack results in: ∆¯ T=κ κair w µ¯ ∇¯ T(20) Summarizing, introducing the thermal diffusion length µas the length scale, tc= 2/ω as the characteristic time scale and the temperature normalization factor T0, the set of dimensional four equations that govern the physics of crack detection by lock-in IR thermography, becomes: 8 3 Theoretical model 3.2 Dimensionless parameters                ∇2T=1 α ∂T ∂t ∂T ∂z   z=0 =2P κπr2 ge −2hx−x0 rg2 +y−y0 rg2icos(2πft) [[ ˙ Q]] = 0 ∆T=κw κair ∇T →                ¯ ∇2¯ T=∂¯ T ∂¯ t ∂¯ T ∂¯z  ¯z=0 =2P µ πr2 gT0κe −2h¯x−¯x0 rg/µ 2 +¯y−¯y0 rg/µ 2icos(2¯ t) [[ ˙ ¯ Q]] = 0 ∆¯ T=κ κair w µ¯ ∇¯ T It has to be mentioned that this entire work deals with the adiabatic problem, that is, it is assumed that there are no heat losses due to convection or radiation mechanisms. This decision is based on the nature of this two mechanisms. Radiation is proportional to T4which means that if the temperature is not high enough its effect is negligible. On the other hand, heat losses by convection are proportional to the difference between the room temperature and the sample temperature. However, IR thermography experiments are carried out a few kelvins above the room temperature so the difference is not large enough to be considered. Furthermore, this simplification is supported by the fact that, except in a few cases, the adiabatic model fits the experimental data relatively well. 3.2 Dimensionless parameters Once the non-dimensional equations have been obtained, the terms that appear in the new equations can be rearranged to find characteristic independent dimensionless parameters which determine the nature of the problem. In particular, the combination of the terms that leads to these parameters can be found looking at the laser boundary condition (14) and the temperature jump condition at the crack (20). In both of them thermal conductivity κappears, so, multiplying and dividing the boundary condition of the laser by the conductivity of the air, this equation becomes: ∂¯ T ∂¯z  ¯z=0 =2Pµ πr2 gT0κair 1 κ κair e −2h¯x−¯x0 rg/µ 2 +¯y−¯y0 rg/µ 2icos(2¯ t)(21) Going further, in this equation the radius of the laser rgappears both in the amplitude and in the exponential terms. So, it seems to be natural to multiply and divide the dimensionless modulation amplitude by the thermal diffusion length: ∂¯ T ∂¯z  ¯z=0 =2P πµT0κair 1 κ κair 1 rg µ2e −2h¯x−¯x0 rg/µ 2 +¯y−¯y0 rg/µ 2icos(2¯ t)(22) Defining the following dimensionless parameters: Π1≡2P πµT0κair ,Π2≡rg µ2,Π3≡κ κair (23) The laser boundary condition reduces to: 9 3 Theoretical model 3.2 Dimensionless parameters ∂¯ T ∂¯z  ¯z=0 =Π1 Π2Π3 e −2 Π2h(¯x−¯x0)2+(¯y−¯y0)2icos(2¯ t)(24) On the other hand, defining the dimensionless parameters associated to the crack geometry, Πw≡w µ,Πl≡l µ,Πd≡d µ(25) being land dthe length and depth of the crack, respectively, the temperature jump condition becomes: ∆¯ T= Π3Πw¯ ∇¯ Tif ¯x∈h−Πl 2,Πl 2i,¯z∈h0,Πdi(26) Therefore, it is observed that, from the manipulation of two of the four equations and the definition of six dimensionless parameters, three associated with the geometry of the crack and three with the experimental conditions, a simplified non-dimensional formulation of the original equations is obtained. The equations to compute are summarized in the following box.              ¯ ∇2¯ T=∂¯ T ∂¯ t ∂¯ T ∂¯z  ¯z=0 =Π1 Π2Π3e −2 Π2h(¯x−¯x0)2+(¯y−¯y0)2icos(2¯ t) [[ ˙ ¯ Q]] = 0 ∆¯ T= Π3Πw¯ ∇¯ T (27) In this work in addition to the parametric analysis associated to the mentioned six parameters, the position of the laser (¯x0,¯y0)and the inclination of the crack θwill also be considered in the analysis in order to check their impact in the results. IR thermography experiments in laboratories are typically carried out with lasers of power on the range of 0.1−10 W, radius at 1/e2of the maximum intensity is typically on the order of 10−4m and are modulated at frequencies on the order of hertz. The cracks to be detected in these experiments are typically on the order of microns. As said before, in this experiments the interest is placed on the oscillation of the temperature so, the most natural choice of temperature scale T0is the maximum of the amplitude over the room temperature. Typically this values are on the order of few kelvins. Considering the mentioned variation ranges, table 1 shows typical values of the dimensionless parameters. In this table it can be seen that the values of the first dimensionless parameter Π1can vary from 103to 104while values on the order of 10−4or 10−1are obtained for Π2. A similar behaviour is found for Π3, where values from 50 to almost 16000 are obtained. This leads to conclude that these three parameters are highly dependent on the material, which is a direct consequence of choosing the diffusion length as longitude characteristic scale. However, the same value of any dimensionless parameter can be obtained for different materials through the appropriate selection of the rest of experimental parameters. 10 3 Theoretical model 3.2 Dimensionless parameters κ(Wm−1K−1) [15] α(mm2s−1) [15] Π1Π2Π3 Cu 397.48 116.0 4685 0.0054 15899.2 AISI 304 14.64 3.68 28210 0.1963 600.0 Pb 34.309 23.3 10454 0.0270 1372.4 Lead Glass 1.13 0.74 58661 0.8490 45.2 Al 225.94 91.0 5289 0.0007 9037.6 Fe 71.965 20.4 11172 0.0307 2878.6 Table 1: Typical values of the dimensionless parameters for different materials. These values have been calculated with P= 1 W, f= 5 Hz, rg= 0.2mm, T0= 2 K. 11 4 Computational numerical simulations 4 Computational numerical simulations There are two key aspects in computational numerical simulations: the mesh and the numerical method that is used to solve the equations. The manner the equations are translated from analytical formulation to computational formulation is an important aspect when using finite element methods because some numerical schemes may not work for some specific problems. In the following a brief description is these is given, even though it is not the scope of this work to go on further detail in this topic. On the other hand, the mesh plays a key role in solving equations by numerical methods because it is, literally, where the equations are solved. This means that, in this case the simulated laboratory sample, has to be as similar as possible to a continuous medium but without having to spend much computational resources and time simulating the problem. The simulations presented in this work have been carried out a workstation with Ubuntu MATE 20.04.4 with an Intel Xeon(R) Gold 5218 CPU @ 2.30 GHz ×64, 192 Gb memory and a graphic memory LLVM 12.0.0, 256b. 4.1 OpenFoam Open Source Field Operation and Manipulation (OpenFoam) [16] is a C++ object oriented library, originally designed for computational fluid dynamics and structural analysis [17, 18]. However, after decades of evolution, it has become a software with multiphysics capacity oriented to a wide variety of physical phenomena, such as: combustion, electromagnetics, heat transfer and others. OpenFoam is used to create executables that fall in two categories: those that allow the manipulation of data, known as utilities, and solvers, which are designed to solve a specific problem in continuum mechanics [16]. Taking into account the nature of the problem to be solved in this work, among the variety of solvers that are implemented in OpenFoam, the solvers of the thermal family have been identified as the most suitable. Between all the options available, solidFoam has been selected as it has been designed for energy transport and thermodynamics on solids. In this solver not only complex boundary conditions, such as the ones of the laser or the crack, can be implemented but it could also benefit a future work in which this model is used to simulate flying spot experiments [19] because it allows to use dynamic meshes. It is worth to mention that OpenFoam presents native parallel calculation capacity, which allows to distribute the nodes of the mesh between the available processes. Hence, this presents a double benefit. On the one hand, several simulations can be done simultaneously. On the other hand, faster calculations can be performed assigning each process a set of nodes. 4.2 Numerical schemes As mentioned before, it is beyond the scope of this work to explain in detail the methods used to calculate the partial derivatives, gradients, etc. since the objective is to build the dimensionless model and these methods are equally valid for equations with and without dimensions. However, it should be mentioned that the schemes used have already been used to carry out simulations in similar physical phenomena with proven capabilities. The finite volume schemes that have been used in this work are 12 4 Computational numerical simulations 4.3 Solution and algorithm control summarized in table 2. Quantity Scheme Time partial derivatives Crank Nicholson Gradients Gauss linear Laplacian Gauss linear corrected Interpolations Linear Table 2: Finite volume schemes used. 4.3 Solution and algorithm control In solution and algorithm control, the solver keyword specifies each linear-solver that is used for each discretised equation, that is, to the method of solving the set of linear equations [16]. On the other hand, the multiple options for preconditioning of matrices in the conjugate gradient solvers (DIC, FDIC, DILU,...) are controlled by the preconditioner keyword [16]. In this work the PCG solver and DIC preconditioner have been used as it has been seen that they suit properly to the characteristics of the problem. In OpenFoam, the matrix solvers are based on reducing the error in the solution over an iterative process, that is, the residual is evaluated by substituting the current solution into the equation and taking the magnitude of the difference between the left and right hand sides [16]. In order to control if that difference is small enough, there are two variables that can be used: 1. Absolute tolerance: measures if the residual is small enough to consider the solution sufficiently accurate. 2. Relative tolerance: limits the relative improvement from initial to final solution. Thus, the solution will be considered sufficiently accurate if the residual is lower than the absolute tolerance, or the ratio of current to initial residuals falls below the relative tolerance. In this work the absolute tolerance has been used as a measure of the accuracy of the simulation and it has been set to 1×10−7. Another aspect to be taken into account is the solution under-relaxation. This is a technique used for improving stability of a computation that can be applied in physical simulations where variables vary that fast that can lead to numerical divergences. Since the nature of the experiment that is being simulated does not imply this kind of problems, the relaxation factor in this work is left as default, that is, equal to 1. The solution and algorithm control configuration used in the simulations is shown in table 3. Solver Preconditioner Absolute tolerance Relative tolerance Relaxation factor PCG DIC 1×10−70 1 Table 3: Solution and algorithm control keywords used. 13 4 Computational numerical simulations 4.4 The mesh 4.4 The mesh As said before, the mesh plays a key role in computational numerical simulations. There are two main strategies to simulate cracks in materials. The first one is to model the problem as a two domain material: the bulk and the air filling the crack. This strategy implies that, close and inside the crack, an extremely fine mesh has to be done which dramatically increases memory resources and calculation time due to the difference in spatial scale between the crack and the bulk [20, 21]. The second option is to model the crack as a contact thermal resistance surface. In this approach, there is no crack volume to mesh, because the latter is modeled as a 2D plane [2, 6], notably reducing the computational resources required. In OpenFoam there are two main options to model a crack. The first one consists in using the function already implemented in OpenFoam, which allows to create the sample as an entire block and then define the geometry of the crack. This method allows the user to simulate as many cracks as wanted without adding complexity to the mesh. However, this strategy has a big drawback: it implies that the geometry of the crack will be rounded to the values to the nearest node of the mesh, meaning that the crack will no longer be as wide, long or deep as it had been modeled. The second option is to reproduce the crack as a contact surface between hexahedral blocks, which implies that the sample has to be divided into multiple domains that later have to be joined. In other words, the sample has to be constructed ‘block by block’. Since the number of blocks to create depends entirely on the type of crack (infinite, semi-infinite or finite), the main disadvantage of following this path is that the complexity of the code increases. On the other hand, following this methodology ensures, by construction, that the crack will have the desired dimensions. As a consequence, in this work, this has been the followed meshing strategy. Another key aspect to be taken into account is the geometry of the mesh. Among all the options, the most used meshes are tetrahedral and hexahedral meshes. These two types of meshes and the individual elements can be seen in figures 4 and 5. Figure 4: (left) Hexahedral and (right) tetrahedral meshes. Figure 5: (left) Tetrahedral and (right) hexahedral mesh elements. Tetrahedral meshes are made up of elements with 4 vertices. These can be rearranged so that a 14 4 Computational numerical simulations 4.5 Mesh resolution control 4.5.2 Mesh refinement Another method to increase the mesh resolution in some regions in OpenFoam is to divide the hexahedrons that are in the region of interest. Unlike the gradient method, mesh refinement increases the computational resources because the number of nodes increases. The implementation of the refinements in OpenFoam consists of: 1. Selecting the nodes of the mesh that are within the region of interest. With this objective, geometrical entities (spheres, cylinders, boxes,...) are defined in the mesh inside which the selected nodes remain. 2. Once the nodes of the region of interest are selected, these are divided until the desired spatial resolution is achieved. As mentioned before, the interest of implementing mesh refinements is that there are regions that are physically more interesting, since they are the ones in which the variables change more rapidly. In particular, during the multiple simulations that have been carried out in this work, several special interest zones have been identified: the crack and the laser spot. In order to mesh these regions with enough spatial resolution, a region proportional to the non-dimensional thermal diffusion length has been carefully refined. With the objective of refining the region close to the crack, the box geometry has been used as it keeps the aspect ratio with the sample. Furthermore, during the simulations that have been carried out, it has been noticed that a single box refinement not only increases significantly the resolution in that region but also reduces both the number of nodes far from the crack and consequently, the computational time. An example of this refinement can be found in figure 13 (a). In addition to the meshing of the crack or the region close to it, an important aspect in IR thermography simulations is to correctly mesh the region in which the laser is introduced. In order to complete this task, the cylinder refinement has been chosen as it keeps the aspect ratio with the circular geometry of the laser spot. During the multiple simulations done in this work, it has been tested that refining the cells of the laser environment three times, each time with cylinders of smaller height and radius, it is enough to capture the laser spot appropriately. This triple refinement can be seen in figure 13 (b). When both of these refinements are combined (see figure 13 (c)), as it is the case of this work, the resulting mesh ends up having on the order of 250000 nodes. (a) (b) (c) Figure 13: (a) Single box, (b) triple cylinder refinements and (c) both of them combined. 21 4 Computational numerical simulations 4.6 Processing 4.6 Processing Once the mesh is done, the equations must be solved. At this stage of the calculation, known as processing, it is worth to notice that even though only the stationary part of the problem is being simulated, for numerical reasons, it may happen that the first oscillations do not represent that part correctly until a few iterations are done. Therefore, in order to find the stable stationary situation without undesired numerical effects a criterion that guarantees the appropriate stationary nature of the calculation is required. In this work this criterion has been selected as the full reproducibility of the obtained thermal waves. Programmatically, the implementation of this condition has been carried out selecting an arbitrary point of the stationary thermal cycle and checking its value over cycling, provided an appropriate time resolution. When the temperature difference between cycles, i.e., the slope of the line joining two consecutive points (as can be seen in figure 14) is found to be below an imposed threshold value, the criterion is satisfied and the next cycle is saved as the stationary oscillation for each node of the mesh. Figure 14: Temperature evaluation of an arbitrary point over cycling until the stationary stability criterion is satisfied. 4.7 Post-Processing As mentioned before, one of the consequences of treating this problem in a simple way, from the mathematical point of view, is that instead of solving the Helmoltz equation [4] which integrates the harmonic nature of the lock-in experiment, the heat equation is solved. Whereas the first provides direct solution for the thermal amplitude, the second doesn’t, resulting the temperature field. The consequence of this is that, after the processing stage, those amplitudes must be found. This is the post-processing stage of the calculation. In this stage, even if the temperature field for the complete sample is calculated, the results presented in this work are limited to a line perpendicular to the crack and crossing the center of the laser spot. This selection has been done as this would be the physically more significant region. Once the region to plot the non-dimensional temperature amplitudes is chosen, the maximum and minimum values of each thermal wave are programmatically searched, leading to the non-dimensional temperature amplitude plot. 22 5 Results 5 Results 5.1 Validation of the model Once the dimensionless numerical model has been established, its validity must be checked. Although the model can be validated against experimental data, in this case it is going to be done against analytical or numerical models, depending on the situation. In this section first, the comparisons of the results of the dimensionless model with the analytical model corresponding to a vertical infinite crack [1] will be presented in figure 15. Second, the non-dimensional results for semi-infinite and finite cracks are going to be validated against another numerical model developed by R.Celorrio et al. [4] (in the following, DG FEM model) written in FEniCS, due to the fact that there is no analytical solution in these cases. The associated results in this case are shown in figures 16 and 17. In order to make the cases as realistic as possible, it has been decided to validate the model with typical experimental parameter values (see table 4) and for two materials: AISI 304 and Cu (see their thermal properties in table 1). P(W) rg(mm) f(Hz) (x0, y0)(mm) 1 0.165 1 (0,0.5) Table 4: Dimensional parameter values that have been used in order to validate the dimensionless numerical model. Aiming to compare the results between models, the dimensional results are transformed to dimensionless following the scale variables previously introduced (see equation (13)). In order to maintain an acceptable sensitivity of the temperature changes, even far from the heat source, instead of the amplitude its natural logarithm is plotted. The dimensionless parameters that have been used can be found in table 5. Figure 15: Comparison between the non dimensional numerical results and analytical solution for infinite vertical cracks. (a) AISI 304 and (b) Cu. 23 5 Results 5.1 Validation of the model Figure 16: Comparison between the non dimensional results and the DG FEM model in the semiinfinite crack case in AISI 304 for two different crack inclinations. (a) 90º and (b) 45º. Figure 17: Comparison between the non dimensional results and the DG FEM model in the finite crack case in AISI 304 for two different crack inclinations. (a) 90º and (b) 45º. In order to quantify the differences between the dimensional (analytical or numerical) and dimensionless models, among all the statistical options, in this work the root mean square error (RMSE) has been chosen due to the ease of interpreting the results. As it can be seen in figures 15, 16 and 17 the agreement between all the models is very satisfactory in the three presented cases, the values presented in table 6, are going to be discussed in detail in section 6.1. 24 5 Results 5.2 Parametric analysis Figure 15 (a) Figure 15 (b) Figure 16 Figure 17 Type of crack Infinite Infinite Semi-Infinite Finite Π123528 4191 23528 23528 Π20.023 0.0007 0.023 0.023 Π3586 15899 586 586 Πw0.00092 0.0082 0.00092 0.00092 Πl- - - 0.92 Πd- - 0.92 0.92 ¯y00.46 0.45 0.46 0.46 θ(º) 90 90 90/45 90/45 Table 5: Dimensionless parameters used to obtain the results shown in figures 15, 16 and 17. Case / Model Analytical DG FEM model Infinite 90º (I) (figure 15 (a)) 0.017 - Infinite 90º (II) (figure 15 (b)) 0.014 - Semi-infinite 90º (figure 16 (a)) - 0.009 Semi-infinite 45º (figure 16 (b)) - 0.012 Finite 90º (figure 17 (a)) - 0.0027 Finite 45º (figure 17 (b)) - 0.022 Table 6: RMSE of the comparisons shown in figures 15, 16 and 17. 5.2 Parametric analysis In this section a catalog of curves in which each non-dimensional parameter is varied, while the others are fixed, is shown. This catalog has two main purposes. First, it can be used as a result guideline for a variety of experimental and material properties parameters. Going more in detail, the second purpose of this catalog, is to provide a parametric analysis showing the influence of each non-dimensional parameter in the resulting thermal amplitude (|¯ T|) plot. Although it is a common practice to normalize the experimental data in order to eliminate the effect of parameters which are difficult to control, such as the absorption coefficient (η) in equation (11), the normalization can also limit the amount of information obtained in the experiment. In order to take advantage of both methodologies, in this work, as general rule, both normalized and nonnormalized plots are shown2. All of the results presented in this section are obtained for a finite crack, as it is the most versatile one. The experimental parameters that lead to these simulations are typical values of experimental setups. It is worth to notice that the power of this non-dimensional formulation lies in the fact that the result is not dependent on each individual experimental parameter, but depends only on their combinations. The dimensionless parameters that remain fixed in each calculation are specified in table 7. 2In the results in which the normalization does not have any impact, only the normalized natural logarithm of the thermal amplitude ln(|¯ Tn|)is shown. 25 5 Results 5.2 Parametric analysis Figure 18: Natural logarithm plot of the thermal amplitude (a) non-normalized and (b) normalized on the transverse sample profile for Π1= 4 ×103,1×105,2×105,3×105. Figure 19: Natural logarithm plot of the thermal amplitude (a) non-normalized, showing the zoomed amplitude jump associated to the crack in the insert and (b) normalized on the transverse sample profile for Π2= 1 ×10−7,2×10−7,3×10−7,1×10−8. 26 5 Results 5.2 Parametric analysis Figure 20: Natural logarithm plot of the thermal amplitude (a) non-normalized and (b) normalized on the transverse sample profile for Π3= 50,500,1000,4000. Figure 21: Natural logarithm plot of the thermal amplitude (a) non-normalized and (b) normalized on the transverse sample profile for ¯y0= 0.24,0.32,0.40,0.48. 27 5 Results 5.2 Parametric analysis Figure 22: Natural logarithm plot of the thermal amplitude (a) normalized and (b) crack region zoomed on the transverse sample profile for Πw= 0.001,0.002,0.005,0.01. Figure 23: Natural logarithm plot of the thermal amplitude (a) normalized and (b) crack region zoomed on the transverse sample profile for Πl= 0.2,0.4,0.5,0.75. 28 5 Results 5.2 Parametric analysis Figure 24: Natural logarithm plot of the thermal amplitude (a) normalized and (b) crack region zoomed on the transverse sample profile for Πd= 0.1,0.2,0.3,0.5. Figure 25: Natural logarithm plot of the thermal amplitude (a) normalized and (b) crack region zoomed on the transverse sample profile for θ= 50º, 75º, 90º, 100º. 29 5 Results 5.2 Parametric analysis Figure Π1Π2Π3¯y0ΠwΠlΠdθ(º) 18 - 0.027 600 0.50 0.001 1 1 50 19 25231 - 600 0.50 0.001 1 1 50 20 25231 0.027 - 0.50 0.001 1 1 50 21 25231 0.027 600 - 0.001 1 1 50 22 25231 0.027 600 0.50 - 1 1 50 23 25231 0.027 600 0.50 0.001 - 1 50 24 25231 0.027 600 0.50 0.001 1 - 50 25 25231 0.027 600 0.50 0.001 1 1 - Table 7: Parameters used to obtain each figure of the catalog. 30 7 Conclusions 7 Conclusions In this work a non-dimensional FEM model for crack characterization in aerospace materials by means of lock-in IR thermography has been developed. In addition to the resolution of the equations, the optimization of the spatial domain to be modeled has been carried out through different meshing strategies, such as non-regular edge distribution or selective refinements of the mesh. These meshing strategies have optimized the calculations in terms of accuracy and computational resources. The developed model shows a very good agreement when compared with other analytical or numerical models, although the formulation of each one has a slight impact on the results. Overall, the comparisons show very low relative errors (at most 0.4%) which confirm the remarkable accuracy of the developed model. The dimensionless model has allowed a general interpretation of the lock-in IR thermography experiment not depending on the parameters of the material or the experimental ones, but on the dimensionless parameters identified. These parameters have been used to perform a parametric analysis where their effect in the thermal amplitude plots has been analyzed. While the effect of Π1is removable by normalization, the effect of Π3cannot be isolated in just one zone of the curve since it affects the morphology of the entire thermal amplitude plot. On the other hand, despite of different morphological reshaping of the curve, the effect of Π2and ¯y0has turned out to be unique, i.e., it is impossible to obtain the same thermal amplitude plot with different values of these. It has been seen that, although the most prominent effect of Πw,Πl,Πdand θis to change the amplitude jump, they also show other secondary effects in other regions such as the slopes far from the heat source. Moreover, except for θ, in the other three parameters those secondary effects can be seen more clearly on the non-illuminated side of the crack. Even though it is beyond the scope of this work, the results obtained for the parametric analysis carried out allow to determine where the sensitivity of the curves is to each of the parameters. This becomes of paramount importance when parametric inversion is performed. Going further in the developed non-dimensional discussion, the degeneracy of the thermographic problem has been addressed. In this line, a non-dimensional reformulation of the problem has been presented. It has been demonstrated that there is no need to distinguish between Π1,Π3and Πw, since different combinations of these parameters lead to the same results. This work has allowed to determine which combinations of experimental parameters and material properties leading to the same amplitude thermograms which are hardly accessible by other procedures. As a potential continuation of the developed investigation the following future work is identified: 1. The study of the degeneracy of the non-dimensional formulation searching other potential parametric combinations. 2. The development of inverse models which make use of the non-dimensional formulation presented in this work. 3. The introduction of more general crack geometries beyond planar. 37 References REFERENCES References [1] N. Pech-May, A. Oleaga, A. Mendioroz, A. Omella, R. Celorrio, and A. Salazar, “Vertical cracks characterization using lock-in thermography: I infinite cracks,” Measurement Science and Technology, vol. 25, no. 11, p. 115601, 2014. [2] R. Celorrio, A. Omella, N. Pech-May, A. Oleaga, A. Mendioroz, and A. Salazar, “Vertical cracks characterization using lock-in thermography: Ii finite cracks,” Measurement Science and Technology, vol. 25, no. 11, p. 115602, 2014. [3] S. Findlay and N. Harrison, “Why aircraft fail,” Materials Today, vol. 5, pp. 18–25, 11 2002. [4] R. Celorrio, A. Omella, A. Mendioroz, A. Oleaga, and A. Salazar, “Advances in crack characterization by lock-in infrared thermography,” International Journal of Thermophysics, vol. 36, no. 5, pp. 1202–1207, 2015. [5] S. Sfarra, E. Marcucci, D. Ambrosini, and D. Paoletti, “Infrared exploration of the architectural heritage: from passive infrared thermography to hybrid infrared thermography (hirt) approach,” Materiales de Construcci´ on, vol. 66, no. 323, pp. e094–e094, 2016. [6] J. Rodr´ ıguez-Aseguinolaza, M. Colom, J. Gonz´ alez, A. Mendioroz, and A. Salazar, “Quantifying the width and angle of inclined cracks using laser-spot lock-in thermography,” NDT & E International, vol. 122, p. 102494, 2021. [7] D. P. Almond and S. Lau, “Defect sizing by transient thermography. i. an analytical treatment,” Journal of Physics D: Applied Physics, vol. 27, no. 5, p. 1063, 1994. [8] J. G. Sun, “Analysis of Pulsed Thermography Methods for Defect Depth Prediction,” Journal of Heat Transfer, vol. 128, no. 4, pp. 329–338, 10 2005. [Online]. Available: https://doi.org/10.1115/1.2165211 [9] J. Laskar, S. Bagavathiappan, M. Sardar, T. Jayakumar, J. Philip, and B. Raj, “Measurement of thermal diffusivity of solids using infrared thermography,” Materials Letters, vol. 62, no. 17-18, pp. 2740–2742, 2008. [10] A. A. Badghaish and D. C. Fleming, “Non-destructive inspection of composites using step heating thermography,” Journal of composite materials, vol. 42, no. 13, pp. 1337–1357, 2008. [11] M. Colom, J. Rodriguez-Aseguinolaza, A. Mendioroz, and A. Salazar, “Imaging real cracks: evaluation of the depth and width of narrow fatigue cracks in and al-alloys using laser-spot lock-in thermography,” in Thermosense: Thermal Infrared Applications XLIII, vol. 11743. International Society for Optics and Photonics, 2021, p. 117430F. [12] G. Nikishkov, “Introduction to the finite element method,” University of Aizu, pp. 1–70, 2004. [13] ´ A. J. O. Mili´ an, “Resoluci´ on num´ erica de problemas inversos en termograf´ ıa activa para la caracterizaci´ on de grietas,” Ph.D. dissertation, Universidad de Zaragoza, 2017. [14] H. P. Langtangen and G. K. Pedersen, “Basic partial differential equation models,” in Scaling of Differential Equations. Springer, 2016, pp. 69–98. [15] Thermtest instruments, “Materials Thermal Properties Database,” https://thermtest.com/ thermal-resources/materials-database, Accessed: 9/5/2022. [16] OpenFoam, “OpenFoam,” www.openfoam.org, Accessed: 11/6/2022. 38 References REFERENCES [17] H. G. Weller, G. Tabor, H. Jasak, and C. Fureby, “A tensorial approach to computational continuum mechanics using object-oriented techniques,” Computers in physics, vol. 12, no. 6, pp. 620–631, 1998. [18] H. Jasak, A. Jemcov, Z. Tukovic et al., “Openfoam: A c++ library for complex physics simulations,” in International workshop on coupled methods in numerical dynamics, vol. 1000. IUC Dubrovnik Croatia, 2007, pp. 1–20. [19] M. Colom, A. Bedoya, A. Mendioroz, and A. Salazar, “Measuring the in-plane thermal diffusivity of moving samples using laser spot lock-in thermography,” International Journal of Thermal Sciences, vol. 151, p. 106277, 2020. [20] N. J. Wallace, N. B. Crane, and M. R. Jones, “Defect measurement limits using flash thermography with application to additive manufacturing,” NDT & E International, vol. 128, p. 102615, 2022. [21] M. Krishnapillai, R. Jones, I. H. Marshall, M. Bannister, and N. Rajic, “Ndte using pulse thermography: Numerical modeling of composite subsurface defects,” Composite Structures, vol. 75, no. 1-4, pp. 241–249, 2006. [22] R. Biswas and R. C. Strawn, “Tetrahedral and hexahedral mesh adaptation for cfd problems,” Applied Numerical Mathematics, vol. 26, no. 1-2, pp. 135–151, 1998. 39