Numerical study using finite element method for the thermal response of fiber specklegram sensors with changes in the length of the sensing zone

The response of fiber specklegram sensors (FSSs) is given as function of variations in the intensity distribution of the modal interference pattern or speckle pattern induced by external disturbances. In the present work, the behavior of a FSS sensing scheme under thermal perturbations is studied by means of computational simulations of the speckle patterns. These simulations are generated by applying the finite element method (FEM) to the modal interference in optical fibers as a function of the thermal disturbance and the length of the sensing zone. A correlation analysis is per-formed on the images generated in the simulations to evaluate the dependence between the changes in the speckle pattern grains and the intensity of the applied disturbance. The numerical simulation shows how the building characteristic of the length of sensing zone, combined with image processing, can be manipulated to control the metrological performance of the sensors.


Introduction
The developments in FSS (Fiber Specklegram Sensor) that have been taking place in recent years, have positioned it as a valuable alternative for sensing multiple variables, given its characteristics such as small size, low comparative cost and its metrological properties that allow configuring different monitoring systems [1 -9]. These types of sensors are conceived as devices that respond to the alterations induced on a multimode fiber by different external physicochemical variables (temperature, pressure, deformation, concentrations) [10 -14] through changes in the intensity distribution in the modal interference pattern at the fiber output. At an experimental level, sensing schemes are developed using monomulti-mono mode fiber configurations, where the third component works as a filter that admits part of the light from the speckle pattern generated by the multimodal fiber, which allows configuring interrogated FSS schemes by variation of optical power. This technique is quite simple and inexpensive, although it misuses part of the spatial information of the speckle patterns associated with the state of the sensing fiber [15 -19].
To evaluate the response of the FSS sensor, various investigations are supported by digital image analysis, whose algorithms allow the analysis of multiple records and estimate the variations induced in the speckle patterns in the face of external disturbances [2, 20 -24]. Therefore, when alterations are presented in variables such as temperature, changes will be induced that can be estimated through correlation analysis carried out on the interference patterns from the sensing system. In this line, studies are reported that analyze the patterns captured by CCD cameras at the output of the multimode fiber, which are performed digital processing and correlation analysis in selected regions of said images [25,26].
An alternative to make the analysis of these optical sensing systems more efficient consists in the generation of multiple speckle patterns, corresponding to modal interference patterns, through computational models of propagation in optical waveguides under external disturbances [27], in order to have a large digital image bank for study, and thus establish the characteristics of the sensor system in the face of induced changes in speckle patterns. In this sense, the finite element method (FEM) is configured as a powerful simulation tool in the generation of optical systems solutions, in particular, for propagation and interference processes in FSS under controlled disturbances [28,29]. In consideration of the above, this work is oriented to study the response of a FSS, throw the analysis of the speckle patterns generated by FEM simulations, when the sensor is subjected to focused thermal disturbances. The evaluation of the sensor response is established by means of a correlation analysis discriminated in regions within these simulated images.

Thermal disturbances representation
To obtain the speckle pattern at the end of a multimode optical fiber (MOF) under a given temperature, the propagation of the electric field E  of each mode must be calculated by solving the vector wave equation: where k 0 = 2 /  0 is the wave number in vacuum,  the wavelength in vacuum, and n is the refractive index of the fiber. This refractive index is given by where n x , n y and n z are found through the change of refractive indices [10]: 0 , with , , In equation (3), is the change in the refractive index due to thermo-mechanical disturbance (changes by material expansion due to temperature variations) and TO i n  is the change in refractive index due to thermooptical disturbance (changes due to variations in temperature). Then, using (3) we obtain: , 2 where C TO stands for the thermo-optic coefficient,  for the coefficient of linear thermal expansion, P 11 and P 12 for the strain-optic coefficients of the material, T 0 for the temperature of the material at the initial moment (room temperature), and n 0 for the refractive index at the initial moment where the material is isotropic. For fused silica at  = 632.8 nm, the C TO = 11.9 × 10 -6 /°C [29,30],  = 0.5 × 10 -6 /°C, P 11 = 0.113, and P 12 = 0.252 [31]. In the core, n 0 = n 0co = 1.4570, calculated using the Sellmeier's equation [28]: For the cladding of the fiber n 0 = n 0cla = 1.4482. This value is calculated using equation (6) with a numerical aperture NA = 0.16 for the simulated fiber: Note that the terms in equations (4), where the strainoptic coefficients are contained, can be neglected because these are much smaller than C TO in the case of fused silica, obtaining the following simplification: Equation (7) shows that the index of fused silica exhibits approximately an isotropic behavior with changes in temperature.

Finite Element Method for numerical calculations of propagation modes
If we assume an optical fiber of infinite length with a refractive index invariant along the propagation axis, a 2.5D approximation of the problem is obtained. This approach consists of finding the state of the system in a plane and then propagating it in the normal direction. In this case, the field can be expressed as Here, the electric field   E , , x y z  of each mode can be calculated by solving equation (1) numerically through the finite element method (FEM) (see more details in [27,28,33]). To calculate the field   E , x y  by FEM, the cross section of the fiber in a reference plane, that is taken as z = 0, is discretized into triangular elements. Equation (1) is solved in each of these elements and then they are assembled to obtain the general solution for each mode (see fig. 1a). In this way, we obtain the vector field   E , x y  and the propagation constants of each mode in the plane z = 0. Afterward, each mode is propagated along a given z, by computing the expression . The vectorial su-perposition of all these modes produces the speckle pattern in the required z-plane [28].

Analysis of specklegrams
Since the changes in intensity due to thermal disturbances are visible at a local scale approximately close to that of the grain size of the speckle pattern, these changes cannot be properly quantified by averaging the intensity of the entire pattern. Therefore, when there are speckle patterns at the output of a MOF, correlation techniques are usually employed to quantify the changes due to these disturbances. This technique is carried out through the algebraic implementation of the correlation coefficient given by [21]: where  I ref  is the average intensity of the speckle pattern at the fiber output in a reference state ( fig. 1b) and  I n  is the equivalent in a disturbed state ( fig. 1c).

Simulated system
The simulated system is an FSS with a mono-multi mode fibers arrangement, with the multimode fiber acting as transducer. The simulated multimode fiber of fused silica is subjected to thermal disturbances and variation in some of its metrological parameters.  The thermal effects and the light propagation for the generation of the speckle patterns in the MOF detailed above are simulated following the procedure described in sections 1 and 2, respectively. Here, a script in Matlab that integrates Comsol Multiphysics is implemented for the numerical solution of Eq. (1) by means of FEM (Finite Element Method), where first, the refractive index changes due to thermic disturbance are calculated through (7). Second, this computed refractive index is replaced in the wave equation (1) which is solved by FEM for the fields   E , , x y z  of the propagating modes support by the MOF in the plane z = 0, that is where the zone of thermal disturbance begins (see figs. 1a, 2). Third, these computed fields are propagated a distance z = Z. Fourth, all these computed modes are vectorially added at the end of the MOF, producing the speckle pattern, as detailed in section 2. Finally, this specklegram at the fiber output are analyzed using the correlation coefficient defined in equation (8). Fig. 3 illustrates the effects that the variation of NA has on the resulting specklegrams. It shows that by increasing the value of NA the number of speckle grains in-creases significantly, as more modes can enter and propagate in the fiber. This is an example of the potential of these simulations to study the behavior of the sensor by modifying certain parameters. However, in this work we are centeaz the influence of the length of the sensing zone on the metrological parameters of the FSS. For this analysis we have selected an NA of 0.16. Fig. 4 shows the correlation curves when changing the thermal disturbance, for different lengths of the sensing zone (Z). Fig. 4a. compares the correlation curves and shows how the correlation coefficient undergoes through a somehow periodic behavior as temperature changes, for five different Z lengths. As the Z length increases, the pseudo period of these curves reduces, showing more sensitivity, i.e., faster changes in correlation. Additionally, these curves seem as having an attenuating envelope.

Specklegram analysis and sensor response
As a plausible interpretation of this behavior, we think that this evolution of the correlation curve is governed by how the changes of the index of refraction as temperature changes affect the different modes propagating in the fiber. Dividing the propagating modes in the fiber into two broad groups of modes, i.e., one group of modes where the refractive index change will only represent a phase change, and another group where index change will represent an exclusion or inclusion as a propagating mode. Thus, modes that are only affected in phase would account for the periodic appearance of the curve, as the intensity of their superposition would be changing periodically and so the correlation curve. On the other hand, the other modes would account for the apparent attenuation on the correlation curve, as their inclusion or exclusion would affect the intensity of the specklegram in an uncorrelated way with respect to the reference.
Also, accordingly to the former interpretation, for an increase in Z, the phase variation of the propagating modes with the change of temperature would be faster because of the cumulative effect, thus reducing the apparent period of the corresponding correlation curve. Fig. 4b. to 4f. show the correlation curves individually. Isolating them allows an easy visualization of the dynamic range and sensitivity for the linear parts of these curves. It can be inferred from these figures that by increasing the length of the sensor zone, not only does the sensitivity of the sensor increase, but as mentioned, it also means a fast change in the overall appearance of the specklegram reducing this zone of linear behavior. This resulting in a decreased dynamic range of the system. Dynamic range, sensitivity, and linearity are metrological parameters of great importance in the definition of sensor performance. These simulations allow us to ascertain these parameters in an easy way. In Table 1 we show the data obtained in the numerical calculations for each length of sensing zone. We applied a linear fitting to the data points in the valid interval and extracted the dynamic range as T = T max -T min , the sensitivity S as the slope of the linear fitting, and R 2 as a measure of the linearity. As mentioned, the data show more quantitatively that the dynamic range and sensitivity are strongly affected by the increase in Z, regularly increasing the sensitivity, but reducing the dynamic range. For Z = 1 mm there is a high sensitivity, even higher than for the two following Z lengths. This shows that these metrological values are also dependent on the processing and selection of zone of correlation calculation. The R 2 values show a good linearity for the valid interval of data.
The correlation measure is dependent on how large the changes in the intensity distribution of the pattern are. Large changes in pattern intensity can result in behaviors such as pseudo periodic shapes shown before or discontinuity of correlation curves. One proposal to avoid these non-monotonic drawbacks is to redefine the areas where the calculation of the correlation curve is performed.

Selection of region of interest (ROI)
This strategy uses reduced areas and avoids the uncontrolled influence of all the pixels in the image, by selecting regions where the behavior of speckles is more monotonical. Since some of the ROIs improve the metrological parameters and others can yield the opposite result, the location and size of these ROIs cannot be predicted a priori and must be established by searching them through the speckle patterns. Fig. 5 shows the result of applying this strategy to the numerical results for 3 different lengths of sensing zone (Z). It is visible from Fig. 5a that there is an improvement in the dynamic range to the point that for the 3 Zs analized, they overpass the larger range in the curves in the Figure 4. It is also evident the change in the sensitivity with the use of ROIs. Fig. 5b illustrates the locations of the correspondent ROIs to each Z length, over the entire area of the speckle pattern. This shows that for an expected behavior, the ROI must be determined individually. Table 2. shows the data of sensitivity and linearity for the curves in Fig. 5a. Although, the values in sensitivity are dependent in the selected ROI, the tendency of higher sensitivities for higher Zs still holds. Also, the use of these ROIs yields good linearities represented in R 2 values over 0.9. As it is shown by the results for Z = 4 mm, the selection of ROIs where low quantity of speckle grains is present can also help to control dynamic range, but with the associated trade-off of low sensitivity.

Conclusions
We showed that the FEM method for the numerical simulation of physical-optical phenomena in optical fibers can be a powerful tool in the study and design of FSSs. The numerical simulation of thermal disturbances applied to an FSS allowed us to show in a convenient way how the building characteristic of the length of sensing zone, combined with image processing, can be manipulated to control the metrological performance of the sensors. Over the correlation curves we showed improvements in dynamic range, linearity, and sensitivity. This work opens the possibility to consolidate a costeffective methodology for designing this kind of sensors.