is needed. Color quantification with a smartphone using red green blue (RGB), hue saturation value (HSV), or lightness and chromatic axes A/B (LAB) color spaces have been reported for this purpose and were reviewed recently in the context of smartphone-based biosensors [17], smartphone-based food diagnostics [18], and quantitative LFA [19]. Roda et al. [17] outlines LAB and HSV color spaces as superior to RGB space for measuring small changes in color. However, which of the individual channels should be used or why HSV or LAB color spaces are better than RGB are not mentioned. Moreover, superior performance of the L channel of LAB is equally reported for the smartphone-based detection of hazelnut allergen with LFAs using carbon black based labeling [11] although no comparison with other channels was presented. In another work [20] paper-based detection of the mycotoxin producing black mold (*Stachybotrys chartarum*) is detected by using the R channel. However, no comparison with the B or G channel or other color spaces was reported. In these previous works [18,19] no recommendations on color space use were given although it has been reported that the R channel instead of the combined weighed RGB values lead to background reduction [21]. In another study, the performance of a ΔRGB or ΔLAB system was compared and ΔRGB outperformed ΔLAB for the quantification of plasmonic-ELISA assays [22] and found to perform well in a colorimetric paper-based assay [23]. To obtain Δ values a resultant of R, G, B or L, A, B vectors relative to control values was calculated [22]. However, no comparison between the single R, G, B or L, A, B channels with ΔRGB was performed in either study. Thus, there currently is no consensus regarding which color space/channel should be used for the smartphone image analyses needed for optimal functioning of smartphone hyphenated colorimetric assays such as ELISA, LFA, and other paper-based systems for food contaminant analyses. However, color space/channel choice clearly effects the performance of smartphone-based assays since several of the studies mentioned above reported substantial differences in assay performance in function of the color space/channel used. As a result, conversions to color spaces other than RGB are perhaps executed unnecessarily or suboptimal channels/color spaces may be chosen as a starting point for smartphone-based image analyses based on incomplete recommendations found in the literature. Other issues with current systems are (i) background illumination variation and (ii) inter-phone channel values variation. A light-shielding box is often used to tackle variation in background illumination [9,17–19], however, this detracts from the opportuneness of smartphone-based analyses if add-on items are required. Inter-phone channel value variation is especially important since it is undesirable to develop any system to be used by consumers if bespoke calibration is required for each phone. Unfortunately, most systems described previously were only tested with a single phone [9,11,17–20,24]. Overall, in the scientific literature there appears to be large variations in opinions regarding the use of color spaces and channel combinations. Moreover, no extensive comparison of the individual channels of RGB, LAB, HSV color space, and ΔRGB values with various phones has been identified.

#### *1.3. Workflow Reported in This Study*

In the present study, prediction accuracy of individual channels of RGB, HSV and LAB color space, and ΔRGB was determined for the quantification of color variation, using pH strips, and color intensity variation, using filter paper with dropped nanoparticles often used for LFA (i.e., gold, latex or carbon black nanoparticles) and nanoparticle and oxidized tetramethylbenzidine (TMB) solutions in ELISA wells. Background correction was shown to avoid measurement variation in the absence of direct sunlight more efficiently than a light-shielding box. Inter-phone (n = 6) variation was limited for color change quantification, permitting the quantification of color with phone A using a calibration curve constructed with phones B–F. The optimized system was validated using various phones with model food application exemplars in the quantification of gluten in buffer, bovine milk in goat milk and goat cheese, and pH determination in soil extracts with commercial pH strips and LFA.

#### **2. Materials and Methods**

#### *2.1. Materials*

Dibasic and monobasic sodium phosphate, sodium carbonate/bicarbonate, HCl (37%), 3,3- ,5,5- -tetramethylbenzidine (TMB), NaOH, HAuCl4, sodium citrate, AgNO3, L-ascorbic acid, blue latex beads (LNPs), horseradish peroxidase (HRP), gluten from wheat (crude, ≥ 75% protein), and grade 1 Whatman filter paper were purchased from Sigma-Aldrich (Irvine, UK). Carbon black (N220) (CB) was obtained from Cabot Corporation (Ravenna, Italy), ZEU Proteon Gluten Express ZE/PR/GL25 and ZEU IC-BOVINO lateral flow assays for gluten and cow milk protein detection and ZEU Proteon ELISA were obtained from Zeulab (Zaragoza, Spain). Cow milk and goat cheese were purchased in the local market. Pure goat milk was home produced. DUS alkaline pH strips (5.00–8.50) were purchased from DFI (Gimhae, Korea). Soil was collected in Queen's University garden (Belfast, UK). UV-Vis measurements were performed using a Tecan Safire IIplate reader. Smartphone measurements were performed with a Huawei P8 Lite (12 megapixels (MP)), iPhone 7 (12 MP), Samsung Galaxy Tab E (5 MP), Xiaomi mi5 (16 MP), HTC One M7 (4 MP) and a Samsung Galaxy J7 (13 MP), and a Smartscope Xscience (Ravensburger, The Netherlands) smartphone loop with 3D printed lenses for magnification where mentioned.

#### *2.2. Nanoparticle Synthesis*

Glassware was cleaned with piranha solution and aqua regia to remove all residues. The Turkevich method [25] was used for gold nanoparticle (GNP) synthesis. Briefly, 500 μL of 100 mM HAuCl4 together with 194.5 mL MQ was brought to boil, while stirring, in a round bottom flask equipped with a condenser. At boiling point, 5 ml sodium citrate solution (1% (w/v)), was added and the mixture was left boiling for 30 min then cooled down gradually. Stock GNP concentration (2.8 nM) was estimated using a protocol detailed in [26] and size estimations were reported elsewhere [27].

#### *2.3. pH System (Color Change)*

Citrate-phosphate buffers (0.1 M; pH 5.0, 5.5), phosphate buffer (0.1 M; pH 6.0–7.5), and carbonate/bicarbonate buffers (0.1 M; pH 8.0, 8.5) were used for calibration (pH 5.0, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5) and prediction (pH 5.5, 6.25, 6.75, 7.25, 7.75, 8.25) curves. For buffered soil extracts, 5 g of soil was added to 25 mL of buffer and vortexed. The pH was adjusted to match the buffer pH described above. Samples were allowed to settle 30 min and directly used for pH measurements. Images taken with the Huawei of pH strips dipped in various buffers are shown in Figure 1.

**Figure 1.** Exemplary images taken with the Huawei of pH strips that were used to build calibration curves.

#### *2.4. Nanoparticle Suspensions and Filter Paper Preparation*

A 40X concentration was obtained from a stock GNP solution by centrifuging (13,000 RCF; 30 min; 20 ◦C). Concentrates were used to prepare a 2/3 dilution series from 37.5 up to 2.195 nM to construct calibration curves and a 2/3 dilution series from 31.25 up to 2.74 nM for predictions. For CB, a 10 mg/mL dispersion was sonicated 30 min then diluted to a 2 mg/mL dispersion and sonicated for another 10 min. Serial dilutions (2X) were made until 0.0078125 mg/mL. For predictions, a 1.5 mg/ml CB solution was made and diluted in 2X steps until 0.0117 mg/mL. LNP stock concentration (2.5%) was 2X diluted until

0.00977% to construct calibration curves. For the prediction set a 1.875% solution was diluted in 2X steps until 0.00732422%. Images taken with the Huawei of the nanomaterials on filter paper are shown in Figure 2.

**Figure 2.** Exemplary images taken with the Huawei of the filter papers with nanomaterials on them at various concentrations. Top row are gold nanoparticles (GNP). Middle row are carbon black nanoparticles (CB). Bottom row are latex nanoparticles (LNP). Particle concentration of the solutions used is indicated in the left bottom corner of each image.

#### *2.5. Liquid Assays Preparation*

Colloid GNPs (200 μL) at varying concentration (8 step, 2X serial dilution from 84 nM for calibration and 7 step 2X dilution from 42 nM for prediction) or 150 μL of TMB with HRP at varying concentration (8 step 2X serial dilution from 60 pM for calibration and 6 step 2X dilution from 40 pM for prediction) were pipetted into transparent 96-well plates. TMB enzymatic reaction was stopped after 30 min with 4N H2SO4 (50 μL/well). Absorbance was read at 450 nm for plates with HRP and 513 nm (plasmon peak) for plates with colloid GNPs. Smartphone pictures were taken thereafter. An exemplary image taken with the Xiaomi of the colloid GNPs is shown in Figure 3.

#### *2.6. Sample Preparation and Picture Capturing*

From GNP, CB, and LNP dilutions 5 μL was dropped on filter paper (n = 3), dried, and photographed. pH strips were immersed in buffer, dried, and photographed after 40 s. LFA test strips were photographed 10 min after exposure to extracts. For liquid assays, a phone with white screen was placed under the 96-well plate to provide counter illumination and avoid reflections in the images. All pictures were taken from 5 cm distance with the flashlight on (Figure 4). For background illumination experiments the following light changes were tried: Dark background in a closed windowless laboratory, normal room light (TL light), indirect sunlight in a windowsill, direct sunlight in a windowsill. A black cardboard box (11 × 11 × 5 cm) was made for the Huawei (Figure 4). The box had a hole precisely in the center for the camera and flash. Prediction images taken under varying background illumination were interpolated on calibration curves constructed at room light conditions to test robustness of the background correction applied. All images taken in the box for prediction

were interpolated on a calibration curve constructed with images taken in the box at identical light conditions (room light).

**Figure 3.** Image taken with the Xiaomi of colloid gold nanoparticles in a 96-well plate. Concentrations vary per row (as indicated). Each row (1–4) contains three replicas (H–F) used to construct the calibration curves. Backlight is provided using the white screen of the Huawei placed under the 96- well plate.

**Figure 4.** Left; an image of the Samsung taking an image of a LFA (the ZEU Proteon Gluten Express) under room-light conditions with the flashlight on at 5 cm distance. Right; an image of the light-shielding box used. The dimensions of the box are 11 × 11 × 5 cm to maintain the standard 5 cm as a distance for the photo capturing. The box had a hole in the center precisely for the camera and flash. The phone displayed on the box is the Huawei.

#### *2.7. Commercial Assays*

LFA for gluten, bovine milk, and cheese detection were used according to the instructions of the manufacturer. Briefly, (for estimation of gluten with the commercial LFA assay; ZEU Proteon Gluten Express ZE/PR/GL25), 1 g of gluten was extracted with 10 mL of the given extraction solution, vortexed and centrifuged for 10 min at 3500 RCF. The supernatant was diluted with analysis buffer to obtain 10, 5, 2.5, 1.25, 0.75, 0.5, 0.025, 0.01, and 0.005 ppm of gluten for calibration curve construction and 2, 0.5, 0.3, 0.15, 0.02, and 0.0075 ppm for predictions. The test strip was immersed in 250 μL of these dilutions for 10 min and photographed. For the estimation of cow milk/cheese in goat milk/cheese the ZEU-IC-BOVINO LFA assay was used. Cow milk was spiked into pure, home produced goat milk at 5%, 2.5%, 1.25%, 0.625%, 0.3125%, 0.15625%, 0.078125%, 0.04%, and 0.02% for calibration curve construction and at 3.75%, 1.875%, 0.46875%, 0.234375%, and 0.117188% for predictions. LFAs were immersed in two drops of these extracts diluted with three drops dilution solution. Purity of the goat cheese used was tested using the RC-BOVINO ELISA kit specific for the IgG of cow milk following the manufacturer's instructions. Next extracts of cow cheese and pure goat cheese were prepared by adding 5 g of homogenized cheese to 10 ml H2O and vortexing. Extracts were centrifuged (10 min; 3000 RCF) and goat cheese supernatant was spiked with cow cheese supernatant (final concentrations 8%, 4%, 2%, 1%, 0.75%, 0.5%, and 0.1% of cow cheese) and used for calibration curve generation. Final concentrations of 6%, 3%, 1.5%, and 0.375% cow cheese were used for predictions. LFAs were immersed in one drop of extract diluted with two drops of the given dilution solution and photographed after 10 min. For all commercial assays duplicates were used. Figure 5 shows an image of a set of LFAs used to build a calibration curve for gluten quantification.

**Figure 5.** A set of LFAs used to build a calibration curve for gluten quantification with various smartphones. The concentrations of gluten used are indicated.

#### *2.8. Scoring System*

A scoring system was used to rate the performance of the channels which allowed to generate calibration curves with an R2 > 0.80 (R, G, B, V, L, H, and V). B of LAB (called C of LAC from hereon to avoid confusion with the B channel of RGB) was not considered since it performed only for pH and even there was outperformed by all other channels. R2 values on the regression functions fitted to the predicted concentrations were ranked from highest to lowest. Lowest being 1 point, highest 7 points. The same scoring system was used for the slopes of the regression fits. The channel with a slope closest to 1 got 7 points the next in line 6, etc. If no regression line could be fitted through the predictions the score was 0 points both for prediction and slope. Slope and R2 scores were summed up to get a final score for each channel.

#### *2.9. Software and Data Treatment*

Standard RGB values were detected with a free app from Google Play (RGB Android) (n = 3 per image) and converted to HSV and Cielab with the open source library of colormine.org (last accessed 11 August 2019). S values were rescaled to have a 1–100 range (instead of 0–1) to match RGB scale for visualization. For A and B values 128 points were added to allow a scale from 0–256 and avoid negative numbers complicating the background correction applied.

ΔRGB was calculated as specified in [22]. The applied image background correction was as follows:

Corrected channel values = (Raw signal value)/(Raw background value) × Raw signal value (1)

For paper-based assays white paper was used as the background except for the LFA where control lines were used instead. For liquid assays, wells filled with water adjacent to the tested wells were used as background. All channel values mentioned in the results are corrected channel values. Data analyses was performed with GraphPad Prism 6 Software. For post-hoc analyses, a one-way ANOVA Tukey's multiple comparison test was used. For two-way ANOVA Sidak's multiple comparison test was used. *p*-values were corrected for multiplicity. For calibration curves, a four-parameter dose-response curve was used unless mentioned otherwise. Mean average error percentile (MAE) was calculated on the totality of predictions over the concentration ranges mentioned above. LOD, IC50, and linear range were obtained by interpolating 90%, 50%, and 20%–80% signal values from fitted normalized curves.

#### **3. Results and Discussion**

#### *3.1. Comparing Channel Performance on a Huawei P8*

Various concentrations (n = 3) of GNP, CB, and LNP nanoparticles were dropped on filter paper and photographed. R, G, B, values of all pictures were extracted and used to calculate L, A, B, H, S, V and ΔRGB values. The values were fitted to either a four-parameter dose-response curve or a twophase decay function, which ever resulted in the best fit in terms of R2 value (Table 1). The fitted calibration curves (with an R2 > 0.8) are depicted in Figure 6, left column. Of the LAB color space (called LAC to avoid confusion of B channels with RGB); only L generated adequate calibration curves in all color systems. A curve in all color systems had an R2 < 0.8 and were not used for predictions. C allowed generating a curve for pH determination only and symmetry was observed which made it only possible to predict pH below pH 7. Thus, C was not used for predictions. However, combining A and C might enable color change quantification, as shown previously [28]. From the HSV color space S never reached R2 values above 0.8 and was not used. H was effective in generating calibration curves for pH (as previously shown [29]) and CB, although errors were large in the latter color system with a mean average error (MAE) on all predictions (n = 18) of 90% ± 113%. The V channel worked well, producing calibration curves with R2 > 0.95 in all color systems except for GNP where R2 was 0.88. Prediction images of various pH, GNP, CB, and LNP concentrations were taken and interpolated using the fitted calibration curves of Figure 6. Predicted values are plotted in a scatter plot (Figure 6; right column). The R2 and slope values of the regressions are shown Table 2. For pH predictions, most channels shown in Figure 6 worked for predictions between pH 8.25 and 6.5. At the lowest pH tested (5.5) only R and ΔRGB allowed good predictions. For variation in nanoparticle concentrations low concentrations caused higher variance in all channels. H performed particularly badly in CB predictions. Larger variation in prediction was observed for the color intensity variations as for color change (MAE for the best channels were typically around 25%–35% ± 20%–40% and 1%–2.5% ± 1%–1.5%, respectively). See Section 3.2 for more detail and MAEs for various phones. This might be explained partly by the greater variance in the hand-made replicas for the latter system and partly by the potentially greater effects of background illumination variation on these tests (which was investigated in Section 3.3). To tease apart differences a scoring system (see Section 2.8) was adapted to rank the performance of the channels for each color system individually as well as for all color systems tested (Tables 2 and 3). B, R, and ΔRGB scored the highest overall score, followed by V, L, G, and finally H. H underperformed since that channel was only effective in pH predictions. Moreover, in three out of four color systems either R, B, or ΔRGB had a top score (albeit shared with L in some cases, Table 2). For LNP predictions, the V channel scored highest. However, the R2 values and slopes of V, L, G, and B were close to each other (Table 2) and differences may only reflect random variation. In any case, this comparison shows that conversion to HSV and LAC color space did not increase performance compared to RGB in all color systems tested.

**Figure 6.** The left column shows calibration curves (with an R2 > 0.8) obtained from background corrected channel values from the RGB, HSV, and LAB (Called LAC to avoid confusion of the B channels of RGB and LAB) color spaces for various pH values, GNP, CB, and LNP concentrations. The right column shows scatter plots of predictions obtained with those calibration curves. The yellow line represents a perfect correlation (slope = 1). Color and symbol codes are indicated.

**Table 1.** Function types used to fit data points of the papers with varying GNP, CB, LNP concentration, and pH strips measured in various buffers. R2 values for each fit are mentioned. The fits with R2 > 0.80 were used as calibration curves in Figure 6. C stands for the B channel of LAB which is called LAC throughout the manuscript to avoid confusion of the B channel of LAB with the B channel of RGB. Sigmoidal stands for four-parameter dose-response.

**Table 2.** Slopes and R2 scores of log–log regression functions fitted to predictions of GNP, CB, LNP concentrations, and pH values shown in the scatter plots of Figure 6. A scoring system was used to attribute scores to R2 and slope values for each channel (see Section 2.8). Scores (Table 3) determined ranks for each color system, as well as total scores per channel. Δ is ΔRGB. Ch is channel.

**Table 3.** Scores of the individual channels in each color system, as well as overall scores for all color systems. Shared scores are in bold.

#### *3.2. Comparing Channel Performance between Phones*

The experiments detailed in Section 3.1 were repeated on a tablet and iPhone to test if channel performance behaved similarly over various digital camera devices. Figure 7 shows a comparison of the total mean average percentile error (MAE) for the predictions (pH n = 15; GNP, LNP, and CB n = 18) over various concentrations. MAEs for pH calculations were low and stayed below 2.5% for R and ΔRGB for all phone models. For color intensity, low-end MAEs were typically 20%–30% ± 20%–30% for GNP and LNP predictions. For CB higher MAEs were observed although best functioning channels (B and ΔRGB) showed MAEs below 50%. Two-way ANOVAs (Table 4) show that for pH the variance in MAE caused by channels was significant (*p* = 0.0002), as well as interaction between channels and phone models (*p* < 0.006) although phone models alone did not cause significant variance in MAE (*p* = 0.24). For (LNP), (GNP), and (CB) predictions the phone model caused significant differences (*p* < 0.0001) as well as channel and interaction (*p* < 0.0001) in all but the GNP dataset, where only the phone model effects on MAEs was significant (*p* < 0.0001). Thus, choosing any specific channel as universally ideal for smartphone colorimetric analyses or even for a color specific system seems challenging, especially for color intensity change. Nonetheless, Sidak's post-hoc simple effects within rows of multiple comparisons showed that most of the different effects between phones was limited to one or another channel of the RGB color system, or the H channel (for CB predictions). For ΔRGB no significant effects were observed in this test and thus shows that an error in one of the RGB channels can be compensated in this model, which makes it an interesting option. L and V equally show little variation between phones. However, variance on the MAE as well as absolute errors in the L and V values, although not significant, was larger when compared to the best functioning channel of RGB or ΔRGB in each color system (for L in pH predictions the highest MAE is 3.8% ± 8.4% versus 1.8% ± 1.7% for ΔRGB and for CB predictions the MAE was approximately 40%–50% ± 30% for B and ΔRGB and 50%–80% ± 30%–70% for L and V). Overall, H channel performed poorly in the color intensity experiments and did not outperform R, G, B, V, or L channels for color change predictions. Thus, H was no longer used in the following experiments.

**Figure 7.** Mean average percentile errors (MAEs) for pH (**a**), (LNP) (**b**), (GNP) (**c**), and (CB) (**d**) predictions using an iPhone (dark green), Huawei (dark red), or Tablet (blue). Significant post-hoc Sidak multi comparisons of two-way ANOVAs are indicated. Stars indicate *p*-values with *p*-value correction for multiplicity. \*\* = *p* < 0.01, \*\*\* = *p* < 0.001, \*\*\*\* = *p* < 0.0001.

**Table 4.** Two-way ANOVA on phone model and channel for GNP, CB, LNP, and pH predictions.

#### *3.3. Comparing Box*/*No-Box E*ff*ects on Predictions in Various Background Settings*

The effectiveness of the internal background correction was tested, with a Huawei, by comparing MAEs for pH and [GNP] prediction under varying background illumination with predictions when a light shielding box was used. Channels R, G, B, ΔRGB, V, and L were used since H has proven only functional for pH measurements and even for that application, performance was suboptimal. Little variation in MAEs was observed for pH estimation using the channels R and ΔRGB, even at illumination in direct sunlight (Figure 8). Moreover, two-way ANOVA was only significant for channels and not for interaction or background illumination (Table 5). Interestingly, for R and ΔRGB predictions variation was bigger using the box compared to all background illumination conditions. This may be explained by extensive scatter and unequal light distribution within the box compared to when no box was used. If the background varies between two images but stays equal throughout the individual images then internal background correction should largely correct for it. However, the error can be introduced if background illumination within a picture is patchy. This may also explain the observation by Masawat et al. that a larger box produces less error on predictions as a small one [24]. For [GNP] predictions, no difference was observed between box and no box over the conditions dark, room light, and indirect sunlight for all channels. However, in direct sunlight the background correction applied ceases to function and two-way ANOVA and post-hoc multiple comparisons were highly significant (*p* < 0.0001) for background illumination and direct sunlight, respectively (Table 5). At all other background illumination conditions, the MAE and variation on the MAE were similar in all conditions, including using the box. Thus, the simple background illumination correction applied eliminates the necessity to use a box for all color systems tested if measuring in direct sunlight is avoided.

**Figure 8.** Left, MAEs for pH predictions (n = 15), right, MAEs for [GNP] predictions (n = 18), at various illumination conditions for R, G, B, ΔRGB, V, and L channels without use of a box and in a box at room illumination conditions. Images under dark (black balls) background were taken in a dark room; room background (red balls) in a windowless laboratory illuminated with a tube light bulb; indirect sunlight (green balls) in a windowsill at a cloudy day; direct sunlight (dark yellow balls) in a windowsill in full sunlight.

**Table 5.** Two-way ANOVA on background illumination and channel for GNP, and pH predictions.

#### *3.4. Channel Performance for Colored Liquids*

Channel performance to quantify color intensity in liquid solutions (as typically done for ELISA) was investigated using various concentrations of colloid GNPs (mimicking plasmonic ELISA) and various amounts of HRP to oxidize TMB (mimicking standard HRP-based ELISA). Color change was measured using a benchtop spectrometer and a smartphone (Xiaomi) using R, G, B, ΔRGB, V, and L channels. All four-parameter dose response curves showed good fits (R2 > 0.8) and were used to predict GNP and HRP concentrations. Predictions and linear regressions are shown in Figure 9a,b. Slopes and R2 values of the calibration curves and R2 and slope values of the linear regression functions are shown in Table 6. For [HRP] prediction, the spectrometer, B, G, and ΔRGB channels were acceptable. For [GNP] prediction the spectrometer and G, B, ΔRGB, and L channels were acceptable (Figure 9a; Table 6). For comparison, the calibration curves of these channels and the spectrometer were normalized (Figure 9c,d). LOD, linear range, and IC50 for both assays are shown in Table 6. Interestingly, G, B, and ΔRGB channels for colloid [GNP] had a ~3X lower LOD as the absorption curve of the spectrometer (~1.5 versus ~7.5 nM, respectively). For the [HRP] curves, B and ΔRGB channels gave lower LODs as the spectrometer (~2.3 versus ~3.4 pM, respectively). Thus, the smartphone-based system was slightly more sensitive than the spectrometer. However, the linear range was reduced compared to the spectrometer. MAEs were calculated for G, B, ΔRGB, and spectrometer [HRP] predictions. For G and B only three concentrations (n = 3 × 3) were used since the linear range was reduced. For ΔRGB and spectrometer four concentrations (n = 4 × 3) were used. One-way ANOVA showed that MAEs for [HRP] prediction using G and B channels did not significantly differ from the MAEs for [HRP] prediction when the spectrometer was used (although with reduced range). ΔRGB had slightly higher MAEs (*p* < 0.05 Tukey post-hoc; Figure 9c inset). This is probably due to the error introduced into ΔRGB from the R channel which varied significantly and could not be used to build a calibration curve. For [GNP] predictions ΔRGB and L channels had significantly higher MAEs as G and B and the spectrometer (*p* < 0.01; Tukey post-hoc; Figure 9d inset). Again, the R channel could not be used to build a calibration curve.

#### *3.5. Channel Performance Comparison Using Commercial LFA for Milk Allergen Detection*

Commercial LFA for the quantification of bovine milk in goat milk and cow cheese in goat cheese (Figure 10a,b) were used to further challenge the smartphone-based quantification using the Huawei. For the LFA used to quantify bovine milk in goat milk (Figure 10a) G, B, L, and V channels showed promising calibration curves (R2 > 0.9). However, the curve determined with the V channel values was flat (min–max difference about 25 corrected value units). Equally, the L channel was quite flat and had slightly lower R2 values (0.91) as B and G (0.93 and 0.94, respectively). ΔRGB did not allow construction of a calibration curve, which was most likely caused by variations observed in the R channel values. Thus, only B, G, and L curves were used for predictions. Linear regressions on predictions were good for each channel (R2 > 0.95). However, the L channel did not predict at 0.23% cow milk (Figure 11). Goat cheese extract was spiked with cow cheese and used to test the applicability of the assay for the identification of cow protein in goat cheese. Here, only the B channel gave a satisfactory calibration curve (R2 = 0.951; Figure 10b). Predictions showed excellent linear regression (R2 > 0.98) (Figure 10b, inset).

**Figure 9.** (**a**) Linear regression functions for oxidation of tetramethylbenzidine (TMB) with horseradish peroxidase (HRP) in ELISA wells for the smartphone analyses using various channels and a benchtop ELISA plate spectrometer. (**b**) Linear regression functions for colloid GNP at various concentrations in ELISA wells for the smartphone analyses using various channels and a benchtop ELISA plate spectrometer. Color codes and channels used are indicated. (**c**,**d**) normalized four-parameter dose-response curve fits for HRP oxidation of TMB (*a*) and colloid GNP (*b*) in ELISA wells. Green and blue circles stand for G and B channels, grey diamonds ΔRGB, magenta triangles absorption values measured by a spectrometer. (**c**,**d** inset) MAEs calculated from the predictions shown in (**a**,**b**). Stars indicate *p*-values from post-hoc analyses with *p*-value correction for multiplicity. \* = *p* < 0.05, \*\* = *p* < 0.01.

**Table 6.** Analytical parameters of prediction and calibration curves for GNP and TMB solution ELISA. R2 and slopes given for prediction curves are from linear regression lines. All calibration curves were

**Figure 10.** (**a**) G (green balls), B (blue balls), L (black triangles), and V (turquoise squares) channel values were fitted to a calibration curve for LFA quantification of cow milk spiked into pure goat milk. ΔRGB values (grey diamonds) did not allow to fit to a curve. (**b**) B channel values fitted to a calibration curve and, (inset), linear regression on predictions for cow milk in goat milk at three concentrations (n = 6).

**Figure 11.** Predictions of cow milk (*x*-axis) in goat milk using L, G, and B channels. Each replica (n = 2) is indicated. Dashed line represents a curve at 45 degrees.

#### *3.6. Channel Performance Comparison Using Commercial LFA for Gluten and pH Strips for Soil pH Prediction*

Next, the ability to use various phones to quantify gluten with a commercial LFA, as well as pH strip quantification in buffered soil was tested (Figure 12). B values performed again optimally for the LFA quantification and were used to construct calibration curves for the quantification of gluten in a buffer using four different phones (Xiaomi, Huawei, iPhone, and Tablet) (Figure 12a). All phones enabled the construction of calibration curves (R2 = 0.93; 0.88; 0.94; 0.86 for Xiaomi, Huawei, iPhone, and Tablet, respectively). Unfortunately, the curves did not overlap, except for the Tablet and Huawei, and cannot be directly used by end-users without calibration of individual phones. Magnification of LFA strips was attempted to correct inter-phone variation. LFA were photographed under the lens of a low-cost instrument (the smartscope; Ravensburger). Here, calibration curves with G channel values showed the highest R2 values (0.87, 0.95, 0.93, and 0.92 for Xiaomi, Huawei, iPhone, and Tablet, respectively) (Figure 12b). Although the distance between the curves was narrower as in the previous set-up, the overlap of all curves was insufficient for the use of various phones to predict values on the same calibration curve. Thus, for color intensity variation using LFA, it seems that two- point calibration or camera calibration is necessary. Some suggested methods exist to obtain such calibration and improve sensitivity [30–33]. One interesting concept is to use a color reference chart in order to stabilize color variations caused by built-in automatic image correction operations [30]. Another option would be adjusting the white balance of the phones to a standardized value. This may be done by using a reference grey card or by locking the exposure and gain while selecting a preset white balance

in a control window that is illuminated with an external phone-independent constant light source [32]. Setting the white balance, gain, and exposure can be done in various manual camera applications (e.g., Open Camera, ProShot) and using such a system in combination with a greycard or external constant light source (flashlights of various phones can be different in light output) may decrease inter-phone variation. Moreover, higher sensitivity may be reached by adjusting the exposure time of the assay in such manual camera applications as was previously shown for other smartphone-based LFA quantification for the detection of bacterial fruit blotch [33]. However, the four-parameter dose-response curves did overlap for all four phones for pH determination of buffered soil samples using ΔRGB (Figure 12c) and R (Figure 12d) using the Xiaomi, Huawei, iPhone, and Tablet (R2 > 0.99) showing that the proposed method worked without additional adjustments for this assay. Overlap was slightly less when using ΔRGB values, thus, the R channel was used from hereon for this application. Another two phone models (HTC and Samsung) were included to check the universality of the system. Four-parameter dose-response curves showed again excellent fits (R2 > 0.99) and overlap with the other curves (Figure 12d). Finally, MAEs of predictions (n = 15) were calculated for each phone using its own calibration curve or a calibration curve constructed with all phone models except the phone model used to take the images for the predictions (Figure 12d; inset). Two-way ANOVA analyses of this data was significant for the type of curve (*p* < 0.001). MAEs increased significantly when predictions were performed on phone A with calibration curves prepared using phones B–F. However, for all phones except the tablet this increase in MAE was below 3% which explains why post-hoc simple effects within phone models and multiple comparisons was only significant for the tablet (*p* < 0.001). Thus, direct quantification of color change by end-users with the pH strips without the need of phone calibration seems possible if slightly increased MAEs are acceptable.

**Figure 12.** (**a**) B channel values fitted to four calibration curves for LFA quantification of gluten spiked into buffer with four phones (models indicated). (**b**) Calibration curves for LFA quantification using the smartscope and corrected G values. (**c**) Calibration curves for pH in buffered soil extract using corrected ΔRGB values. (**d**) R channel values fitted to calibration curves for pH estimation of buffered soil extracts with six phones (models indicated). Inset, the results of post-hoc analyses on the MAEs for pH predictions (n = 15) using the phone's specific calibration curve or a mixed curve of all phone models except the one used for prediction. Stars indicate *p*-values with *p*-value correction for multiplicity. \*\*\* = *p* < 0.001. All calibration curves were made using a four-parameter dose-response function.

#### **4. Conclusions**

The ability of all channels of the RGB, HSV, and LAB (named LAC here) color spaces to quantify color intensity and color change variations was tested in paper based and liquid assays. Channel performance was compared under varying background illumination using various phones. A channel of LAC proved not suitable for color quantification and C channel performance was suboptimal. L channel showed good performance in most systems tested (as did V) but both never outperformed the best RGB channel for a specific test. S performed suboptimal in all systems tested. H channel performed satisfactorily for color change but not for color intensity quantification. Moreover, H channel never outperformed the best RGB channel for a specific test. Overall, R functioned best for color change while B and G worked best for color intensity variation. ΔRGB values did show some robustness towards errors in an individual channel of RGB thus providing more universal applicability. However, ΔRGB stops to perform if an error in an individual channel of the RGB color space is too high. Thus, all RGB channels and ΔRGB should be initially plotted for assay development to determine the optimum channel. However, conversions to LAC or HSV color space are unnecessary. These results are specific for colorimetric assays and the implications might not hold for fluorescence or chemiluminescence-based assays where using luminance values seems a more logical choice. This being said, for such assays the light measured will equally pass the Bayer filter of the phone and is converted to RGB values. Mathematically converting these RGB values to LAB values may lead to similar results as reported here. However, a detailed characterization of the performance of RGB, LAB, and HSV color spaces/channels of such luminance-based assays with various phones and fluorophores is needed to test this assumption. Using a light-shielding box to prevent the error caused by background variation was less efficient than internal background correction if images were not taken in direct sunlight. Thus, use of a box is superfluous for colorimetric analysis. For color intensity variation it was shown that images taken with phone A could not be used on a calibration curve taken with phone B. Magnification did improve this situation but did not completely resolve the problem. Thus, camera calibration, white balance adjustment, and exposure time adjustment should be considered for LFA quantification. For color change quantification (pH determination of buffered soil extracts) calibration curves of the six phone models tested overlapped significantly. Thus, color change quantification by end-users without a light-shielding box or phone specific calibration seems feasible. Future research will focus on reducing prediction error and inter-phone variation for color change and intensity-based assays by combining channels from RGB, HSV, and CieLAB color spaces with machine learning algorithms.

**Author Contributions:** J.L.D.N. was the main author executing the experiments, data curation, and writing. L.B. assisted in the execution of experiments and data curation. Y.Z. and K.M.B. assisted in executing the experiments. C.T.E. and K.R. provided advice and assisted in reviewing and editing the final document. K.C. contributed to the conceptualization, supervised this work, and provided assistance in reviewing and editing the manuscript.

**Funding:** This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No 720325.

**Acknowledgments:** We further want to thank Master Sparrow Nelis for providing his smartscope for image magnification and Ilyana Nelis for providing home produced pure goat milk.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **An Innovative Ultrasonic Apparatus and Technology for Diagnosis of Freeze-Drying Process**

#### **Chin-Chi Cheng 1,\*, Yen-Hsiang Tseng 2 and Shih-Chang Huang 1**

Received: 17 April 2019; Accepted: 8 May 2019; Published: 11 May 2019

**Abstract:** The freeze-drying process removes water from a product through freezing, sublimation and desorption procedures. However, the extreme conditions of the freeze-drying environment, such as the limited space, vacuum and freezing temperatures of as much as −50 ◦C, may block the ability to use certain diagnostic sensors. In this paper, an ultrasonic transducer (UT) is integrated onto the bottom of a specially designed frozen bottle for the purpose of observing the freeze-drying process of water at varying amounts. The temperatures and visual observations made with a camera are then compared with the corresponding ultrasonic signatures. Among all of the diagnostic tools and technologies available, only ultrasonic and visual records are able to analyze the entire progression of the freeze-drying process of water. Compared with typical experiment settings, the indication of drying point for water by the amplitude variations of ultrasonic L3 echo could reduce the process period and energy consumption. This study demonstrates how an innovative frozen bottle, an integrated ultrasonic sensor and diagnostic methods used to measure and optimize the freeze-drying process of water can save energy.

**Keywords:** freeze-drying process diagnosis; ultrasonic transducer (UT); freezing/drying point; drying period

#### **1. Introduction**

Dehydration extends food's usage period longer than that of fresh food by preserving it in a stable and safe condition [1,2]. The conventional methods of drying include solar drying, air drying, spray drying [3], microwave drying [4], infrared drying, fluidized bed drying [5], spouted bed drying, vacuum drying and freeze-drying [6]. Drying methods can be separated into natural and artificial categories. Artificial methods are more advantageous than the natural methods [7] because they can remove large amounts of moisture efficiently by being able to control the different parameters involved such as the temperature, drying air flux and time of drying and so forth. [8]. For long-term storage of food, drug and biopharmaceutical products, most manufacturers utilize the freeze drying process, due to its advantages including better stability, easy handling and storage, as well as better overall product quality [9,10].

The freeze-drying process is comprised of the freezing, primary drying and then secondary drying stages [11] to remove the water from a product. The water contained in food is cooled down and becomes ice during the freezing stage. This stage governs the sublimation and desorption rates and the quality of the lyophilized product [12]. Then, in the primary drying stage, the air in the vacuum chamber is exhausted and the chamber pressure is reduced below the vapor pressure of ice. Meanwhile, the shelf temperature increases gradually to sublimate the ice. The residual water inside the food will be desorbed thoroughly in the secondary drying stage [13].

In the freezing stage, the accuracy of the freezing point data is related to the water activity, frozen water [14], freezing and the thawing of frozen food. The freezing point is also important in order to estimate the freezing time, the end point of freezing and the fraction of unfrozen water in food. Then, in the primary drying stage, the chamber pressure and food temperature [15] decides the sublimation rate of the ice, primary drying period, dried pore structure and product quality. A lower product temperature and the corresponding lower vapor pressure of the ice can result in longer primary drying times and higher manufacturing costs [16]. However, if the product temperature increases above the "critical formulation temperature," this may lead to losing the pore structure, to shrinkage or to fully collapsing. In the secondary drying stage, the product temperature controls the rate of desorption and the obtainable moisture level. Further reduction of the chamber pressure is not necessary in this stage [17]. Secondary drying times are usually designed to achieve a reduction of moisture content within the dried product to less than 1% [18]. Due to the importance of these parameters and its wider applications to the freeze drying process, the estimated freezing/drying points and the ability to model these properties is crucial in food processing (freezing and drying) and food stability during storage [10,19].

Being able to perform real-time process diagnosis will be beneficial to understanding the complex interplay between the different elements during the freeze-drying process to better enhance quality control procedures. However, the extreme conditions of the freeze-drying environment, such as the limited space, vacuum and freezing temperatures as low as −50 ◦C, may block the ability to use certain diagnostic sensors [6,20]. The typically used diagnostic tools for real-time measurement of the freeze-drying process are temperature and pressure sensors [20,21]. In order to provide temperature information, a temperature sensor is put in direct contact with a frozen sample [21]. However, inconveniences arise because it is often not easy to remove the sensor after the freeze-drying process. Additionally, using this diagnostic method on drug and biopharmaceutical products may cause undesired contamination due to a sensor being in direct contact with the sample. The chamber pressure is measured by the pressure sensor or vacuum gauge. The measured chamber pressure value cannot directly indicate the physical variation of frozen product real-time during the freeze-drying process [20,22]. In this paper, an idea for diagnosing the physical property variations of a product in the frozen bottle to provide additional information for a machine system is presented, as shown in Figure 1.

**Figure 1.** Controlling structure of the freeze-drying process by using the innovative frozen bottle method.

In order to circumvent the limitations of the mentioned sensors and provide real-time process diagnosis, the ultrasonic technique is one of the most widely known non-destructive and non-intrusive methods [23]. The fundamental properties of ultrasonic signals include reflection and transmission coefficients, velocity, attenuation and scatter signals given off from the materials. These signatures have specific relationships with properties of materials [24], process variations [25] and sample quality [26]. Ultrasonic echoes can also detect the temperature [27,28] and pressure [29] of a material. Passot et al. utilized ultrasonic technology to decrease the sublimation time during the primary drying stage by

controlling the nucleation temperature [30]. The authors presented an innovatively designed frozen bottle with an ultrasonic transducer and ultrasonic technology to be utilized for the diagnosis of the freezing point during the frozen process in our previous paper [22]. In this study, the ultrasonic technique will be applied for real-time diagnosis of the freeze-drying process. An ultrasonic transducer (UT) is integrated into the bottom of a frozen bottle. During the freeze-drying process at various amounts of water, the process is analyzed by using ultrasonic technology to evaluate the freezing/drying points and the drying period to optimize the process and save energy.

#### **2. Design of Frozen Bottle with Ultrasonic Transducer**

In order to solve the aforementioned limitations regarding the transmission of signals in a vacuum, as well as the limitations of being unable to view the process using the steel freezing bottle apparatus or the opaque samples, an innovative freezing bottle apparatus was designed as shown in Figure 2. The innovatively designed container used for freezing is comprised of two parts: a steel plate as the bottom holder and a transparent polymethylmethacrylate (PMMA) tube as the side wall of the container. Figure 2a shows the side view of the frozen bottle. A steel plate with a diameter of 25 mm and a height of 6mm was designed for better heat transmission with the heating/cooling shelf. A transparent PMMA tube measuring 31 mm high, with a diameter of 21 mm and a thickness of 2 mm was designed for better visibility. The steel plate and PMMA tube were then glued together. Figure 2b shows the top view of the apparatus. Figure 2c shows the bottom view of the apparatus. The UT is comprised of Pb(ZrxTi1−x)O3 (PZT) material and integrated into the cavity of steel plate using the sol-gel spray technique, as described in previous publications [22,24,31]. The operational temperature range of the UT is between −100 and 400 ◦C.

**Figure 2.** Photographs of the designed apparatus from (**a**) the side, (**b**) the top and (**c**) the bottom view.

To understand the transmission paths of the ultrasonic signals within the apparatus, a schematic view of the apparatus placed upright is displayed in Figure 3. A temperature sensor (Type T thermocouple) was installed in the middle of the apparatus for measuring the water/ice temperature. The temperature can be recorded and then compared with the ultrasonic signals received during the freeze-drying process. As shown in Figure 3a, an ultrasonic transducer (UT) was integrated onto the bottom of the steel plate. When electric pulses are applied to the piezoelectric film, the ultrasonic signals are transmitted into the steel plate. Ln (n = 1, 2, ... ) denotes the nth round trip of longitudinal ultrasonic echoes reflected from the interface of the steel plate/water or ice and Lw is the echo reflected from the water/air interface. The L3 and Lw echoes are used to monitor the freeze-drying process and water/ice state. The height of the water level and the thickness of the bottom of the apparatus are denoted as H and D, respectively.

**Figure 3.** (**a**) Schematic drawing of cross-section of the designed apparatus with an ultrasonic transducer (UT) and thermocouple and (**b**) shows the typical ultrasonic echoes acquired by the UT.

When the water was fed into the apparatus, the typical ultrasonic signals acquired by the UT are presented in Figure 3b. It was observed that the L1 echo, which reflected from the steel plate/water or ice interface, appeared as 0.81μs and remained so during the entire process. The operating frequency of the ultrasonic transducer was 9.13 MHz and the −6 dB bandwidth was 8.0 MHz. The signal to noise ratio (SNR) for the ultrasonic L1 echo was 38.5 dB. The Lw echo, which propagated in the water and reflected from the water/air interface, was observed at 20.39 μs. Following the L1 and Lw echoes, there were several echoes reflected from the steel plate/water or ice interface. The desired ultrasonic echoes were detected by the acquiring windows of the acquiring system. The time delay difference between the L1 and Lw echoes is denoted as Δt.

#### **3. Experiments**

In order to confirm the hypothesis that the apparatus was able to perform improved diagnostics of the freeze-drying process in a vacuum, a 4 L air-cooled split-type shelf freeze-drying machine, equipped with vacuum chamber, shelf, control units, refrigeration system and vacuum pump, was used to carry out all experiments, as shown in Figure 4a. The vacuum chamber was used for freezing and drying the samples under temperature conditions ranging from +50 to −40 ◦C and pressures from 760 to 0.05 Torr. The apparatus was set on the shelf of vacuum chamber for the freeze-drying process. The control unit was composed of a programmable logic controller (PLC) and a human-machine interface (HMI) for controlling the freeze-drying process. The measured temperature data from the thermocouple were recorded by the PLC control unit every second during the freeze-drying process.

The ultrasonic signal is triggered and received by the pulse/receive unit (5072PR, Olympus, Japanese), as shown in Figure 4b. It is equipped with a broadband negative spike pulser and receiver, which can be operated in reflection or transmission mode. All the experiments in this study were carried out using the ultrasonic pulse/receiver unit's pulse-echo mode. The ultrasonic signals were acquired every second during the entire cycle. The ultrasonic signals received during the freeze-drying process were monitored and recorded with the digital storage USB oscilloscope (DSO-U2400, Perytech, Taipei, Taiwan), as shown in Figure 4c.

To compare the temperatures and ultrasonic signals, a charge-coupled device (CCD) camera, digital video recorder (DVR) and PC were used to capture dynamic images of the apparatus during the freeze-drying process, as shown in Figure 4d–e. The CCD camera (CV-M10BX, JAI, Japan) is a progressive scan camera with a standard interlaced video output at a resolution compatible with VGA or SVGA formats. The camera was equipped with a 35 mm lens for better resolution. The DVR (08KD, Kingnet, New Taipei City, Taiwan) is an 8-channel recorder that uses the H.264 image compression format.

**Figure 4.** Photographs of (**a**) the split-type freeze-drying machine with vacuum chamber, controller unit, refrigeration unit and vacuum pump, (**b**) the ultrasonic pulse/receive unit, (**c**) a PC-based digital oscilloscope, (**d**) Charge-coupled device (CCD) camera, (**e**) digital video recorder (DVR) and PC.

In this study, different amounts of water at levels of 5, 10 and 15 mm were put into the apparatus for the duration of the freeze-drying process. The settings of a typical freeze-drying procedure are outlined in Table 1. The shelf temperature was set at −30 ◦C to freeze the water and the freezing period was 120 min. During the freezing period, the water in the apparatus/container was cooled down and it began to freeze. After enough cooling and freezing, the chamber was exhausted and the chamber pressure was reduced from 760 to 0.17 Torr within several seconds. Then the primary and secondary drying processes were started and the ice sublimated into gas continuously. The shelf temperature was incrementally increased by 10 degree intervals from −30 to 10 ◦C at heating and pausing periods of 60 and 120 min, respectively. The pausing period of 0 ◦C is 300 min to sublimate ice further. During the sublimation of the ice, the chamber pressure was reduced further to 0.036 Torr.

**Table 1.** The typical procedure used for the freeze-drying process.

#### **4. Results and Discussions**

For exploring the physical phenomena that occur during the freeze-drying process and searching for the indicators of freezing or drying completeness, the measured temperature and pressure in the vacuum chamber, ultrasonic echoes and velocity of the water/ice and visual observation during the freeze-drying process are analyzed in the following sections.

#### *4.1. Temperature and Pressure Variations in the Vacuum Chamber during the Freeze-Drying Process*

Temperature and pressure variations during the freeze-drying process have a close relationship with the state of the water/ice and both affect the size of the ice crystal nuclei and the freeze-dried product's quality. In the experiment, the measured room, shelf, cold trap, water/ice temperatures and pressures for the freezing and drying periods of the freeze-drying procedure of the water at a level of 15 mm are shown in Figure 5a,b, respectively. In Figure 5a, the room temperature was 32 ◦C during the entire process and the chamber pressure was kept at 760 Torr for the experiment's duration of 120 min. The shelf and cold trap temperatures were set at −30 and −40 ◦C and reached the desired values at 40 and 15 min, respectively. The water temperature was 32.1 ◦C in the beginning and was cooled down to a supercooled state of −0.4 ◦C, 16.1 min into the experiment. After the latent heat of the water was released completely, the phase change from water to ice was completed at 43.7 min and the ice temperature was −16.5 ◦C. Therefore, the cooling (ΔPCT) and freezing (ΔPFT) period determined by the sample temperature could be defined as the period from the start of cooling to the supercooled state and from the supercooled state to the phase change end, respectively.

**Figure 5.** Temperature (room, shelf, sample, cold trap) and pressure (chamber) variations during freeze-drying process at (**a**) the freezing stage and (**b**) drying stage.

At the 120-min mark of the experiment, the drying procedure started and the chamber pressure was reduced from 760 to 0.17 Torr. The shelf temperature was increased in 10 ◦C increments, from −30 to 10 ◦C, according to the designed schedule in Table 1. The cold trap temperature was cooled down further to condense the vapor exhausted from the chamber. The chamber pressure was reduced further to 0.036 Torr during the drying process. At the 370-min mark, the ice temperature increased suddenly from −36 to −10 ◦C due to the sublimation of the ice and gradually increased to 23 ◦C by the end of the experiment. This abnormal situation was due to the detachment of the thermocouple from the ice and will be confirmed in the following section. Consequently, the ice temperatures during the second drying period is not available. The drying end could not be estimated from the temperature information and the typical experiment procedure would end at the 1020-min mark of the experiment.

#### *4.2. Amplitude Variation of Ultrasonic Signals during the Freeze-Drying Process*

The application of freeze-drying technology can guarantee the product's quality. However, the process diagnosis and control will benefit by guaranteeing product quality and energy conservation. The conventional tools for monitoring this process may be limited to visual observation and temperature measurement. According to the authors' knowledge, utilizing ultrasonic technology to monitor the freeze-drying process is rare. In order to investigate the correlation between the observed ultrasonic signals and the freeze-drying process of water, the amplitudes of the L3 and Lw echoes in Figure 3b with respect to the process time were acquired and presented in Figure 6. The details of freeze-drying process measured by the amplitude variations of ultrasonic L3 and Lw echoes were demonstrated as follows:

**Figure 6.** Amplitude variations of ultrasonic L3 and Lw echoes during freeze-drying process at (**a**) freezing stage and (**b**) drying stage, indicating the water-in, freezing point (ice crystal), phase change end and drying end.

#### *4.3. Variation of Ultrasonic Velocity during the Freeze-Drying Process*

During the freezing process, the slow cooling speed caused bigger ice crystal nuclei than did the fast cooling time. Ultrasonic velocity has a close relationship with temperature and it is one of the indicators that shows the cooling speed of water, the phase change end and the start of sublimation. The ultrasonic velocity in the water is calculated according to the following equation:

$$\mathbf{v}\_{\rm W} = 2\mathbf{H}/\Delta\mathbf{t} \tag{1}$$

where H is the height of the water level in Figure 3a and Δt is the time delay between the ultrasonic L1 and Lw echoes in Figure 3b. The result is shown in Figure 7.

In Figure 7, when the water was poured into the container, the ultrasonic velocity was first measured at the 3.0-min mark under the water temperature of 27.2 ◦C. At that moment, the ultrasonic velocity in the water was 1588.5 m/s. From the 3.0 to the 17.4-min mark, when the water temperature decreased to 1.0 ◦C, the ultrasonic velocity decreased to 1434.5 m/s, indicating that the temperature of the water had decreased during the cooling process. Comparing Figure 7 with Figure 5a, the relationship between ultrasonic velocity and water temperature is expressed by the following equation.

$$\text{tr}\_{\text{W}} \left( \text{T}\_{\text{W}} \right) \;= 1.42418 \times 10^3 + 1.03190 \times 10 \times T\_w - 0.15453 \times T\_w^2 \tag{2}$$

where vw (m/s) is the ultrasonic velocity in water and Tw is the water temperature. The ultrasonic velocity decreased with the reducing water temperature. The ultrasonic velocity presented by authors in this study is higher than that presented by Bilanuik and Wang [34], due to some minerals contained in the water. Therefore, the cooling speed of water can be determined by measuring the ultrasonic velocity in the water. From the 17.4 to 38.3-min mark of the experiment, the ultrasonic velocity disappeared and the amplitude of the Lw echo also decreased in noise level due to the formation of ice crystal nuclei. After the 38.3-min mark, the ultrasonic Lw echo appeared again and the ultrasonic velocity gradually increased from 3586.4 to 3597.2 m/s, due to the phase change end of the water and the completion of the freezing process. From 43.0 to 120.0-min mark, the ultrasonic velocity maintained a steady level of 3597.2 m/s for the pausing period. Therefore, ultrasonic velocity can be used to indicate the cooling speed, the phase change end and completion of the freezing process of water.

**Figure 7.** Ultrasonic velocity variation in the water during the freeze-drying process, which indicates the water-in, the freezing point (ice crystallization) and the phase change end.

#### *4.4. Visual Observation of the Container during the Freeze-Drying Process*

The dynamic physical phenomena that can be observed in the container shows the phase change from water to ice and the sublimation detail of the ice during the freeze-drying process in great detail. These phenomena could also be evidence of other measured parameters. Figure 8a–l are the photographs of the container with water/ice during the freeze-drying process at the 0–888-min mark, respectively. The freeze-drying process observed in Figure 8a–l is described as follows:

**Figure 8.** Photographs of the water/ice in the container during the freeze-drying process along with the corresponding process times of (**a**) 0 min, (**b**–**d**) 15 min 54~56 s, (**e**) 22 min, (**f**) 40 min, (**g**) 240 min, (**h**) 300 min, (**i**) 480 min, (**j**) 540 min, (**k**) 840 min and (**l**) 888 min.

The freeze-drying process can be observed clearly in Figure 8. These results can provide additional evidence for the previous measuring methods. However, this observing method is only suitable for the transparent container and plain solution. For the unobservable container and the opaque material, such as metal container and colored samples, the temperature and ultrasonic techniques are more capable to diagnose the mentioned phenomena. From the measured results, the thermocouple detached from the ice and further temperature information of the ice was not available. In this viewpoint, the ultrasonic technique is a more reliable and capable method for diagnosing the freeze-drying process.

#### *4.5. Freeze-Drying Processes of Various Water Levels*

Before applying the ultrasonic technology to the freeze-drying process, the linear relationship between ultrasonic signatures and the water level is the fundamental requirement. To verify this relationship, the containers, which were filled with water in the height of 5, 10, 15 mm, passed through the freeze-drying process under the shelf temperature setting of −30 ◦C. In this research, the variation during the drying period will be evaluated. During the drying stage, the thermocouple cannot measure the samples' temperature profiles due to the detachment of thermocouple from the sample. The ultrasonic signals are able to diagnose the drying phenomena. The experimental results of the amplitude of the ultrasonic L3 echo are shown in Figure 9. Figure 9a–c are the amplitude variations of ultrasonic L3 echo in water levels of 5, 10, 15 mm, respectively, during the drying stage. The timings of the drying end indicated by ultrasonic L3 echo are 329.8, 526.7 and 885.6-min mark for the water levels of 5, 10 and 15 mm, respectively. The drying periods of various water levels are also indicated. The corresponding timings of the drying end indicated by ultrasonic L3 echo are illustrated in Table 2. It seems that there is a linear relationship between the timings of drying end and the water level.

**Figure 9.** Amplitude variations of ultrasonic L3 echo with respect to various water levels at drying stage. Water level: (**a**) 5, (**b**) 10 and (**c**) 15 mm.

**Table 2.** Timings and temperatures of water in, freezing point, phase change end and drying end indicated by temperature and ultrasonic L3 echo for water levels of 5, 10, 15 mm.

The measured drying periods are compared with the water level for evaluating the mentioned linear relationship. The results are shown in Figure 10. The estimated error of drying period under the experimental conditions is less than 5%. In the water level range from 5 to 15 mm, the average drying periods indicated by ultrasonic L3 echo increase from 185.1 to 742.8 min. The slope of the fitting line is 55.8 min/mm. The drying period is expressed as:

$$
\Delta P\_{DlIT} = -112.8 + 55.8 \,\text{\*} \, H \tag{3}
$$

where ΔPDUT is the drying period in Figure 9 and H is the water level in Figure 3. Even though there are no temperature and visual information, one can still estimate the required drying period according to the filled water level based on Equation (3). In the future, this ultrasonic sensor and evaluating technology would be installed into freeze drying machine to detect the phase change and feedback the information to machine control system, as shown in Figure 1. Comparing with the typical experiment settings of 1020 min, the ultrasonic technique can clearly indicate the dynamic phenomena and completion of cooling/freezing/drying stages at each water level before the end of the experiment for reducing the freeze-drying process and saving energy.

**Figure 10.** Drying periods measured by ultrasonic signal with respect to various water levels. Water level: 5, 10, 15 mm.

#### **5. Conclusions**

In this paper, an ultrasonic transducer (UT) is integrated onto the bottom of a specially designed container to analyze the freeze-drying process of water at various water amounts. The measured ultrasonic signatures are compared with the temperature and visual records. Among these three measured methods, ultrasonic and visual records are able to document the entire progression of the freeze-drying process, including the water-in, freezing/drying points, the phase change end of the water and the cooling/freezing/drying periods. The ultrasonic velocity in the water also indicates the cooling tendency of the water. The drying period increases with the water level. The increase rate, which is evaluated by the amplitude of the ultrasonic L3 echo, is 55.8 min/mm. During the drying stage, the thermocouple cannot measure the entire temperature profile of the ice due to the detachment of the thermocouple from the ice. Only the ultrasonic signals and visual records are utilized for diagnosing the drying phenomena. However, the observing method is only suitable for the transparent container and plain solution. Comparing with the typical experiment settings, the indication of drying ends for water at level of 15 mm by the amplitude variations of ultrasonic L3 echo could reduce the process period of 134.4 min/cycle and save 13% of consumed electricity. Therefore, this study demonstrates the use of a specially designed container, integrated ultrasonic sensor and technology for analyzing and optimizing the freeze-drying process of water for saving the process cost and energy.

**Author Contributions:** C.-C. Cheng initiated the idea, provided the draft and completed the revising of the paper. Y.-H. Tseng prepared the integrated sensor and practiced the experiments of freeze-drying process. S.-C. Huang analyzed the data and assisted to revise paper.

**Funding:** This work is funded by Ministry of Science and Technology of Taiwan, ROC under Contracts No: MOST 104-2221-E-027-121, MOST-105-2221-E-027-103 and MOST 107-3113-E-008-003.

**Acknowledgments:** The authors would like to provide special thanks to Hong-Ping Cheng (Department of Energy and Refrigerating Air-Conditioning Engineering, National Taipei University of Technology) for providing the useful discussions, Rei-Zhu Hsueh (Tai Yiaeh Enterprise Co., Ltd., New Taipei City, Taiwan) for providing useful maintaining assistance of freeze drying machine.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Temperature and Strain Correlation of Bridge Parallel Structure Based on Vibrating Wire Strain Sensor**

### **Lu Peng 1,\*, Genqiang Jing 1, Zhu Luo 1, Xin Yuan 2, Yixu Wang 1 and Bing Zhang 1**

Received: 29 October 2019; Accepted: 21 January 2020; Published: 24 January 2020

**Abstract:** Deformation is a ubiquitous phenomenon in nature. This process usually refers to the change in shape, size, and position of an object in the time and spatial domain under various loads. Under normal circumstances, during engineering construction, technicians are generally required to monitor the safe operation of structural facilities in the transportation field and the health of bridge, because monitoring in the engineering process plays an important role in construction safety. Considering the reliability risk of sensors after a long-time work period, such as signal drift, accurate measurement of strain gauges is inseparable from the value traceability system of high-precision strain gauges. In this study, two vibrating wire strain gauges with the same working principle were measured using the parallel method at similar positions. First, based on the principle of time series, the experiment used high-frequency dynamic acquisition to measure the thermometer strain of two vibrating wire strain gauges. Second, this experiment analyzed the correlation between strain and temperature measured separately. Under the condition of different prestress, this experiment studied the influencing relationship of temperature corresponding variable. In this experiment, the measurement repetitiveness was analyzed using the meteorology knowledge of single sensor data, focused on researching the influence of temperature and prestress effect on sensors by analyzing differences of their measurement results in a specified situation. Then, the reliability and stability of dynamic vibrating wire strain gauge were verified in the experiment. The final conclusion of the experiment is the actual engineering in the later stage. Onsite online meteorology in the application provides support.

**Keywords:** parallel position; bridge structure; temperature; vibrating wire strain sensor

#### **1. Introduction**

Bridge health monitoring and diagnostic discriminant models have always been a key challenge for the transportation sector worldwide. In a previous study, on a regular and irregular basis, construction workers used different monitoring instruments to test some components of structure, analyzed the data, and finally evaluated the performance of structure. In 2018, Mao by monitoring dynamic characteristics of Sutong cable-stayed bridge (SCB), including acceleration and strain responses, as well as modal frequencies, are investigated through one-year continuous monitoring data under operating conditions by the structural health monitoring system. One-year continuous modal frequencies of SCB are identified using the Hilbert–Huang transform method. Variability analysis of the structural modal frequencies due to environmental temperature and operational traffic is then conducted. Results show that temperature is the most important environmental factor for vertical and torsional modal frequencies. The traffic load is the second critical factor especially for the fundamental vertical frequency of SCB [1]. In 2019, Wang reported that the real-time monitoring data collected from a long-span cable-stayed bridge is utilized to demonstrate the feasibility of the improved BDLM-based method. In particular, the present BDLM-based method allows for probabilistic forecasts, offering substantial information about the target TIS response, such as mean and confidence interval. Results show that the improved BDLM is capable of capturing the relationship between temperature and TIS. Compared to the AR model, multiple linear regression (MLR) model, and BDLM without the AR component, the improved BDLM shows better forecasting performance in modeling and forecasting the TIS of a long-span bridge [2]. Such measurement methods have the drawback of large errors and discontinuities. The ideal health monitoring system should accurately reflect the change in grassroot structure under the influence of factors such as temperature, humidity, and other environmental factors, installation and deployment methods, and the sensor's own error, and establish an effective health assessment and prediction model. For example, Mao proposed according to one-year continuous monitoring of strain data recorded by the structural health monitoring system of the Sutong cable-stayed bridge, the lifetime fatigue reliability of three welded details of the orthotropic steel deck was investigated, detailed analysis of the separated components of the raw strain data was first conducted, and included the slow-varying trend and the dynamic component. The strain dynamic component was mainly induced by the local vehicle axle loads. Rainflow counting was used to obtain the stress range histograms, which were then used to calculate the equivalent stress range according to the lognormal-fitting method. Finally, a time-dependent fatigue reliability evaluation of the described three welded details was conducted using one-year monitoring strain data. Results showed that the fatigue performance of the two welded details, RTDD and DTD, remained satisfactory after 100 years of operation because the failure probabilities were both lower than 10−5. The designed cutout of the diaphragm was applied to the RTD weld at the welded connection between the U rib and diaphragm. This cutout was validated as a means to help achieve better fatigue resistance for the RTD weld [3]. For bridge strain monitoring sensors, most experimental procedures use resistive strain gauges, vibrating wire strain gauges, fiber grating strain gauges, etc. Although the development of fiber grating technology in recent years has led to a new round of equipment upgrades for structural monitoring methods, vibrating wire strain gauges are most widely used in stage bridge health monitoring systems.

In recent years, studies on the heat of vibrating wire sensors at home and abroad have continued unabated, and domestic and foreign experts and scholars have conducted a lot of research to improve the sensor performance. The research on sensors is gradually extended towards high precision, large scale, small volume, and multiple applications. During the improvement in the performance of the sensor, higher requirements are also imposed on the measurement accuracy. To improve accuracy, in 2010, He proposed a low-voltage excitation with the feedback method. By pre-excitations, the vibration frequency of the sensor can be used as the output of the driver. The feedback signal frequency is very close to the real frequency of the sensor, so the wire can reach resonance state quickly. The optimal excitation strategy was verified by the new designed detection circuit. The data perform as short-time excitations with large resonance amplitude, therefore the anti-interference ability got enhanced at lower cost of signal processing circuit, finally increasing the measurement time and improving the accuracy of measured frequency [4]. In 2010, Wen proposed a method for the frequency measurement of vibrating wire sensors with LM3S6965 as the control core. An equal precision measurement was used to effectively improve the measurement accuracy of the system [5]. In 2016, Tian et al. proposed a design method of nonlinear compensation that is provided for the nonlinear relation between the measured force and the frequency of vibrating wire sensor. Data density is increased by interpolation based on the principle of cubic spline interpolation. The F-f curve of vibrating wire sensor is fitted and revised by MATLAB based on the least-square method, more accurately the smooth curve according to engineering practice is obtained and a reasonable conclusion is obtained. Experiments show that the method achieves the accurate fitting to the F-f curve of vibrating wire sensor with the limited data. The F-f function of vibrating wire sensor which is long buried can be rapidly determined by this procedure [6]. In 2017, Chen et al. proposed an

adaptive vibration measurement method based on fast Fourier transform, which uses digital Fourier filtering to automatically filter noise interference and then uses Quinn algorithm for high-precision frequency calculation. Based on the STM32 processor platform, the method was used in vibrating a wire sensor frequency measurement system. The test results show that the relative error of system frequency measurement is less than 0.01% in the case of no noise interference; in the case of severe white noise interference (signal to noise ratio is −20 dB), the system's frequency measurement relative error is less than 0.3%. Compared to other frequency measurement methods, it was found that the method has better noise suppression and frequency measurement accuracy [7]. In terms of static metrology and calibration, the experiment proved that a vibrating wire strain gauge can be widely used in the field of geotechnical engineering, and the accuracy requirement for variable measurement is relatively low. However, in the application of traffic engineering, the measurement accuracy is correspondingly improved, accompanied by the needs of metrology and calibration. Many scholars used a strain sensor calibration device consisting of a calibration frame and digital display to calibrate certain metering characteristics of strain gauge. For example, in 2009, Xu invented a vibrating wire strain sensor calibration device that uses a dial gauge as a gauge for deformation, with a deformation range of up to 3 mm and a strain measurement resolution of about 0.1%. The results show that the strain in the middle of the sensitive grid of strain gauge is the largest and gradually decreases to zero at both ends when measuring. The force has no effect on the sensitive gate base layer and bonding layer at both ends and there is a strain transition zone between the base layer and the sensitive gate, the bonding layer, and the base layer. The longer the horizontal width of bonding layer is the thinner the thickness and the larger the shear modulus will be, and the strain transfer will be more efficient if the active zone of the sensitive grid is longer [8]. In 2011, Zhang et al. designed a vibrating wire strain sensor calibration device. They used a grating scale as the measuring standard for deformation. The maximum calibration distance was 300 mm, and the measurement resolution was 2 μm. The degree is about 0.2% [9]. In 2016, Mai et al. invented a vibrating wire strain gauge calibration device. A dial gauge was used as the measurement standard for deformation amount, so that the resolution was 0.01 mm [10]. During engineering application, many in-depth studies have been conducted on the dynamic calibration of strain gauges. In 2016, Bai et al. affixed FBG sensors on standard beams. Loads were applied to the standard beams to deform them, and the beams were measured. The surface strain was used to calibrate the sensitivity coefficient of FBG strain gauge, and the measured strain range could be analyzed up to 2000 με [11]. In 2017, Zhang proposed a method for the indirect calibration of fiber grating strain sensors using a fiber grating temperature sensor. The lossless calibration of strain sensors was achieved [12]. During the detection, monitoring, and metering of vibrating wire strain gauges, the temperature deformation caused by temperature changes has been the research direction of many scientists. Chen et al. monitored data in engineering applications based on the working principle of vibrating wire strain gauges. The relationship between temperature and strain was analyzed using the relevant data of bridge strain monitoring, and the relationship between temperature influence and strain of the string itself and the temperature field of structural section was evaluated [13]. Bai et al. analyzed that when the temperature changes significantly, a variety of mathematical models were used to fit and calibrate the monitoring data, eliminating the temperature drift of strain gauge, and the experimental results reflect the original characteristics of deformation [14]. Agostiono investigation of the sensing features of the long-period fiber gratings (LPGs) fabricated in hollow core photonic crystal fibers (HC-PCFs) by the pressure assisted electric arc discharge (EAD) technique. In particular, the characterization of the LPG in terms of shift in resonant wavelengths and changes in attenuation band depth to the environmental parameters: Strain, temperature, curvature, refractive index, and pressure is presented. Results show that LPGs in HC-PCFs represent a novel high performance sensing platform for measurements of different physical parameters including strain, temperature and, especially, for measurements of environmental pressure. The pressure sensitivity enhancement is about four times greater if comparing LPGs in HC and standard fibers. Moreover, differently from LPGs in standard fibers, these LPGs realized in innovative fibers, i.e., the HC-PCFs, are not sensitive to the

surrounding refractive index. During the online calibration study of strain sensors [15], offline removal and reinstallation of sensing elements pose a risk to the continuity and consistency of monitoring data. At present, the load test method is mainly used, but the load test method cannot eliminate the instability of structure such as the stiffness and strength of the structure itself [16,17]. Therefore, a more effective method at this stage is to install a traceable high-frequency dynamic monitoring sensor in parallel to the sensor to perform parallel measurement, thus achieving the online calibration of a long-term monitoring sensor. At present stage, the measurement calibration of the strain sensor system in various industries only adopts static calibration before installation or carries out a static and dynamic load test for calibration, and does not carry out online calibration during use. The static calibration cannot determine the influence of the error and noise, temperature change, and prestress effect on the strain value. Static and dynamic load tests cost a lot of manpower and materials, and cannot provide reference for sensor calibration under complicated passive excitation. In the online calibration research of strain sensors, because the offline disassembly and reinstallation of the sensing elements pose a risk to the continuity and consistency of the monitoring data, the load test method is mainly used at this stage, but the load test cannot eliminate the rigidity and strength of the structure itself.

In the research process, we carried out a number of relevant tests. For example, the calibration of fixed excitation by a simply supported beam model is shown in Figure 1. The strain correlation test of a small simply supported beam at constant temperature in a high and low temperature box is shown in Figure 2. Fatigue test verification is shown in Figure 3. In order to further verify the test, a test was carried out on Jiujiang Bridge to verify Figure 4. Most of the above studies adopted common source excitation schemes, mainly to verify the calibration of the strain monitoring system under common source excitation by parallel measurement methods, which provided important verification support for the coupling relationship between temperature, prestress, and strain output in the online calibration of dynamic strain measurement in this study.

**Figure 1.** Schematic diagram of simply supported beam test model.

**Figure 2.** High and low temperature dynamic strain test. (**A**)Test equipment, (**B**) Assembly drawing of experiment beam.

**Figure 3.** Verified by fatigue tests on bridge model.

**Figure 4.** Test verification of Jiujiang Bridge. (**A**) Jiujiang Bridge, (**B**) Test process.

During this study, two vibrating wire strain gauges with the same working principle were measured using the parallel method at similar positions. The feasibility of a parallel measurement scheme was verified [18]. Under the condition of different prestresses, this experiment studied the influencing relationship of temperature corresponding variable [19,20]. We focused on researching the influence of temperature and prestress effect on sensors by analyzing differences of their measurement results in a specified situation. The measurement repetitiveness was analyzed using the meteorology knowledge of single sensor data, and then the reliability and stability of a dynamic vibrating wire strain gauge were verified in the experiment. The relevant measurement performance of a dynamic vibrating wire strain gauge was demonstrated.

#### **2. Working Principle of Vibrating Wire Strain Gauge**

A vibrating wire sensor is tested by steel string vibration. An experiment is conducted to characterize the force according to the variation in vibration frequency. In an actual output frequency signal, there is no strain gauge; it must be field calibration, signal drift, long distance transmission, and long time. The problem of poor durability was used, and the robustness was good [21]. This experiment solves the shortcoming of unstable strain gauges for a long time, and this product can be widely used in bridge monitoring.

A vibrating wire sensor has good measurement characteristics; it can achieve nonlinear characteristics of less than 0.1%, sensitivity of 0.05%, and less than 0.1%/10 ◦C temperature error.

After the strain gauge is manufactured, its steel string has a certain initial tensile force *T*0 and thus has an initial frequency *F*0. When the strain gauge is installed, the tensile force of vibrating wire changes with deformation, and the strain can be measured by the tensile force change of vibrating wire. Set the tension of the vibrating wire to *T* and the natural frequency to *f*. The relationship between tension and frequency can be expressed as Equation (1):

$$T = \mathbb{K}f^2\tag{1}$$

where K is related to the length of string, and the mass per unit length can be expressed using Equation (2):

$$
\Delta T = T - T\_0 = K(f^2 - f\_0^2) \tag{2}
$$

Assuming that the strain increment of strain gauge can be set to the strain increment of vibrating wire, Equation (3) can be derived as follows:

$$
\varepsilon\_h = \varepsilon\_\mathcal{S} = \frac{\Delta K}{EA} \tag{3}
$$

When *EA* is the axial stiffness of steel string, it can be derived from Equation (4) as follows:

$$
\varepsilon\_{\rm li} = \frac{K}{EA} (f^2 - f\_0^2) = k\_{\rm li} (f^2 - f\_0^2) \tag{4}
$$

A mathematical model of the vibrating wire sensor can be expressed using Equations (5) and (6) as follows:

$$\mathbf{F} = \mathbb{K}(f^2 - f\_0^2) \tag{5}$$

$$\mathbf{F} = A(f^2 - f\_0^2) + B(f - f\_0) \tag{6}$$

When the length is such that the fine string of mass m is subjected to tension F (Figure 5), the natural frequency *f* can be expressed as Equations (7)–(9) as follows:

$$f = \frac{1}{2} \sqrt{\frac{F}{ml}}\tag{7}$$

$$f = \frac{1}{2l} \sqrt{\frac{ES\Delta l}{\rho l}} = \frac{1}{2l} \sqrt{\frac{E\Delta l}{\rho\_v l}}\tag{8}$$

$$f = \varphi(F) \tag{9}$$

**Figure 5.** Working principle diagram of vibrating wire strain gauge.

Sensitivity can be derived using Equation (10) as follows:

$$f^2 = \frac{1}{4l^2} \frac{E\Delta l}{\rho\_v l} = \text{K}\varepsilon \tag{10}$$

After differentiation, Equations (11) and (12) can be deduced as follows:

$$2fdf = \mathbb{K}d\varepsilon\tag{11}$$

$$\mathbf{k} = \frac{df}{d\varepsilon} = \frac{\mathbf{K}}{2f} \tag{12}$$

The material coefficient can be calculated using Equation (13):

$$\begin{cases} K = \frac{1}{4l^2} \frac{E}{\rho\_v} = \frac{ES}{4l^2 \rho} \\ \varepsilon = \frac{\Delta l}{l} \end{cases} \tag{13}$$

The above formula shows that sensitivity *k* is directly proportional to material coefficient *K* and inversely proportional to the vibration frequency of string.

After many experiments, the initial frequency is *f*0 when the measured tension is *F*0, and the vibration frequency is *f*1 when the measured tension is *F*1 = *f*0 + *f*. The nonlinear error of vibrating wire strain gauge can be expressed using Equations (14) and (15) as follows:

$$f\_1 = \frac{1}{2} \sqrt{\frac{F\_0 + \Delta F}{ml}} = \frac{1}{2} \sqrt{\frac{F\_0}{ml}} \sqrt{1 + \frac{\Delta F}{F\_0}} = f\_0 \sqrt{1 + \frac{\Delta F}{F\_0}} = f\_0 \sqrt{1 + \varepsilon\_F} = f\_0 \left(1 + \varepsilon\_F\right)^{\frac{1}{2}}\tag{14}$$

$$f\_1 = f\_0(1 + \frac{1}{2}\varepsilon\_F - \frac{1}{8}\varepsilon\_F^2 + \frac{1}{16}\varepsilon\_F^3 - \cdots) \tag{15}$$

At that time *F*2 = *F*0 − Δ*F*, leading to Equation (16):

$$f\_2 = f\_0(1 + \frac{1}{2}\varepsilon\_F - \frac{1}{8}\varepsilon\_F^2 + \frac{1}{16}\varepsilon\_F^3 - \cdots) \tag{16}$$

Its quadratic nonlinearity error can be expressed as Equation (17):

$$\frac{\left|\frac{1}{8}f\_0\varepsilon\_F^2\right|}{\frac{1}{2}f\_0\varepsilon\_F} = \frac{1}{4}\varepsilon\_F\tag{17}$$

The above formula shows that the larger ε*F*, the larger δ*m*. At the same time, the ambient temperature mainly influences the frequency stability. The bulk density ρ*v* and Δ*l* caused by *F* do not change with the ambient temperature, so the frequency stability can be expressed using Equation (18):

$$
\gamma\_f = \frac{df}{f} = \frac{dE}{2}E - \frac{3}{2}\frac{dl}{l} \tag{18}
$$

When the ambient temperature changes, the vibrating wire strain gauge material and the measured structural material have different linear expansion coefficients, and the sensor is subjected to additional stretching or compression. The additional strain can be expressed as Equation (19):

$$
\varepsilon \mathbf{T} = (\alpha - \beta) \bullet \Delta \mathbf{T} \tag{19}
$$

where:

ε*T*: Additional strain caused by temperature effect;

α: Linear expansion coefficient (◦C−1) of the structural material to be tested;

B: The coefficient of linear expansion (◦C−1) of a steel string of a vibrating wire strain gauge;

ΔT: Temperature change amount.

In actual application, the vibrating wire strain sensing system generally adopts software compensation. After the thermistor is set to collect the working environment temperature in the electromagnetic coil, the optimized temperature compensation algorithm is combined with the software to compensate in the demodulation instrument. The temperature strain compensation of the vibrating wire strain gauge mostly utilizes the two-dimensional regression method, polynomial fitting method, and neural network method. After the test results of the vibrating wire strain gauge are linearly fitted, the strain calculation method can be expressed as Equation (20):

$$\boldsymbol{\varepsilon} = \mathbf{a} \times \left( f\_i^4 - f\_0^4 \right) + \mathbf{b} \left( f\_i^2 - f\_0^2 \right) + \mathbf{k}\_T \left( T\_i - T\_0 \right) \tag{20}$$

where:

ε: The dependent variable of the current time relative to the initial position (10−6);

*k*: The steel string strain gauge minimum reading 10−6/(kHz2);

*f* 2 *i* : The Steel string strain gauge current output modulus kHz2;

*f* 2 0 : The Steel string strain gauge initial output frequency modulus kHz2;

*kT*: The steel string strain gauge temperature correction factor 10−6 ◦C;

*Ti*: The steel string strain gauge current time temperature value (◦C);

*T*0: Temperature value when measuring *f*0 ( ◦C).

When using polynomial fitting, the coefficients a and b were calculated using the least squares fitting method, and the strain calculation method for integrating temperature changes is shown in Equation (21):

$$\boldsymbol{\varepsilon} = \mathbf{a} \times (f\_i^4 - f\_0^4) + \mathbf{b}(f\_i^2 - f\_0^2) + \mathbf{k}\_T (T\_i - T\_0) \tag{21}$$

The abovementioned various theoretical calculation methods and vibrating wire strain gauge measurement characteristics can better measure the strain, eliminate the strain generated by environmental influence, more accurately understand the mechanical strain of structure, and analyze the stress state of structural facility.

#### **3. Test Plan**

Two sensors were installed in parallel on the 45th steel tooling for free acquisition at the same frequency. The device (Figure 6) shows the material properties of No. 45 steel (Table 1). The acquisition device is equipped with a high-frequency dynamic acquisition device, which is a dynamic vibrating wire automatic acquisition system. This system uses a nonsweeping technology scheme to prevent the steel string vibration from attenuating, and an embedded mirror oscillation circuit to ensure excitation frequency. By matching the phase and true motion of steel cord and simultaneously detecting the resonant frequency of steel string through several waveform periods, noise immunity and resolution of measurement have largely improved compared with the flat-domain periodic averaging method.

**Figure 6.** Caption for figure schematic diagram of the test fixture.

**Table 1.** Forty-five steel material characteristics.

The measurement module used the patented 8,671,758 products, and the dynamic measurement rate is 20 to 333 Hz. At the same time of dynamic measurement, the module also performs auxiliary measurement, performs static measurement at 1 Hz, provides finer measurement resolution, and better anti-interference to external noise source performance. The thermistor input signal of each vibrating string acquisition channel was measured at a high resolution of 24 bits at 1 Hz. The performance of thermistor parameters is shown in Table 2. The excitation module has a resolution of 26 mV and a dynamic measurement rate of 20, 50, 100, 200, and 333.33 bHz. The range of sensor resonance frequency is shown in Table 3. The measurement frequency accuracy is ± (0.005% reading + measurement resolution). The noise level corresponding to the measurement resolution b at different sampling rates is shown in Tables 2–4.

**Table 2.** Thermistor performance parameters.

**Table 4.** Measurement resolution b (typical value for a 2.5 kHz resonant frequency sensor).

In the measurement, both steel plate and strain gauge will produce different deformations due to temperature changes. To minimize the acquisition error, a correlation curve between temperature and strain was analyzed. During the test, welding was used. To solder the two sensor holders on the tooling and for better verification, the effect of temperature response, and correlation, when collecting the zero point, one sensor was set to the free state, one sensor passed the fixture. The force was applied such that the sensor zero point acquisition and the first sensor's zero point acquisition strain difference is about 140 με.

#### **4. Analysis of Test Data**

The relationship between temperature and strain measurement shows that the strain changes correspondingly when the temperature changes. The strain reversely sways as the temperature changes. We got a good corresponding relationship. The temperature and strain data collected using the two sensors A and B were compared, and the correlation between the two data was analyzed. During the analysis, there are usually two ways of correlation: Autocorrelation and cross-correlation. The autocorrelation function is known as the autocorrelation equation and used to describe the correlation between the correlation functions of related data at different times as shown in Equation (22):

$$R\_f(\tau) = f(\tau) \* f'(-\tau) = \int\_{-\infty}^{\infty} f(t + \tau)f'(t)dt = \int\_{-\infty}^{\infty} f(t)f'(t - \tau)dt\tag{22}$$

At the same time, cross-correlation or cross-covariance were also used to represent a measure of similarity between two signals. Cross-correlation mainly analyzes the degree of correlation between two time series. Cross-correlation is essentially similar to the convolution of two functions. For discrete functions *fi* and *gi*, it can be defined as Equation (23):

$$(f\*g)\_i = \sum\_j f\_j^\* g\_{i+j} \tag{23}$$

If the continuous signal is set to two sets of *f*(*x*) and *g*(*x*), then the cross-correlation is defined as Equation (24):

$$(f\*g)(x) = \int f^\*(t)g(x+t)dt\tag{24}$$

During this test analysis, the data collected using the two sensors A and B were compared, and the correlation between the two sets of data was analyzed. Their similarities were derived to verify the accuracy and reliability of dynamic acquisition of data obtained using sensors A and B. The correlation coefficient between the temperature data of sensors A and B was calculated to be 0.9983, and the correlation degree is highly correlated. The correlation coefficient between the strain data of sensors A and B is 0.9895, and the correlation degree is highly correlated. Cross-calculation can be obtained. The correlation between temperature and strain of sensor A is −0.6683, and the correlation between temperature and strain of sensor B is −0.5573. The experimental results show that the values are significantly correlated (Figure 7).

(**C**) Temperature change of 1.8 °C

**Figure 7.** Free state strain values of strain sensors A and B in a temperature change of 0.2, 0.6, and 1.8 ◦C.

During the test, to better use the knowledge of metrology to analyze the dynamic RMS change output, verify the reliability of linear relationship between temperature and strain monitoring, and repeat the calculation of temperature measurement and strain measurement, the experiment measured the comprehensive reflection of various random influencing factors, including the instability of tooling materials used, random error caused by the laying process, instability of sensing instrument, environmental conditions, and other factors, as well as the actual measured randomness. The measured object also affects the dispersion of measured values, especially when the random variation in the measured object during the bridge monitoring process is large. Therefore, the dispersion of experimentally measured values is typically slightly greater than the dispersion introduced by the sensor static calibration standard itself. To be less affected by outliers, all the deformations and the corresponding temperature were measured. During the calculation of repeatability of the experiment, the arithmetic mean value of measured value should be calculated first, and the calculation formula can be expressed as Equation (25):

$$\overline{F} = \frac{\sum\_{i=1}^{n} F\_i}{n} \tag{25}$$

where:

*Fi*—The strain measurement indication of the *i*-th measurement, με;

*n*—The number of measurements.

The measurement repeatability can be quantitatively expressed using the experimental standard deviation Sr(y), and the calculation formula can be expressed as Equation (26):

$$\text{Sr (y)} = \sqrt{\frac{\sum\_{i=1}^{n} \left(F\_i - \overline{F}\right)^2}{n-1}} \tag{26}$$

where:

*Fi*—The measured strain measured at the *i*-th measurement, με; *F*—The arithmetic mean of indications of strain measurement, με; *n*—number of measurements.

Using the above formula, under the synchronous acquisition test conditions, sensors A and B were repetitively used for the temperature and strain measurements, respectively. The repeatability of sensor A strain indication is 0.16071, and the temperature measurement repeatability is 0.09971. The repeatability of strain measurement of sensor B is 0.11743, and the repeatability of temperature measurement is 0.08209. The repeatability shows that the stability of sensor B is better, and the experiment can be performed under the condition that the initial state of sensor is small. The measured values are more stable and more clearly characterized by the corresponding relationship between temperature and strain.

At the same time, the relationship between strain and temperature was analyzed during a temperature variation of 0.2, 0.6, and 1.8 ◦C. Figures 4–6 show that the correlation did not change with the change in temperature. A stable negative correlation curve relationship was maintained. Among them, sensor A has a relatively weak induction of temperature and deformation under the prestrain of 140 με, which can be characterized as shown in the figure. In the absence of external excitation, the change in experimental numerical temperature is less than the perception of sensor B. At the same time, the corresponding strain produces a small change in the output. Sensor A temperature and the strain relationship map shows that in the free state of acquisition device, the experimental temperature will also have a corresponding change during small changes. However, both sensors A and B can consistently respond to the corresponding output temperature and strain signals (Figure 8).

(**A**) Strain-temperature relationship at 0.2 °C

(**B**) Strain-temperature relationship at 0.6 °C.

**Figure 8.** Strain-temperature relationship between 0.2, 0.6, and 1.8 ◦C intervals for strain gauges A and B.

After calculating different temperature changes, the experimental values show the correlation between strain and temperature. By analyzing the strain-temperature correlation function, it can be found that the temperature change is less than that of sensor B when there is no external excitation. At the same time, the change in output corresponding to the strain is small. In different temperature variation ranges, the experimental data were used to analyze the relationship between the temperature and strain of sensor A. It was observed that the acquisition device is in a free state, and the temperature is small. The effect of corresponding change will be produced. According to statistical analysis, the strain-temperature correlation can be expressed as Δμε = Δ*T* + 510, and the strain measurement result has a negative correlation with the temperature change. From the relationship between the temperature and corresponding strain curves of sensor B in the free state, it can be observed that the strain-temperature correlation can be expressed as Δμε = Δ*T* in the state where sensor B is more free under the same conditions as sensor A. Δμε = Δ*T* + 367; the strain measurement results have a stable negative correlation with temperature change. During data acquisition, it was found that strain gauge A is larger at the same temperature due to prestressing, the strain output is larger than the output of sensor B, and the correlation between temperature, peak, and valley is more prominent. It was assumed

that 0.2 is a one-fold change in temperature. During the three-fold change and nine-fold change, the steady-state strain curve is stable, and the correlation function is not fluctuating. However, in the state where the prestrain of sensor A is 140 με, the induction of temperature and deformation is relatively weak, and the difference is about 3–6 με, which is about 2–5% of the prestress, as shown in Figures 7 and 8. In the later test, more theoretical and experimental verifications are needed. Combined with the artificial intelligence model training method, the expansion effect caused by the material itself after the vibration of vibrating string was analyzed, and the environmental impact error was maximized. The tooling surface was optimized for temperature conductivity. Since this test is the calibration test for the maximum fitting, instrumentation error, etc., the welding process was used, so that the above effects can be neglected. At the same time, referring to Table 4, the noise of the device is 9.775 × 10−9 με, which can be ignored. Through comprehensive evaluation, it was found that the method of synchronously collecting sensors at similar positions can achieve the online comparison of the corresponding variable measuring sensors. The later test can be calibrated by high-frequency dynamic acquisition to calibrate the low-frequency acquisition and combined with the relevant optimization algorithm for calculation. The maximum optimization restores the original deformation magnitude, and at the same time reduces the instability of constitutive performance parameters such as the stiffness and strength of structure. The experimental results show the calibration of the corresponding variable monitoring sensor.

#### **5. Conclusions**

In engineering, the strain test method is usually used to determine the actual stress state of structure and monitor the structural deformation increment, and the structural strain test is an important method to solve the structural strength problem. During strain detection and monitoring, the calibration of strain gauge plays a significant role in the accuracy of data. To ensure the accuracy of traceability of magnitude, we should study the corresponding relationship between the temperature and strain of strain gauge and the online calibration test method, which is important for the online metrological traceability system of strain monitoring sensor. We focused on researching the influence of temperature and prestress effect on sensors by analyzing differences of their measurement results in specified situation. In this study, a dynamic vibrating wire measurement module was used to dynamically collect the temperature and strain of two vibrating wire sensors in the free state, and the correlation between temperature and strain was analyzed. The correlation coefficient between the temperature data of sensors A and B were statistically analyzed. The degree of relevance 0.9983 is highly correlated. The experimental results show that the correlation coefficient between the strain data of sensors A and B is 0.9895, and the correlation degree is highly correlated. After cross-calculation, a significant correlation was observed between temperature and strain indications. An experiment was carried out to analyze the measurement repeatability by using the knowledge of metrology. After analysis, it was found that the repeatability of temperature and strain values is within 0.2, and the consistency between the measured results of the same measurement is good. At the same time, the relationship between strain and temperature in different temperature variation intervals was analyzed. It was found that the temperature change has a small effect on the strain in the prestressed state, but the relationship between temperature and strain is relatively stable. Finally, the experiment verified the feasibility of parallel measurement scheme corresponding to the online calibration of variable monitoring sensor and analyzed the effect of temperature corresponding to the variable measurement under different prestress conditions, providing theoretical and experimental verification support for establishing an online metrological verification standard. This study can also better confirm the accuracy of bridge health monitoring data, effectively reduce the monitoring error caused by temperature drift, effectively eliminate the impact of environmental load, and better feedback the essential deformation of bridge structure. Considering the reliability risk of sensors after a long-time work period, such as signal drift, the approach proposed in this paper can be used as a substitute to such strain monitoring system in order to determine whether these sensors should be replaced. We expect our approach to be expandable

to various applications. It can provide important support for the effective monitoring data extraction of the bridge health monitoring system.

**Author Contributions:** Conceptualization, L.P., G.J., Z.L., Y.W., and B.Z.; Data curation, L.P., X.Y., and Y.W.; Formal analysis, Y.W.; Funding acquisition, L.P.; Writing of original draft, L.P.; Writing—review and editing, L.P. All authors have read and agree to the published version of the manuscript.

**Funding:** This research was funded by the Project of National Key R & D Program of China grant number 2017YFF0206305 and The Key Projects of Public Scientific Research Institutes (2018-9031).

**Acknowledgments:** The authors would like to acknowledge the support of National Key R & D Program of China under grant number 2017YFF0206305 and the Key Projects of Public Scientific Research Institutes (2018-9031).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Train Hunting Related Fast Degradation of a Railway Crossing—Condition Monitoring and Numerical Verification**

#### **Xiangming Liu \* and Valéri L. Markine \***

Department of Engineering Structures, Delft University of Technology, 2628 CN Delft, The Netherlands **\*** Correspondence<EMAIL_ADDRESS>(X.L<EMAIL_ADDRESS>(V.L.M.)

Received: 4 March 2020; Accepted: 15 April 2020; Published: 17 April 2020

**Abstract:** This paper presents the investigation of the root causes of the fast degradation of a railway crossing. The dynamic performance of the crossing was assessed using the sensor-based crossing instrumentation, and the measurement results were verified using the multi-body system (MBS) vehicle-crossing model. Together with the field inspections, the measurement and simulation results indicate that the fast crossing degradation was caused by the high wheel-rail impact forces related to the hunting motion of the passing trains. Additionally, it was shown that the train hunting was activated by the track geometry misalignment in front of the crossing. The obtained results have not only explained the extreme values in the measured responses, but also shown that crossing degradation is not always caused by the problems in the crossing itself, but can also be caused by problems in the adjacent track structures. The findings of this study were implemented in the condition monitoring system for railway crossings, using which timely and correctly aimed maintenance actions can be performed.

**Keywords:** railway crossing; wheel-rail impact; train hunting; numerical verification; railway track maintenance

#### **1. Introduction**

In the railway track system, turnouts (switches and crossings) are essential components that allow trains to pass from one track to another. A standard railway turnout is composed of three main parts: switch panel, closure panel, and crossing panel, as shown in Figure 1. In a railway turnout, the crossing panel is featured to provide the flexibility for trains to pass in different routes.

**Figure 1.** Standard left-hand railway turnout and the definition of the passing routes.

For rigid crossings that are commonly used in conventional railway lines, the gap between the wing rail and the nose rail usually results in high wheel-rail impacts in the transition region where the wheel load transits from the wing rail to the nose rail (vice versa, Figure 2), which makes the crossing a vulnerable spot in the railway track. In the case of crossings that are mainly used for the through route traffic (e.g., crossings in the crossover), there is no specific speed limit [1] and trains can pass through the crossings with a high velocity of up to 140 km/h. The high train velocity makes the wheel-rail impact more serious. In the Dutch railway system, around 100 crossings are urgently replaced every year [2] due to unexpected fatal defects, which not only result in substantial maintenance efforts, but also lead to traffic disruption and can even affect traffic safety.

**Figure 2.** Wheel-rail interaction in the railway crossing for through route traffic.

In contrast to a switch panel, wherein sensors are instrumented for condition monitoring [3,4] and remaining useful life prediction [5], monitoring in a crossing panel is usually absent. As a result, the real-time information on the condition of railway crossings is limited. The present maintenance activities are mainly reactive and based on the experience of the contractors. In this case, the root causes of the crossing degradation are not always resolved by the maintenance actions, and the crossings are likely to be operated in a degraded condition. To improve this situation, necessary guidance for maintenance actions is highly required.

Proper crossing maintenance usually relies on condition assessment and degradation detection, which can be realized through field monitoring. In recent years, condition monitoring techniques have been frequently applied in the railway industry. Aside from the above-mentioned instrumentation on the turnout switches, vehicle-based monitoring systems have been applied in track stiffness measurement [6] and estimation [7], track alignment estimation [8], hanging sleepers detection [9], and track fault detection [10], etc. Compared with the normal track, the current studies on railway crossings are mainly based on numerical simulation. Typical contributions include wheel-rail interaction analysis [11–21], damage analysis [16,17,22,23], and prediction [18,24,25] as well as crossing geometry and track stiffness optimization for better dynamic performance [16,26]. Field measurements are mainly used for the validation of numerical models. The monitoring of railway crossings for condition assessment and degraded component detection is still limited.

In the previous study, key indicators for the crossing condition assessment based on the field measurement were proposed [27,28]. Additionally, a numerical vehicle-crossing model was developed using a multi-body system (MBS) method to provide the fundamental basis for the condition indicators [29]. In this study, the condition indicators, as well as the MBS model, were applied in the condition monitoring of a fast degraded railway crossing. The main goals of this study were to investigate the root causes of the crossing degradation as well as to assess the effectiveness of the current maintenance actions.

Based on the objectives, this paper is presented in the following order. The experimental and numerical tools, including the crossing condition indicators, are briefly introduced in Section 2. The measurement results and the crossing degradation analysis as well as the effectiveness of the current maintenance actions are presented in Sections 3 and 4. Based on the measurement results and field inspections, the root causes for the fast crossing degradation were investigated with the assistance of the MBS model, as presented in Section 5. In Section 6, the verification of the effectiveness of the maintenance actions is given. Finally, in Section 7, major conclusions are provided.

#### **2. Methodology**

In this section, the experimental tools for the crossing condition monitoring, as well as the indicators for the crossing condition assessment, are briefly introduced. The MBS vehicle-crossing model for the verification of the experimental findings is also presented.

#### *2.1. Experimental Tools*

The experimental tools mainly consisted of the in-site instrumentation system modified from ESAH-M (Elektronische System Analyse Herzstückbereich-Mobil) and the video gauge system (VGS) for wayside monitoring, as briefly described below. Both tools have already been introduced and actively applied in previous studies. Detailed information regarding the installation and data processing can be found in [27,30].

#### 2.1.1. Crossing Instrumentation

The main components of the crossing instrumentation are an accelerometer attached to the crossing nose rail for 3-D acceleration measurement, a pair of inductive sensors attached in the closure panel for train detection as well as train velocity calculation, and the main unit for data collection. An overview of the instrumented crossing is shown in Figure 3.

**Figure 3.** Crossing instrumentation based on ESAH-M.

The main outputs of the crossing instrumentation were the dynamic responses of the crossing nose, including the wheel-rail impact accelerations and locations, etc. All these responses were calculated within the transition region, which can be obtained through field inspection [29]. Based on these measured responses and the correlation analysis between the responses [28], two critical condition indicators related to the wheel impact and fatigue area, respectively, were proposed.

The wheel impact is reflected by the vertical accelerations, which were obtained from the crossing and processed through statistical analysis. This indicator is mainly based on the magnitude of the impacts due to each passing wheel (Figure 4a), and the changes in time indicate the different condition stages of the crossing (Figure 4b).

**Figure 4.** Indicator for the wheel impact. (**a**) Procedure for the obtainment of wheel impacts. (**b**) Example of the variation of the wheel impacts in different condition stages.

The fatigue area is defined as the region where the majority of wheel impacts are located on the crossing, and where ultimately the crack initiates (Figure 5a). In practice, the fatigue area can be simplified as the confidence interval of [*a* − σ, *a* + σ], where *a* is the mean value of the wheel-rail impact locations, and σ is the standard deviation. The location and size of the fatigue area are critical values for the assessment of crossing wear and plastic deformation. A wide fatigue area usually represents well-maintained rail geometry. As demonstrated in Figure 5b, when the crossing condition was degraded from "Worn" to "Damaged", the fatigue area was dramatically narrowed and shifted further from the theoretical point (TP) of the crossing. More information about the fatigue area can be found in the previous study [27].

**Figure 5.** Demonstration of the crossing fatigue area detection. (**a**) Definition of the fatigue area. (**b**) Example of the fatigue area changes in different crossing condition stages.

#### 2.1.2. Wayside Monitoring System

The VGS for wayside monitoring is a remote measurement device based on digital image correlation (DIC). It uses high-speed digital cameras to measure the dynamic movements of the selected targets in the track. The system, set up together with the targets installed on the crossing rail next to the instrumented accelerometer, is shown in Figure 6a, and the demo of the displacement measurement is shown in Figure 6b. The main outputs are the vertical displacements of the tracked targets with a stable sampling frequency of up to 200 Hz.

**Figure 6.** Wayside monitoring. (**a**) System setup. (**b**) The screen of displacement measurements.

Due to the limitation of the experimental conditions, the wayside monitoring system is usually set up close by the side of the track, which will introduce extra noise in the measured displacement results. To improve the accuracy of the measurement, the noise part needs to be eliminated. The noise mainly comes from the ground-activated camera vibration, which can be manually activated by hammering the ground near the camera. The measured camera vibrations in both the time and frequency domains are given in Figure 7.

**Figure 7.** Ground activated camera vibration. (**a**) Time domain signal. (**b**) Frequency domain responses.

Despite the differences in the displacement responses in the two monitored crossings, the main resonance of the camera vibration was around 15–45 Hz. In the previous study [30], the main components in the displacement signal were elaborated. The train-track components related to displacement responses are mainly distributed below 10 Hz, which do not overlap with the camera vibration introduced noise. The noise part due to camera vibration can then be reduced through low-pass filtering, as shown in Figure 8.

**Figure 8.** Examples of the measured rail vertical displacement.

The magnitude of the dynamic vertical displacement of the rail directly reflects the intensity of the track movement due to the passing trains. By comparing the measured rail displacement with the reference level, which can be obtained from numerical simulation using the parameters in the designed condition, the ballast settlement level of the monitored location can be estimated. The MBS model for the crossing performance analysis is described later in this section.

#### *2.2. Multi-Body System (MBS) Vehicle-Crossing Model*

The numerical model for the crossing performance analysis was developed using the MBS method VI-Rail (Figure 9a). The rail pads, clips, and ballast were simulated as spring and damping elements (rail busing and base busing, Figure 9b). In the vehicle model, the car body, bogie frames and the wheelsets were modeled as rigid bodies with both the primary suspension and secondary suspension taken into account (Figure 9b). The track model was a straight line with the crossing panel (Figure 9c) situated in the middle of the track. The rail element for the acceleration and displacement extraction was the lumped rail mass located 0.3 m from the TP of the crossing (Figure 9d), which is consistent with the setup of the field measurements (Figures 3 and 6a).

**Figure 9.** Multi-body system (MBS) model. (**a**) Vehicle-track model. (**b**) Flexible connections in the model. (**c**) Crossing profiles. (**d**) Rail element for acceleration extraction.

The detailed model development, experimental validation, and numerical verification can be found in the previous study [29]. Corresponding to the condition indicators, the main outputs of the MBS model are the wheel impact acceleration, transition region and wheel-rail contact forces. Using the MBS model, the condition of the monitored crossing, as well as the detected track degradations, can be verified.

#### **3. Field Measurements and Analysis**

The monitored crossing was a cast manganese crossing with an angle of 1:9. As part of a crossover, trains mainly pass the crossing in the facing through route (Figure 2) with a velocity of around 140 km/h. The on-site view of the crossing is shown in Figure 10a. According to the maintenance record, this crossing was suffering from fast degradation with the service life of only around three years (18 years on average [2]). At the beginning of the condition monitoring, the damaged crossing was completely renovated.

**Figure 10.** Overview of the monitored crossing. (**a**) On-site view. (**b**) Sketch view.

Figure 10b gives a sketch view of the crossing, including the setup of the monitoring devices and the layout of the adjacent structures, especially the small bridge in front of the crossing. Considering that the bridge is located quite close to the monitored crossing, the performance of the crossing might be affected by the bridge, which will be discussed later.

The measurement results from the crossing instrumentation were based on multiple train passages in one monitoring day. For the wayside monitoring, one sufficient train passage is enough to estimate the ballast condition. To maximally reduce the influence of the vehicle-related variables, the selected results were restricted to the commonly operated VIRM trains with velocities of around 140 km/h.

#### *3.1. Wheel Impacts*

Based on the estimated transition regions, the wheel impact accelerations were calculated. The distribution of the wheel impacts due to multiple wheel passages is shown in Figure 11a. The sample size, in this case, was 78 passing wheels. It can be seen that the wheel impacts presented a bimodal distribution. Around 80% of the wheel impacts were below 50 g, while the remaining 20% of the wheel impacts were extremely high with a mean value of around 350 g. Such a polarized distribution of impacts indicates the highly unstable wheel-rail interaction in this crossing. It was demonstrated in a previous study [29] that for this type of railway crossing, the average level of the wheel impact is around 50 g, meaning that the 20% of high impacts of the monitored crossing are already more than seven times higher than the average impact level. It can be imagined that such high impacts will dramatically accelerate the degradation procedure of the crossing.

**Figure 11.** Vertical acceleration responses of the monitored crossings. (**a**) Distribution based on multiple train passages in one day. (**b**) Example of impacts due to one bogie.

An example of the impact acceleration response in the time-domain due to the first bogie of a VIRM train is shown in Figure 11b. It can be seen that for the two passing wheels from the same bogie, the impacts can be quite different. The impact due to the front wheel was up to 350 g, while the rear wheel activated vertical acceleration was only 20 g. It has to be noted that the high impacts were not always introduced by the front wheel, but appeared to have random occurrences. Such results further confirmed the instability of wheel-rail interaction at this crossing.

#### *3.2. Fatigue Area*

The measured fatigue area of the monitored crossing is presented in Figure 12. It can be seen that the wheel impacts were widely distributed at 0.22–0.38 m from the TP with the fatigue area size of 0.16 m. According to the previous study [28], the transition region (Figure 2) for this type of crossing is around 0.15–0.4 m. The fatigue area widely covered 64% of the transition region, which can be considered to be in line with the expectation of a new crossing profile. Such results further confirmed that the crossing rail was not worn or deformed.

**Figure 12.** Measured fatigue area of the monitored crossing.