From Thermal Inspection to Updating a Numerical Model of a Race Bicycle : Comparison with Structural Dynamics Approach

Carbon fiber bicycle frames are complex-shaped structures and are prone to delaminations and difficult to inspect. The use of finite element model updating is common in structural dynamics but not so common in active thermography inspection. However, there are many advantages to using thermography when inspecting bicycle frames. These include the fact that the inspection can be performed in situ, can cover large areas, and is a quantitative method. In this paper, a numerical model of a bicycle frame will be updated and optimized by the surface temperature distribution captured with pulsed thermography. These results will be compared and benchmarked against frequency response function (FRF) measurement data as a reference. The chosen temperature decay measurements to be used as reference data will be of key importance. The goal of this manuscript is to compare both measurement results and model predictabilities after performing finite element model updating with respect to accuracy and speed.


Introduction
The application of composite materials in performance components used daily such as bicycle frames has increased due to its beneficial mechanical properties with respect to metals.Such materials are strong, lightweight, and fatigue-and corrosion-resistant when combined with the proper resins [1].Composite materials are technically complex structures due to their anisotropic behavior, which makes quality control particularly challenging.Thus, structural health monitoring is essential in the production and lifetime stages of such composite structures.This work is concerned with carbon bicycle frames [2].To detect and analyze defects and guarantee an extended lifetime, composite materials have been the subject of a tremendous amount of non-destructive inspection (NDI) research [3,4].Non-destructive evaluation (NDE) techniques are broadly used to evaluate material properties and to inspect the integrity of structures without damaging the structure [5].Infrared (IR) thermography has remarkable advantages with respect to conventional NDE techniques, such as ultrasound C-scanning, X-ray radiography, tap-testing, and visual testing [6].The technology is contactless, safe, fast, and reliable [7].The study of active thermography has become an important NDE technique for damage detection and material updating in metallic structural components [6,8,9] and in carbon fiber reinforced polymer (CFRP) and glass fiber reinforced polymer (GFRP) composites [10].As the application of composite laminates such as CFRPs and GFRPs increases, the thermal conductance inside the components becomes anisotropic, and quantitative prediction becomes complex.The common form of active thermography makes use of heat pulses emitted to the surfaces and is called pulse thermography (PT) [11].PT is used not only to identify defects but also to determine the thermal diffusivity of the material, inspect the structural integrity, and monitor the production process of composite structures.
High-performance carbon fiber bicycle frames are strong and lightweight.However, they are susceptible to damage caused by low-energy impact loading, which have many possible causes: small stones can be projected by nearby cars or the front wheel of the bicycle, the bicycle may fall and hit the ground, and defects can occur during manufacturing [12,13].In addition, material thickness inconsistency can occur during manufacturing and cause false positive defect hits.
Finite element (FE) modeling is common in the design of large structures but generally remains in the design stage.Two of the largest drawbacks of the use of FE beyond the design stage are the approximated material property parameters and the high computational cost of building an accurate FE model of the full structure.Finite element model updating (FEMU) can deliver a solution to these two problems, whereby the model is updated with non-destructive testing data, mainly done with ultrasound vibration C-scanning [14][15][16].The foundation of the technique has been described by Friswell and Mottershead in 1995 [14].
The main idea is to solve a set of initial numerical models with variable parameter values for the uncertain properties.By comparing the resulting response values (mostly done on vibration modes and/or eigenfrequencies) with experimental measurements, the variable parameters are modified iteratively until the response values correlate with the experimental verification data.Especially for complex-shaped objects and non-linear behavior, the methodology delivers faster and more accurate results with the inclusion of higher-order terms in the equations [14].
Unfortunately, the computational cost of an FEMU routine is even more expensive than that of a regular numerical simulation, as for each iteration at least one new numerical model needs to be computed.Therefore, a methodology based on the technology developed by Steenackers has been used to replace the FE model with an approximative regression model in order to reduce the calculation times [17].When the relationship between the input parameters and the requested output information on the mesh nodes can be modeled with significant accuracy, a major part of the time-consuming FE computations can be bypassed during the model updating routine by implementing a so called "metamodel."Some further computational intelligence techniques have been used within FEMU routines for structural dynamics and compared by Marwala in 2010 [15].Some general challenges remain.Due to inaccuracy in the model and to measurement errors, there is a danger that wrong assumptions are made.Proper validation, physical understanding, and adequate modeling are essential to verify the value of the performed FEMU routines.
These updating methods are widely used in structural dynamics for modal vibration inspection [14,[18][19][20][21][22], but their use is limited in heat transfer optimization problems related to non-linear convection and radiation boundaries [23,24].FE model updating can be used to optimize the complex parameter set of active thermographic setups by correlating the numerical results with NDT data [25].These methods are widely used in modal vibration inspection [20][21][22], but their usage is limited in thermal optimization problems [4,23,24] due to the difficulties in estimating the diffusivity.
In this paper, we compare 'easy-to-perform' FEMU on thermal non-destructive measurements, with respect to vibrational analysis, with regular finite element model updating.We therefore limit our experimental measurements to non-contact measurement techniques from a single view point.In practice, it is generally difficult to measure the same structure from different positions around the structure.The goal is to evaluate the influence of this limitation on the measurement results in order to evaluate local material thickness.

Thermal Experiments
A carbon fiber reinforced plastic (CFRP) bicycle frame was inspected in order to estimate the tube thickness and inspect for defects such as cracks and delaminations using an MWIR FLIR Phoenix thermal camera (3-5 µm, 640 × 512 pixels, 25 mm lens) and an optical excitation source consisting of two 3.2 kJ Balcar flash excitation units.The bicycle was painted in black and yellow, and the tubular shapes and metallic inserts made evaluating the structure challenging.A representation of the test setup is shown in Figure 1.

Structural Dynamics Experiments
The same CFRP bicycle frame was evaluated by measuring frequency response functions (FRFs) and using state-of-the-art modal analysis techniques with a PSV 400 Polytec scanning laser doppler vibrometer system and an electrodynamic modal shaker.The bicycle frame was suspended freely, approximating free-free boundary conditions, and accelerations were measured on different positions over the structure, which is shown in Figure 2. From the measured FRFs, eigenfrequencies and corresponding mode shapes were estimated as reference data to be used for modal updating.In total, the setup contained 290 measurement points.Within all points, the deflection along 1 degree of freedom (DOF) was measured.With respect to the optimal excitation direction, we used the DOF perpendicular to the left side of the bicycle frame.FRF measurements were performed in free-free boundary conditions, using a swept sine excitation over a frequency band of 1-6000 Hz with a resolution of 0.5 Hz, and the modal parameters (modal frequencies, damping, and mode shapes) were estimated [26].

Numerical Simulation
The finite element model of the CFRP bicycle frame was built from a 3D scan measured with a Creaform 3D scanner with an accuracy of less than 0.1 mm using the structured light principle.The resulting .stlwas cleaned and used to build a 2D structured mesh of CQUAD elements within the Siemens Simcenter 11 software (Siemens PLM, Plano, Texas, US), and an overall element size of 5 mm was used.The tubes were formed by an extruded 3D mesh from the 2D mesh basis.A representation of the 2D mesh is shown in Figure 3.The laminated properties of the bicycle frame were approximated using a five-layered laminate consisting of material properties shown in Table 1:   A symmetric quasi-isotropic laminate was represented as the stacking recipient.The precise stacking and layer properties are unknown, so it was assumed that the represented stacking was sufficiently accurate as an input for the material parameter optimization process, where the global thermal diffusivity properties in the thermal case and the mode shapes in the structural case are updated based on reference measurements.A detailed picture of the 3D mesh with the different layers is shown in Figure 4.The mesh consisted of 28,623 elements with a structured material orientation according to the laminate properties.The metal inserts were neglected, as they were not positioned in the vicinity of the evaluated measurement areas.In contrast to the mesh generation, numerical analysis was performed using Comsol Multiphysics 5.0 to facilitate interaction with the optimization routine using Matlab.
The excitation source was modeled as a combination of two point sources.This excitation source results in a non-homogeneous spatial distributed heat excitation.The thermal time constant results in, aside from a spatial non-homogeneity, a different distribution related to the type and history of the heat source, which needed to be considered in the numerical model.The used flash excitation units are visualized in Figure 5.It can be clearly seen that the excited radiation is delayed from the initiated trigger signal due to the thermal inertia.To correlate the numerical model accurately with the experimental data, this thermal inertia will be considered as a correction (weighting) of the input signal with the thermal response.

Methodology
Both the quality of the experimental data and the FE model data will contribute to the evolution of the optimization algorithm.The implemented optimization algorithm is based on the adaptive response surface method described in [24].With respect to modal analysis updating, the original optimization routine is used, developed by Steenackers et al. [27].In this section, we will start with a short description of the experimental data acquisition and post-processing of the thermal and modal data.The numerical model is briefly described in Section 2.3, which ends with a brief description of the used optimization routine for finite element model updating.

Thermal Experiments
The thermal experiments were performed using pulsed thermography (PT), schematically represented in Figure 6, via a flash excitation source and a thermal camera.This active-thermography NDE technique is based on an analysis of the response of the surface temperature of the structure to an emitted thermal wave initiated by an excitation pulse [28].It is one of the most popular methods in active thermography due to its simplicity and fast operation.The thermal data were stored for several time steps.A 3D matrix with the XY-direction containing the temperature of each pixel and the Z-direction capturing each time frame was delivered [4,29].
Local material properties, thickness, or experimental set-up parameters, such as excitation power and distance, could be defined [29].The measurement accuracies were dependent on ambient temperature, emissivity, reflecting temperature, and humidity.These parameters were fixed when measurements were taken.An empirical rule of thumb was used to determine the minimally detectable diameter/depth ratio (should be equal to 2) [30].
Measurements were post-processed using pulsed phase thermography (PPT), which is a popular signal transformation technique for PT measurements.Considering its relation with thermal wave theory, as described in [7], the (dirac) pulse excitation of a PT measurement can be decomposed into a multitude of individual sinusoidal waves with different amplitudes and frequencies [31].The benefit of a short dirac equivalent pulse is found in the flat frequency distribution of the discrete Fourier transformed equivalent spectrum.This method is a combination of the advantages of lock-in thermography with the speed of PT.As the dirac pulse is a combination of multiple sine waves, a fast examination can be performed for multiple depths through the structure.By computing the blind frequency f b and using Equation (1), it is possible to quantitatively estimate the defect depth.
f b is the blind frequency, and C 1 is a constant value typically equal to 1.82 [32].
The maximum penetration depth in a material is defined by the thermal diffusion length µ = α π• f b , which can be translated to the measured depth according to Equation (1).Defects are defined by their blind frequency f b .The lower the f b value is, the deeper into the material one can probe for defects.

Modal Analysis Theory
Multiple-degree-of-freedom (MDOF) systems can be expressed by the following equation [33]: with M, C, and K are the dynamic mass, damping, and stiffness matrices, and x(t) is the displacement of the MDOF system as a function of time.
Transforming Equation (2) to the Laplace domain yields the transfer function matrix H(s) between displacement and force vectors, X(s) = H(s)F(s): where common-denominator polynomial d(s) = det(Ms 2 + Cs + K), also known as the characteristic polynomial.
When the damping is small, the roots of the characteristic polynomial d(s) are complex conjugate pole pairs, λ m and λ * m , m = 1, ..., N m , where N m is the number of modes in the system.The transfer function H(s) can be rewritten in a pole-residue form, i.e., the so-called "modal" model: with complex conjugate residues R r , R * r , and poles: The pole is a complex number with numerical value of its real part (σ), the damping value, and its imaginary part, the damped natural (modal) frequency (ω).
Figure 7 shows a typical frequency response function measured on the bicycle frame.For a frequency range between 1 and 6000 Hz, 10 modes are selected.The modal parameters are estimated [26].The estimated resonance frequencies of the performed measurements can be found in Table 2.

Finite Element Model Updating
Making use of the modal model described in [14], the frequency response function H ik (ω) between input DOF k and output DOF i at angular frequency ω can be written as where N m is the number of modes, Q r is the modal scaling factor, λ r is the pole number, and r, ψ ir and ψ kr are the mode shapes for mode r and * for the complex conjugate.We will use eigenfrequencies and mode shapes to correlate and update the model parameters to the real, optimized values.
In this section, the automated updating method [22] used to update the finite element frequencies based on two updating parameters will be briefly explained, using poles (resonance frequencies) and corresponding mode shapes as reference data.A non-linear least squares problem is solved in order to minimize an objective function consisting of the sum of the squares of the differences in corresponding resonance frequencies and mode shapes between the test modes and FE modes.In the objective function, the resonance frequencies and mode shapes are both measured but are weighted differently, taking their importance and value into account [34].
The parameter updating problem is stated as minimize (ϑ) where f i and φ ji are the i th natural frequency and associated mode shape at measurement point j, respectively.W f and W φ are the weight coefficients of the frequency and mode shape, respectively; nm denotes the number of modes and np the number of measurement points.Subscripts FE and E denote the items from the FE model and experiments, respectively, ϑ denotes the normalized parameter vector, which allows for variation within the interval [0,1], and (ϑ) is a real valued scalar objective function.

Adaptive Response Surface Method
The optimization routine makes use of adaptive response surface optimization to optimize numerical models based on a multitude of measured data points.The routine is designed to handle multiple-output time series data for the thermal case and frequency data for the modal case [24,27].The optimization procedure can be divided into the following steps: 1.
Starting reference simulation points are selected, and an adapted object function is defined by taking into account the difference between the FE model values and the optimal target values.

2.
The FE model is replaced by an interpolation-based meta-model of response surfaces to decrease optimization time, but this replacement remains an accurate approximation.

3.
The optimization routine is run on the specific object function to find the minimum of the meta-model.An optimization routine is used according to a genetic algorithm computation to find the optimum speed of the complex numerical model.

4.
The estimated parameter values are used as input parameters for an improved FE model that corrects the meta-model.

5.
Only the points that are closer to the minimum are used to form the response surface (pan & zoom function).
The optimized response model is not built from a pre-defined number of design experiments but is adapted and refined during the optimization routine by the pan & zoom command [24].The method automatically calculates the parameter influence sensitivity and sequentially resolves multiple sets with major influencing parameters and fixed parameter values for minor influencing parameters.Afterward, the second parameter set is optimized, and an influence check is then performed for all parameter sets until the results converge.

Thermal Case
The thermal excitation power, excitation duration, and the thermal camera frame rate are used as optimization parameters.The laminate thickness is used as an objection function, and there are physical boundary conditions as well.These limits are shown in Table 3.

Modal Case
The local thickness of each element is optimized using the modal frequencies (global information) estimated from the structural FRFs as well as the mode shapes (local information) measured using the LD-vibrometry measurements.It should be noted that the measurements only apply in one degree of freedom, as the flexure amplitude is only measured perpendicular to the left side of the bicycle frame.

Objective Function
The objective function locates the minimum of the response surface.The minimization is defined in Equation (8), where (ϑ) is defined as the objective function, and T is the response parameter measuring the relative temperature or the resonance frequencies.The response parameter is dependent on the optimized parameter ϑ.The objective function is tested if the boundary conditions are not violated. (ϑ

Results and Discussion
With the numerical model, it is possible to estimate the exposed heating flash, optimally define the heating power for the thermal case, and adjust the mass and structural material properties for the modal case.In this way, the tube thickness of the frame can be measured by the response, considering an anisotropic fiber distribution in the bicycle structure.
The numerical model used in the optimization routine is shown in Figure 8.The model consists of three different material mesh collectors: one for the cover paint with a thickness of 0.11 mm and two for the CFRP plies with a consecutive 90 degree angle between each ply and with the respective global orthotropic ply material properties.In the initial assumption, the full bicycle frame has a continuous thickness of four CFRP plies in a regular stacking of [0/90] s .The material properties and heat exchange properties are defined using the numerical updating process as described and compared with a preliminary static measurement.

Thermal Case
One time instance of the thermal results is shown in Figure 9, where an inter-and intralaminar crack through all plies including the surrounding delaminations is shown.The thickness of this area was difficult to compute due to the crack and type of delamination.The crack showed a thickness (crack depth) of 0 mm, which is below the lower limit.In addition, the delamination was nearly as thick as the minimal thickness.In reality, this thickness is basically the local thickness of the tube, which was further investigated.
The local thickness of different points was examined in comparison to certain benchmark points that could be physically measured.The measurements were performed with an accuracy of 20 µm and are shown in Table 4.It can be seen that the estimated thickness at the three inspected points was accurate, and this was validated by physical measurements.This confirms the two values that cannot be verified without destroying the bicycle and confirms the assumption that a measurement inspection routine can be defined using advanced path planning in combination with thermal FE routines for complex-shaped composite materials.

Model Case
The measured FRF of the used mesh elements of the saddle tube is shown in Figure 7.The global resonance frequencies of the bicycle frame are listed in Table 2.The first natural bending resonance frequency occurred at 102.8 Hz.The excitation shaker was fixed in the middle of the frame at the saddle tube.After the model updating routine was performed, the updated material thicknesses of the evaluated points was changed according to the results in Table 4.
As a result of the updating procedure, an average error of 4.2% in relative difference was achieved.By taking the mode shape in the objective function into account (Equation ( 7)), three modes were switched in the updating process, demonstrating the usefulness of the mode shapes and not just the resonance frequencies in the objective function.If this value were not included, the wrong frequency value differences for these modes would have been minimized.

Comparison
By comparing the updated thickness values of the thermal and structural updating (Table 4), one comes to the conclusions that a thermal model update delivers more precise results than a structural model update.However, some remarks should be made:

•
The mode shapes are only measured in one degree of freedom instead of three, which results in missing mode shape data along the two other degrees of freedom.The main reason for this consists in the complexity of 3D non-contact vibration analysis.However, this problem does not occur when thermal data is used as a reference.• By performing model updating using vibrational data to evaluate local material thickness, it can be seen that the combination of local measurement and global responses biases the local thickness measurement.

•
Although estimating the thermal properties as structural properties is more complex, the thermal results appear to be better matched.

•
The thermal updating routine is far more time-consuming in contrast to modal updating.However, the thermal measurements on the structure are easier to perform.
Modal updating versus thermal updating appears to be comparatively complex when local parameters such as material thickness are evaluated.Considering complex-shaped structures with a high variation in thickness, measurement time during thermal inspection can be feasibly reduced by investing in extensive computation time.This is especially true when a series of similar products is inspected, as the thermal properties are less variable than structural properties during production.The main reason for this is the use of local temperature measurements in model updating.This was demonstrated by the use of local thermal material properties to evaluate local material thickness (Table 4).When the structure contains a wide variety of materials containing high-temperature materials with little conduction, different conclusions may be drawn.Further investigation is needed to make conclusions about those type of structures.

Conclusions
We can assume that the use of correlation and updating routines is of major importance in improving the predictive quality of numerical simulations.It has been demonstrated that techniques similar to numerical model updating can be used for both thermal updating and modal updating.However, an automatic updating routine based on mode shape correlation is thermally not available.The results demonstrate that similar accuracy can be achieved using thermal measurements for updating purposes.In addition, the link between temperature response and material thickness can be approximated using an analytic equation, which facilitates the definition of normalized parameter values in the objective function.

Figure 1 .
Figure 1.Photograph of the thermography measurement setup.

Figure 2 .
Figure 2. Photograph of the modal analysis measurement setup.

• a thin uniform
layer of paint 0.11 mm thick; • a CFRP layer 25% of the overall thickness with a global fiber orientation of 0 • ; • a CFRP layer 25% of the overall thickness with a global fiber orientation of 90 • ; • a CFRP layer 25% of the overall thickness with a global fiber orientation of 90 • ; • a CFRP layer 25% of the overall thickness with a global fiber orientation of 0 • .

Figure 3 .
Figure 3. Overview of the 2D mesh of the bicycle frame.

Figure 4 .
Figure 4. Close-up of the 3D mesh with the different plies schematically represented by different colors.

Figure 7 .
Figure 7. Frequency response function of the expermimental LDV measurement.

Figure 8 .
Figure 8. Numerical finite element model result after 1 s.Temperatures in • C.

Figure 9 .
Figure 9. Thermal response of the crack after PCT (Bin 4).

Table 1 .
Material properties used within the numerical model.