Multi-physics Coupling Optimization for Fixing Cubic Fabry-Pétro Cavity Based on Data Learning

1. Introduction

High-finesse Fabry-Pérto (FP) cavity is one of the most vital components of the Ultra-stable laser (USL), which has been widely used in several space missions, such as LISA (Laser Interferometer Space Antenna) [1], Taiji Program in Space [2], DECIGO (DECI-Hertz Interferometer Gravitational wave Observatory) [3], Post-GRACE (Gravity Recovery and Climate Experiment) [4], etc., and other high-precision fields [5,6,7,8,9,10].

The Pound–Drever–Hall (PDH) method is extensively employed in USLs, by which the continuous-wave lasers are locked to the resonance frequency of the FP cavity via high-speed and wide-band electronic control system [11,12]. The instability of a laser’s fractional frequency is entirely defined by the length instability of the FP cavity, the perturbations to this length must be minimized for spectral purity [13]. Therefore, the designation of the FP cavity must reach the demand of ultra-high precision coupling the multi-physic (Optical/Thermal/Mechanics) fields, including its shapes (cylindrical, cubic, spherical, multiple-bore, midplane, etc.), fixing strategy (expecting the deformation as small as possible), and vacuum temperature control (getter and ion pumps, 3 or more layers of heat shields, etc.) [14,15,16,17,18].

The cubic FP cavity has the advantage of being force-insensitive and thus appropriate for space missions [19]. Based on this idea, follow-up studies have developed several ground-transportable USLs [20]. A typical USL system based on a cubic FP cavity is shown in Figure 1 [21], The FP cavity made of ultra-low-expansion (ULE) glass is attached to a titanium bracket by four titanium screws, which support the truncated vertices symmetrically. Gaskets made of PEEK are placed between the screws and the spacer for heat and vibration insulation. In the center of the top and bottom surfaces, two mirrors (one plane and one concave) are optically attached. Subsequently, the bracket was affixed to the inner heat shield, and three heat shields altogether were employed, as depicted in Figure 1b, between which PEEK gaskets are inserted to prevent heat transfer. Finally, heat shields are secured to the vacuum chamber using titanium screws. The ion and getter pump maintained a high vacuum level below $2\times {10}^{-6}$ Pa. As evident from the system composition above, The FP cavity’s performance is significantly affected by the shape parameters of the cavity and the temperature fluctuation. However, as far as we know, there is currently no research that specifically focuses on optimizing the design parameters of the FP cavity in consideration of multi-physic field coupling, which this paper precisely addresses.

The paper is organized as follows. In Section II, the multi-physics coupling theory and the finite element method are introduced for problem modelling. The results of the cubic FP cavity performance are discussed in Section III, including displacement distribution and cavity length change under the different experiments. Section IV consists of the mechanical optimization process, through which the optimal combination of key design parameters is obtained. Finally, Section V presents the discussions and concluding remarks.

2. Multi-physics coupling theory and finite element method for FP cavity

2.1. Multi-physics coupling theory

Multiphysical field coupling plays an essential role in understanding and optimizing complex systems involving multiple interactive physical phenomena. It helps make informed decisions for improved performance and reliability. According to the conservation of energy, the total energy consists of the conductive heat and the radiation heat. The heat conduction equation for the temperature distribution of optical components over time is [22,23]

$$\rho c\frac{\partial T}{\partial t}=\kappa \left[\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}+\frac{{\partial }^{2}T}{\partial z^2}\right]+q_a$$

(1)

$\rho$, c and $\kappa$ are material density, specific heat capacity and thermal diffusivity, respectively. t is the time, $q_a$ is the laser energy absorbed per unit volume and can be represented as

$$q_a(x,y,z,t)=\beta (x,y,z)I(x,y,t)$$

(2)

Where $\beta (x,y,z)$ is the volume absorption rate, $I(x,y,t)$ is the distribution of laser intensity

$$I(x,y,t)=P_{0}f(\mathbf{o},\mathbf{e})$$

(3)

$P_0$ is the beam power, $f(\mathbf{o},\mathbf{e})$ is the distribution of the disposed beam, which is often considered Gaussian

$$f(\mathbf{o},\mathbf{e})=\frac{1}{2\pi {\sigma }^2}exp\left(-\frac{d^2}{2{\sigma }^2}\right),d=\frac{\parallel \mathbf{e}\times (\mathbf{x}-\mathbf{o})\parallel }{\parallel \mathbf{e}\parallel }$$

(4)

Where $\mathbf{o}$ is the beam origin point, $\mathbf{e}$ is beam orientation, $\mathbf{x}$ is point position, d is the squared distance from the beam axis, $\sigma$ is the standard deviation.

The laser energy absorbed by the surface is used as the thermal boundary condition, considering the thermal radiation of the surface, then yield

$${\left.\kappa \frac{\partial u(x,y,z,t)}{\partial z}\right|}_{\begin{array}{c}\partial \Omega \end{array}}=\eta I(x,y,t)+ϵ\sigma A(T_{amb}^4-T^4)$$

(5)

assuming that the initial temperature distribution is uniform

$${\left.u(x,y,z,t)\right|}_{t=0}=T_0$$

(6)

Where $\eta$ is the surface absorptivity, $T_0$ the initial distribution, $\partial \Omega$ the boundary of the part, $ϵ$ surface emissivity, $\sigma$ the Stefan–Boltzmann constant, A the surface area, $T_{amb}$ the ambient temperature, T the temperature of heat radiation sources.

Combine the thermal conduction equation 1 with the following Equation 7, finally the deformation of the component can be obtained.

$${\nabla }^2\mathbf{u}+\frac{1}{1-2\nu }\nabla (\nabla ·\mathbf{u})=\frac{2(1+\nu )}{1-2\nu }{\alpha }_T\nabla T$$

(7)

where $\mathbf{u}$ is the displacement field, $\nu$ is Poisson’s ratio, ${\alpha }_T$ is linear thermal expansion coefficient, T is the internal temperature distribution, ∇ is Hamilton operator.

2.2. Finite element method and model establishment

The finite element method (FEM), which origin in the early 1960s, is by far the most commonly used approach in numerical analysis and engineering fields [24]. Take a multidimensional steady-state heat conduction process as a model problem, which can be described by the Poisson equation with homogeneous boundary conditions

$$\begin{array}{c}-{\nabla }^{2}u=f\mathrm{in}\Omega ,\\ u=0\mathrm{on}\partial \Omega ,\end{array}$$

(8)

with domain $\Omega \subset {\mathbb{R}}^d$. where u is the unknown function to be solved, f is the known function (source term). The weak form of Equation (8) can be achieved by choosing a function v from a space U of smooth functions, then forming the inner product of both sides with v, i.e.,

$$-〈{\nabla }^{2}u,v〉=\langle f,v\rangle$$

(9)

Assume that besides having the necessary smoothness, the functions which to be the solutions also satisfy the boundary conditions, the space U of test functions is of the form $U=\left\{v:v\in C^2(\Omega ),v=0\right.$ on $\left.\partial \Omega \right\}$. More concretely, let d = 2, the weak form of Equation (8) is given by

$${\int \int }_{\Omega }\nabla u·\nabla vdxdy={\int \int }_{\Omega }fvdxdy.$$

(10)

To obtain a numerical method, U is required to be finite-dimensional with basis $\{u_1,\cdots ,u_n\}$. Then the approximate solution $u^h$ of Equation (8) can be represented as

$$u^h=\sum _{j=1}^n{c}_j{u}_j.$$

(11)

Once a basis has been chosen for the approximation space U, the next is to determine the coefficients $c_j$ in Equation (11). By inserting $u^h$ into the weak form Equation (10), and selecting as trial functions v the basis functions of U, a system of equations is obtained

$${\int \int }_{\Omega }\nabla \left[\sum _{j=1}^n{c}_j{u}_j\right]·\nabla u_{i}dxdy={\int \int }_{\Omega }f{u}_{i}dxdy,i=1,\dots ,n$$

(12)

Which is known as the Ritz-Galerkin method and can be written in matrix form, $Ac=b$, where A is the stiffness matrix

$$A_{i,j}={\int \int }_{\Omega }\nabla u_j·\nabla u_{i}dxdy$$

(13)

3. Results and analysis of multi-physics coupling

3.1. Displacement distribution of FP cavity under the fixing force

For the purpose of modelling the interaction between different physical fields, COMSOL multiphysics was adopted as the simulation platform. Taking into account the thermal effect of the laser beam in optical and thermal coupling, the magnitude of beam power $P_0$ in Equation (3) was set to ${10}^{-4}$ Watt. In thermal and mechanical coupling, the heat radiation to the environment (default ${20}^{\circ }$C) and heat conductivity were taken into consideration. Besides, assume that the temperature fluctuation of the vacuum chamber was ${10}^{-3}$ Kelvin. Then the finite element model was assigned materials, meshed, and imposed constraints in solid mechanical analysis. A fixed force of 200 N was applied to each of the four truncated symmetrical vertices. Finally, a virtual prototyping multi-physic coupling model was established and the results could be calculated by FEM.

The displacement and Von Mises stress distribution of the FP cavity under the fixing force is shown in Figure 2, which color map is scaled logarithmically for better visualization. It is evident that the displacement and stress amplitude are higher at the four supporting vertices, which have symmetrical distributions as a result of the given fixing strategy. The performance index that should be put the most emphasis is the cavity length change, that is, the displacement of the cavity along the Z axis. Therefore, the Z component displacement of the plane mirror, the concave mirror, and the cavity length change, i.e., the sum of the above two, are shown in Figure 3. It can be seen that the displacement of the two mirrors is symmetrical about the line of ${45}^{\circ }$ and -${45}^{\circ }$, respectively. In order to further study the variation characteristics of the cavity length change $\Delta L$, two stripe regions of data in Figure 3c with a width of 2 mm along the ${45}^{\circ }$ and -${45}^{\circ }$ are extracted separately. Each of them consists of a total of 40000 data points and is converted to a one-dimensional signal.

A noise-assisted data analysis method named Ensemble Empirical Mode Decomposition (EEMD) is introduced, which consists of an ensemble of workers, each worker performing Empirical Mode Decomposition (EMD) on a copy of the input signal with added noise. The mean of these EMD results is considered as the final result [25]. EMD recursively decomposes a non-stationary signal into data-dependent basis functions termed Intrinsic Mode Function (IMF) [26]. In the EEMD process, the white noise with a standard deviation of 0.1 is added for the calculation, and the number of realizations is set to 100. The EEMD results of cavity length change along ${45}^{\circ }$ and $-{45}^{\circ }$ are shown in Figure 4. It can be concluded that the data comprises a fundamental half-sinusoidal component with a low frequency (IMF6) and other components of higher frequencies. The spectral kurtosis (SK) is a common dimensionless time series statistic for detecting and characterizing non-stationarities in a signal that can reflect the random distribution of time series data. A high SK level corresponds to a high level of nonstationary or non-Gaussian behavior [27,28]. The spectral kurtosis of cavity length change $\Delta L$ along ${45}^{\circ }$ and $-{45}^{\circ }$ are computed and visualized by kurtogram respectively, as shown in Figure 5, for detecting and characterizing non-stationarities in a signal. The Kurtogram uses normalized frequency, i.e., the sample rate is set to 40000 Hz for time normalization to 1 second. The Kurtograms reveal that the maximum K value is higher in the $-{45}^{\circ }$ direction compared to the ${45}^{\circ }$ direction, with both directions exhibiting low levels at $K_{max}$, which indicates that the impact of random signal fluctuations with high frequencies is negligible.

3.2. Analysis of cavity length change under the random vibration experiment

During aerospace missions, it is imperative to perform vibration tests on the FP cavity, which is one of the weakest components of a USL. Table 1 displays the parameters for random vibration testing using the Long March 5 rocket platform [29]. The test is formed by three successive stages, i.e., rising, holding, and falling. Moreover, the parameters are more stringent than actual launch conditions to guarantee payload reliability.

The displacement power spectrum densities of cavity length change responses at different points along the circle of radius 5mm and 2.5mm are shown in Figure 6, respectively. Various PSD responses along the logarithmic frequency curves at eight points every ${45}^{\circ }$ around the circumference are presented, which first increase, then remain constant, and then decrease as the test progresses. The PSD value at ${180}^{\circ }$ is the highest of radius 5mm, while the lowest value occurs at ${315}^{\circ }$. By contrast, the PSD values along radius 2.5mm at ${45}^{\circ }$, ${135}^{\circ }$, ${225}^{\circ }$, ${315}^{\circ }$ are far more than others (almost zero), which indicates a more complicated $\Delta L$ distribution under the vibration test. Therefore, the parameter of 100 Hz which corresponds to the maximum PSD value of the FP cavity is chosen for the next analysis.

Figure 7 shows the displacement and Von Mises stress distribution under the frequency of 100 Hz, whose color map is logarithmically scaled. Interestingly, the variation of the distribution is non-uniform, i.e., they have local maxima in addition to the four supporting areas, especially on the mirror surface. The Z component displacement of the plane mirror, the concave mirror, and the cavity length change are shown in Figure 8. It can be observed that the cavity length displacement has local maxima and minima around the circumference at 2.5mm and 5mm radii. Similarly, two stripe regions of data in Figure 8c with a width of 2 mm along the ${45}^{\circ }$ and -${45}^{\circ }$ are extracted separately and converted to a one-dimensional signal.

The EEMD results of cavity length change along ${45}^{\circ }$ direction and $-{45}^{\circ }$ direction under the vibration test are shown in Figure 9. In the EEMD process, the white noise with a standard deviation of 0.6 is added for the calculation, and the number of realizations is set to 200. The IMFs are more complex and suffer from frequency aliasing (IMF3-6) compared to the fixed force condition’s cavity length change. The pattern of EEMD result is almost identical along both directions, since the two sets of data appear to be highly similar.

The Kurtogram of cavity length change $\Delta l$ data along ${45}^{\circ }$ and $-{45}^{\circ }$ under the vibration test are computed and visualized in Figure 10, respectively. There is a higher maximum K value in the -${45}^{\circ }$ direction compared to the ${45}^{\circ }$ direction, with both directions also exhibiting low levels at $K_{max}$, indicating that transient fluctuations of the signal have little influence even if it appears highly nonstationary.

4. Mechanical optimization of fixing cubic FP cavity

4.1. Determinations of design spaces and performance indexes

In Section 3 the variation regularity of cavity length is calculated and visualized through multi-physics coupling and signal processing method. The cavity length change is expected to be minimal in all circumstances, consequently, the max displacement of the cavity length $w_F$ under the supporting force, and the max displacement of cavity length $w_{vib}$ under the random vibration test are determined as performance indexes. Besides, the volume of cubic cavity V is decided as a performance index for the purpose of lightweight.

Apparently, the cavity length change is determined by mechanical structure. However, there are only a few parameters that can be varied to optimize the cubic FP cavity. As shown in Figure 11, the truncating length of the FP cubic edge l, the diameter of the center hole d, and the magnitude of the supporting Force F applied on four vertices are determined as the key design variables.

Considering a reasonable range of the design space in practical engineering, the optimization problem of the FP cavity can be expressed as

$$\begin{array}{cc}\hfill \underset{x}{min}& O\left(x\right)={[V\left(x\right),w_F\left(x\right),w_{vib}\left(x\right)]}^T\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& x_1=d\in [5,13]mm,\hfill \\ \hfill & x_2=l\in [20,45]mm,\hfill \\ \hfill & x_3=F\in [100,400]N.\hfill \end{array}$$

(14)

4.2. Orthogonal experiments design and implementation

The orthogonal experimental design, which has been widely used in the engineering optimization field [30], is an efficient and economical technique that uses uniform orthogonal combinations of the design parameters to represent the characteristics of the entire design space. For the optimization of the FP cavity, an orthogonal experiment containing 49 sets of the design parameters is conducted according to the design space in Equation (14), and the corresponding results are shown in Table 2.

4.3. Establishment and comparison of data learning models

In order to establish the nonlinear relationship between the design variables and the performance indexes, the neural network (NN) is initially selected as the fitting model. The working mechanism of a single neuron can be described by

$$h_{\theta }\left(x\right)=g(\sum _{i=1}^n{\omega }_i{x}_i+b)$$

(15)

Where $x_i$ is the input of the neural network, b is the bias, and $\omega$ is the corresponding weights, which are the parameters to be learned. $g(·)$ is the activation function, of which Sigmoid and ReLU are commonly used [31]. A network of neurons organized hierarchically is called a neural network. To minimize the loss function, the gradient descent method is applied by using forward propagation and backpropagation, as well as repeated iterations according to the chain rule.

The samples in Table 2 are divided into 90% for training and 10% for testing. The layer size is determined as 40 and Bayesian regulation is adopted as the training algorithm. Figure 12 illustrates the training performances of the neural network model after 1000 epochs. Clearly, the neural network with a very small fitting error (${10}^{-11}$) exhibits excellent performance.

To further compare the performance level of different fitting models, two surrogate models, i.e., the quadratic response surface model and the Kriging model were selected and established for comparison [30]. The complete quadratic response surface model (RSF) with 10 unknown parameters is defined by

$$y_{RSF}\left(x\right)=a_0+\sum _{i=1}^3{a}_i{x}_i+\sum _{j=1}^3{a}_j{x_j}^2+\sum _{i=1}^3\sum _{j=1,j\ne i}^3{a}_{ij}{x}_i{x}_j=\left[1,x,x^2,x_1{x}_2,x_1{x}_3,x_2{x}_3\right]·C$$

(16)

Where $y_{RSF}\left(x\right)$ is the quadratic response surface function, the regression coefficients $a_0$, $a_i$, $a_j$, $a_{ij}$ are the constant term, the primary term, the quadratic term, the cross terms of the complete quadratic response surface model, respectively. C is the coefficient matrix of the quadratic response surface function, and the response surface coefficients can be fitted by the least square method.

The Kriging response function $y_{KRG}\left(x\right)$ is comprised of two parts: a regression model and a Gaussian correlation model

$$y_{KRG}\left(x\right)=\beta f\left(x\right)+z\left(x\right)$$

(17)

Where $f\left(x\right)$ is the global basis function, $\beta$ is the regression coefficient. $z\left(x\right)$ is the Gaussian correlation function whose mathematical expectation and covariance are given by

$$\begin{array}{c}E\left(z\right(x\left)\right)=0\hfill \\ Cov(z\left(x_i\right),z\left(x_j\right))={\sigma }^{2}R(\theta ,x_i,x_j)\hfill \\ R(\theta ,x_i,x_j)={\displaystyle \prod _{d=1}^q}exp(-{\theta }^d{\left|x_i^d-x_j^d\right|}^{p^d})\hfill \end{array}$$

(18)

Where ${\sigma }^2$ is the variance, $\theta$ is the unknown parameter of association, $R(\theta ,x_i,x_j)$ is the correlation function of point $x_i$ and $x_j$. For a point x, the predictive value calculated by Kriging model is

$$\stackrel{\wedge }{y}\left(x\right)=\stackrel{\wedge }{\beta }f\left(x\right)+r^T\left(x\right)R^{-1}(Y-\stackrel{\wedge }{\beta }F)$$

(19)

Where $r^T\left(x\right)$ is the covariance matrix between the unknown and known point, Y is the target value matrix of the sample point, and F represents the basis function matrix of the known points.

The test set of 5 points is used for model validation, and the comparisons between NN and surrogate models are shown in Figure 13. The three models perform quite well overall, yet the prediction accuracy of the three models decreases for volume, displacement under fixing force, and displacement under vibration test. This is likely because the volume variation along the design parameters is the most linear among the three indexes, followed by $w_F$, and lastly $w_{vib}$. As the variation of the index goes more non-linear, predicting with models becomes increasingly challenging.

The root mean square error (RMSE) is adopted as performance indicator

$$RMSE=\frac{1}{m{\overline{y}}^{Exp}}\sqrt{\sum _{i=1}^m{(y_i^{Exp}-y_i^{RSM})}^2}$$

(20)

Where m is the sample number of the sample training set, $y_i^{Exp}$ is the experiment result, $y_i^{RSM}$ is the prediction at i-th point of response surface model, ${\overline{y}}^{Exp}$ is the average of experiment points. The fitting accuracy is higher when the RMSE of the fitting model is closer to 0.

The RMSE results of three prediction models are shown in Table 3. It is evident that the neural network model provides higher accuracy compared to other models, of which the KRG model has the largest error. Predicting model performance directly impacts the algorithm’s accuracy in searching for the optimal solution, therefore the neural network model will be used for the next optimization process.

4.4. Evolutionary algorithm optimizing and performance verification

An evolutionary algorithm (EA) is a type of stochastic optimization method inspired by the natural process, especially natural selection in biological evolution, such as genetic algorithm, particle swarm algorithm, ant colony algorithm and so on [32]. EA is used to solve complex optimization and search problems where traditional optimization techniques may struggle. As a typical and proven effective evolutionary algorithm, Non-Dominated Sorting Genetic Algorithm II [33] (NSGA-II), which uses a fast non-dominated sorting algorithm, sharing, elitism, and crowded comparison, is adopted to search for the optimal combination of design variables.

Specifically, the NSGA-II configuration starts with an initial population of 1000 and uses the previously well-trained NN model to continuously evaluate the population’s performance at each generation. Figure 14 displays the Pareto optimal front after 100 generations of evolution.

The solution set of the Pareto optimal front surface consists of optimal points that are mutually non-dominant, which means that each point in the Pareto optimal front may be the optimal solution. In most cases of engineering, practicality and convenience must also be taken into account. After subtly selection, the final optimal combination of design variables is determined as $\{5,32,250\}$. The multi-physic model with the optimal parameters is established and evaluated. Moreover, the comparison of performance before and after optimization is presented in Table 4. It is satisfying that the performance indices are all better than they were before optimization by fine-tune the value and combination of design variables, especially when it comes to the displacement under the fixing force and vibration test, both of which are decreased by more than 60%. The ultimate outcome validates the effectiveness of the proposed multi-physic coupling optimization strategy based on data learning.

5. Conclusion

In this paper, a multi-physics coupling optimization strategy for fixing cubic Fabry-Pétro cavity based on data learning is proposed. A multi-physic model for the cubic FP cavity is established and the performance is obtained by finite element analysis. The key performance indices(V, $w_F$, $w_F$) and key design variables (d, l, F) are determined considering the features of the FP cavity. The 49 sets of data by orthogonal experiment are acquired for the establishment and comparison of different data learning models (NN, RSF, KRG). The result turns out that the neural network has the best performance. Based on NSGA-II, the Pareto optimal front is obtained and the optimal combination of design variables is finally determined as $\{5,32,250\}$. The performance after optimization has proven to be a great improvement, of which the displacement under the fixing force and vibration test are decreased by more than 60%. The method we proposed shows great application potential, in the future, we anticipate applying this optimization strategy in various situations and fields.