Transient optical simulation by coupling elastic multibody systems, ﬁ nite elements, and ray tracing

. Transient dynamical – thermoelastic – optical system simulation is an important expansion of classical ray tracing through rigid, resting lenses because the operating performance of high-precision optical systems can be in ﬂ uenced by dynamical excitations or thermal gradients. In this paper an approach for an integrated optical system simulation using the coupling of elastic multibody system simulations, thermoelastic ﬁ nite element analysis and gradient-index ray tracing is presented. Transient mechanical rigid body motions and elastic deformations, thermally induced refraction index changes, and thermal elastic deformations can be considered simultaneously in the ray tracing using the presented method. The calculation of the dynamical and thermal disturbances, the data transfer and coupling, and the gradient index ray tracing method are introduced. Finally, the approach is applied on a transient triplet lens optical system and some investigation results are shown.


Introduction
High performance optics, like lithographic objectives, have a high susceptibility to mechanical vibrations or temperature changes. In order to accurately predict the behavior of optical systems in operation, transient dynamical and thermal disturbances should be considered in the optical system design. Dynamical excitations, e.g. from the waver motion system or from vibrations in the environment, lead to rigid body motions, elastic deformations, and mechanical stresses which influence the imaging performance. Besides the mechanical disturbances, thermal effects should be determined and taken into account in optical system simulation. They occure due to a heating of the optical elements resulting from the energy absorption of the propagated light. Especially with the increasing power of the used lasers, the consideration of thermal effects is becoming more important. Temperature dependent changes of the refraction index in the used optical materials, thermal elastic deformations and thermomechanic stresses result from the thermal distortions. These effects have to be calculated and considered in the optical simulation. In this work a holistic simulation method for the analysis of transient dynamical-thermoelastic-optical systems is introduced. It considers time dependent rigid body motions, elastic mechanical and thermal deformations and thermally induced refraction index changes in the optical simulation where gradient index (GRIN) ray tracing is used. Additionally, a numerical example is presented which shows the application of the method on a transient dynamicalthermoelastic-optical system.
The behavior prediction of dynamically and thermally distorted optical systems requires the mechanical and the thermoelastic system analysis before the optical analysis can be performed. The mechanical analysis is often performed using finite element (FE) analysis [1,2]. The presented method uses elastic multibody system (EMBS) simulation for the calculation of the transient rigid body motions and the elastic deformations. Thereby, the system can be divided in rigid and elastic parts and model order reduction (MOR) methods can be applied [3]. This leads to a more efficient dynamical simulation with still a sufficient accuracy compared to transient mechanical FE simulations.
For the thermoelastic analysis an FE simulation is used. This is described in many papers and the resulting refraction index changes are then transferred to the optical system either using shape functions [4,5] or optical path difference maps [6]. Thereby, mostly commercial programs are used for the analysis. In contrast to this, the introduced method does not use a commercial program for the thermoelastic FE analysis. This enables more control over the coupling of the thermal and the mechanical effects. Besides that, in the presented method cylindrical polynomial

EOSAM 2021
Guest editors: Concita Sibilia, Alessandro Belardini and Gilles Pauliat approximations are used to describe the spatial refraction index distribution. Thereby, arbitrary three dimensional refraction index fields can be described and the deformations resulting from the EMBS analysis can be considered.
Subsequently to the dynamical analysis, the thermoelastic analysis, and the result data transformation, the presented method includes a GRIN ray tracing procedure. Thereby, the previously calculated system distortions are considered in the optical simulation.
The novel contribution of the presented method is the combination of the EMBS simulation with the thermoelastic FE analysis and the GRIN ray tracing and the detailled explanation of all necessary result data transformation. This allows the simultaneous consideration of transient dynamical and thermal elastic deformations in the ray tracing. Additionally, it enables the efficient consideration of transient mechanical deformations in the calculation of the thermally induced position dependency of the refraction index and its transient changes. A workflow for the presented method is given in Figure 1.
In the following, the transient system simulation approach is introduced. Thereby, the GRIN ray tracing, the dynamical analysis and the thermoelastic analysis together with the respective result data transformation are presented. Afterwards, an exemplary transient holistic optical system simulation is performed using the introduced method.

Methods and numerical example
In the design of optical systems, thermal and dynamical disturbances should be considered in order to enable accurate performance prediction over the operating period. In the following, a holistic simulation procedure for transient optical systems which considers dynamical and thermoelastic disturbances is presented in detail. First the optical simulation method is introduced. Then, the dynamical and the thermoelastic simulation methods together with the result data transformation are explained. Finally, a numerical example for the application of the presented method is introduced.

Ray tracing in GRIN media
According to Fermat's Law a light ray always takes the path s opt with the shortest time of travel between two points. This can be transferred into a minimum for the light path between two points. This so called optical path length along a trajectory S with a location dependent refraction index n(r) is defined as In homogeneous media, the refraction index is independent from the ray position. Thus, the ray path is a straight line which only changes its direction at the entry in a medium with another refraction index. In optical systems with homogeneous media but deformed optical surfaces, the light rays can be traced from surface to surface from the object to the image plane. At each deformed surface S i a ray position r i and a new ray direction d i can be calculated using and with the distance a i between the previous and the current intersection point, the refraction index of the previous medium n iÀ1 and of the current medium n i , the normal vector n i of the deformed surface, and the intermediate  8]. The calculation of the ray intersection point r i with the current deformed surface S i needs to be realized using an iteration to compute the parameter a i of equation (2). The smaller the difference between the deformed and the ideal surface the faster the intersection point is found. A ray in a medium with an inhomogeneous refraction index n(r) changes its direction within the medium. The partial differential ray equation describes the ray propagation. As shown in [9], this equation can be derived from equation (1) and therefore fullfills Fermat's Law. To track a ray through a medium with a gradient in the refraction index, equation (3) must be solved numerically. This so called GRIN ray tracing can be performed using a fourth order Runge-Kutta method [10]. A more detailed description of the ray tracing in homogeneous and in GRIN media is given in [7,8,11].
In order to trace a ray through a transient thermally and dynamically excitated optical system, accurate knowledge about the optical surface positions, the surface normals, and also the surface shapes over the time is essential. Besides that, the position dependent refraction index must be known for all time steps. These information is needed for the correct calculation of the ray positions and directions in equations (2) and (3). The transient surface shapes and the spatial distribution of the refraction index changes should be described by continuous functions in order to avoid discretization errors. The calculation of the dynamical and the thermoelastic system behavior and the transformation of the resulting data to a useful form for the optical simulation is introduced in the following.

Dynamical analysis using elastic multibody system simulation
Dynamical excitations of optical systems can occur, e.g. if mechanical ground vibrations exist instead of a completely static environment. This affects the performance of optical systems. Therefore, dynamical effects have to be considered in optical system analysis by simulating the system mechanically and transforming the results to the optical simulation. In the following, a method for the efficient numerical analysis of transient mechanical, i.e. dynamical effects in optical systems using EMBS simulation is presented. Thereby, transient rigid body motions, and elastic deformations are calculated. The mechanical stresses resulting from the deformations are neglected in this paper, because their influence is small compared to the rigid body motion and the elastic deformation. Possibilities of the consideration of mechanical stresses are given in [11][12][13].
In the beginning of the EMBS simulation, the optical system is divided into parts where only rigid body motion without deformation is expected and parts where deformations must be considered. The optically sensitive elements like lenses must be considered as elastic while parts of the holdings can be considered as rigid bodies. For each elastic body an FE model is created in a commercial FE program and a modal analysis is performed with this model. So, the system matrices, the eigenvalues and the eigenvectors are accessible. With the mass matrix M, the damping matrix D, the stiffness matrix K, and the vector of the time dependent nodal displacements qðtÞ 2 R M the linear equations of motion for the elastic body read The vector f contains the constraints and the mechanical excitations which act on the body. Accurate results of the mechanical analysis require a fine FE mesh, which can lead to an enormous number of elements and nodes and so to a high number of degrees of freedom M. This results in high computation times for the time integration of equation (4). Thus, model order reduction methods are applied in order to calculate the dynamical behavior with a much smaller equation system. Within a MOR a projection matrix V 2 R M Âm enables the approximation of the nodal displacements q by the reduced coordinates " The so called reduction basis V is also applied on all system matrices and the excitation vector. Thus, we get a system of equations with m ( M degrees of freedom.
There are various methods of identifying the reduction matrix [8,14,15]. The best-known method is the modal truncation. Thereby, the eigenfrequencies with their eigenmodes v are ordered with respect to rising frequencies and the m lowest modes are used in The eigenmodes are also called mode shapes, because they represent the bodies shape if it oscillates with its eigenfrequency. With the modal MOR, the elastic deformations of the elastic body are approximated for each time step by a weighted superposition of the eigenmodes where the calculated reduced coordinates " qðtÞ are the weights. After the reduction the elastic bodies, the rigid modeled bodies and the applied forces are assembled into a global EMBS using the floating frame of reference approach [16]. With a numerical integration of the EMBS in the time domain, the dynamical behavior of the simulated system results. The transient motion and deformation for each body is then described by a rigid body motion q r ðtÞ of the reference frame K and a linear elastic deformation " q e ðtÞ relative to the body reference frame.
In order to get the results of the dynamical simulation in an efficient usable form for the ray tracing, the rigid body motion and the elastic deformation must be transferred in coordinates relative to the previous optical surface S i and the deformation must be represented by a continuous description. The necessary transformation steps are visualized in Figure 2. For the transformation of the rigid body motion and the elastic deformation into the relative surface description, transformation matrices are applied on the rigid body motion and on the undeformed FE mesh of each body. The transformation matrices describe the translational and rotational correlation between the coordinate systems K and S i . After that, the elastic deformation is evaluated using equation (5). The nodal deformations z sag for a specific time step relative to the surface coordinate system S 0 i result. The sum of the transformed rigid body motion and the nodal deformation describes the deformed surface for a specific time step. The transformation is shown as step one in Figure 2. As described in Section 2.1, it is advantageous if the difference between the reference and the deformed surface is minimal. Therefore, an additional rigid body motion q 00 r of the undeformed surface is calculated. As shown in step two in Figure 2, the new reference surface minimalizes the nodal deformations z sag and enables a faster calculation of the ray-surface intersection point r i . The calculation of the reference surface is shown in [7]. The last step in the transformation of the dynamical analysis results is the approximation of the deformation z sag by a continuous function. Therefore, Zernike polynomials [17,18] are useful but other two-dimensional polynomials like linear combinations of Chebyshev polynomials [19] can be used as well. These polynomials can then be evaluated efficiently during the ray tracing. Consequently, with the use of the EMBS simulation and the transformation of the result data, time-dependent rigid body motions and elastic deformations can be considered in the simulation of optical systems.

Thermoelastic analysis using finite element methods
Thermal excitation of optical systems occurs due to the absorption of the propagating light by the optical elements like lenses. This results in an inhomogeneous temperature distribution within the lenses, which leads to inhomogeneous changes in the refraction index and to thermal deformations. These effects should be considered in the ray tracing for accurate optical simulation results. Thermomechanic stress effects within the optical elements are neglected in this work, since the other thermal effects predominate [20]. The calculation of the thermoelastic behavior of optical elements and the transformation of the results in a useful form for the ray tracing is presented in the following. For the thermoelastic simulation of an optical system, an FE model for each thermally relevant body is built. Usually, the used thermoelastic FE mesh is coarser than the mechanical FE mesh. This saves computing time and the results are sufficiently accurate, nevertheless. According to [21,22], the time-dependend nodal temperatures #ðtÞ and the thermal deformations uðtÞ can be calculated with the fully coupled thermoelastic equation system ð6Þ It consists of the equation of motion for elastic bodies, cf. equation (4), and the equation for heat conduction in solids. The thermal system matrices are the specific heat matrix D tt and the thermal conductivity matrix K tt . The thermal excitation is defined by the heat flux vector /. The equations are coupled by the thermoelastic stiffness matrix K st and the thermoelastic damping matrix D ts . The thermoelastic stiffness matrix describes the volumetric expansion of the material due to temperature changes and the thermoelastic damping matrix describes a change in the temperature due to the deformations in the material. The latter effect is usually extremely small and therefore D ts = 0 can be assumed. With a numerical time integration of equation (6) the time-dependent temperatures and deformations for each node result. In order to consider these results in the optical ray tracing, they must be transformed into relative surface coordinates and into a continuous description. The thermal deformations u, similarly to the mechanical deformations, can be transformed and approximated by continuous polynomial functions, see Section 2.2. The thermal surface deformations are then considerable via a polynomial evaluation during the ray tracing.
The eigenvalues of the thermal and the mechanical system are in completely different frequency ranges. Therefore, it is valid to calculate the thermal and the dynamical surface deformations separately and simply add the results when the deformed surface is calculated during the ray tracing. With this addition during the ray tracing both, thermal and dynamical deformations, can be considered in the calculation of the intersection between a ray and a deformed optical surface.
The inhomogeneous temperature within optical elements due to the energy absorption affects the refraction index and should therefore be considered in the optical J. Eur. Opt. Society-Rapid Publ. 19, 10 (2023) simulation. According to [23] the change in the refraction index with the temperature in optical glass can be described by Thereby, D 0 , D 1 , D 2 , E 0 , E 1 , and k k are material dependent constants, k is the wave length of the light, and n 0 is the refraction index at the reference temperature of 20°C. From equation (6), the transient temperatures of each node result and can be used for the evaluation of equation (7). So, the spatially disturbed refraction index changes within the optical elements are known for all dynamically undeformed nodes and all time steps. The consideration of the nodal refraction index changes in the dynamically and thermally disturbed optical system requires the spatial description of the nodal changes in the deformed system for each time step. As said at the beginning of this section, the mechanical FE mesh differs from the thermoelastic mesh. So, the mechanical deformations can not be added on the thermal nodes directly. An inverse distance weighting is used to apply the mechanical mode shapes given in the mechanical FE mesh on the thermoelastic FE mesh. Thereby, for each thermal node the nearest ten nodes of the mechanical mesh are determined. The corresponding parts of the mode shapes are then extracted and with the thermal-mechanical node distances as weights new shape functions for the thermoelastic mesh are calculated. With equation (5) the mechanical deformations can then be evaluated so that the transient mechanical deformation of the thermoelastic mesh results. Consequently, the spatial distribution of the refraction index changes are given within the dynamically deformed mesh. Now, the thermal deformations can be added and the transient refraction index changes in the thermal and dynamical deformed system result. In order to consider them in equation (3) in the GRIN ray tracing, the discrete results must be described by a continuous function. Therefore, the refraction index changes are approximated using cylindrical polynomials. They can be created by the linear combination of two-dimensional polynomials for the radial component and one-dimensional polynomials for the height component. For circular lenses, the combination of Zernike polynomials with Chebyshev polynomials can be suggested [7]. From the polynomial approximation a cylinder results within which the refraction index changes are described continuously. This is shown in Figure 3 for an example lens. In the shown case, the approximation cylinder has the height of the lens and the diameter of the transmissive part of the lens. No rays can run through the white lens parts of Figure 3b because they would be cut by the holdings. The polynomial description can be evaluated efficiently during the GRIN ray tracing.
With the use of transient thermoelastic FE analysis and the result data transformation, thermally induced timedependent deformations and inhomogeneous refraction index changes can be considered in the simulation of a dynamically and thermally disturbed optical system.
Altogether, using the procedure described in Sections 2.1-2.3, the simulation of the transient behavior of an optical system with the consideration of time-dependent mechanical rigid body motions, mechanical elastic deformations, thermal elastic deformations, and thermally induced refraction index changes is possible.

Implementation
The introduced method for the holistic simulation of transient dynamical-thermoelastic-optical systems is implemented mainly in MATLAB based software tools. The MOR is performed using MatMorembs 1 , Neweul-M 2 2 is used for the EMBS simulation, TheelaFEA [21] for the thermoelastic FE analysis, and OM-Sim 3 enables the data transformation and the optical simulation. These tools were developed at the Institute of Engineering and Computational Mechanics at the University of Stuttgart and are constantly being further expanded. Only the mesh generation and the calculation of the system matrices of equation (4) are performed using Ansys. The workflow for the holistic simulation is shown in Figure 1. During the development of TheelaFEA the generated system matrices of different example systems were compared to ABAQUS and Ansys in order to ensure the reliability of our softwaretool. Hexahedral elements with eight nodes and linear shape functions [21,22,24,25] are implemented in the softwaretool. Furthermore, TheelaFEA was validated by thermoelastic simulations of simple geometries which were compared to analytical solutions, and thermoelastic analyses for different meshings of lens geometries which were compared to thermoelastic FE Analyses with ABAQUS. Overall, the matrices and the results matched if no excessively distorted elements were used in the analyzed geometries.

Numerical example
In the following, the example system, shown in Figure 4, is introduced. On this system the presented transient dynamical-thermoelastic-optical system simulation approach is applied. The system consists of three lenses which are fixed in a holding cylinder. The lenses are connected to the cylinder via six beams in case of the outer two lenses and via three beams for the middle lens. At the points of contact between the lenses and the holdings no relative motion is possible for the mechanical model. Nearly the whole system is considered being mechanically elastic, only the bottom end of the cylinder is rigid, visualized by the black ring in Figure 4. There, the system is dynamically excited with a spring force. This is visualized with three springs over the system diameter. In Figure 4, two of the three springs are visible. In the beginning the system is deflected from the equilibrium in all spatial directions and then it oscillates freely. For the thermoelastic model, a constant heat flux excites the system thermally and deformations are blocked at the nodes of contact between the lenses and the holdings. For a better interpretation of the results, the heat flux is only applied on the middle lens. This is the thinnest lens and so there is the largest influence of the heat flux. The lens model consists of 1265 thermal and 3795 mechanical degrees of freedom in 936 linear hexahedral elements. The heat input is applied on all nodes within an imaginary cylinder in the middle lens which includes 145 of the nodes. The cylinder height corresponds to the lens thickness, the radius is one third of the lens radius and the center of the cylinder is one quarter lens radius shifted from the middle of the lens. Besides the heat flux, another thermal boundary condition is, that the temperature at the contact points between the holdings and the lens are fixed at room temperature. For the undisturbed system, a ray tracing result and the result of the geometrical image simulation are shown in Figure 5. It can be seen, that the rays focus on the image plane and that the image simulation leads to a regular grid.
The presented method for a transient simulation of dynamically and thermally disturbed optical systems is applied on this system and the results are discussed in the following.

Results and discussion
The results of the holistic optical system simulation with the consideration of transient dynamical and thermal distortions is presented in the following. In Figure 6, the excited system is shown for a single point in time. 4 On the left the result of the EMBS simulation is shown for a specific time step. The elastic deformations and rigid body motions are scaled, such that they become easily visible. In the shown time-step the middle lens is subject to large elastic deformations, as it can be seen on the left side of the lens. For the lower lens the rigid body motion dominates as it is nearly undeformed but rotated as a whole. The thermoelastic analysis result is shown on the right side of Figure 6. It shows the heating of the lens from the decentral heat input cylinder. The thermal excitation leads to the refraction index changes shown in Figure 3b.
In order to analyze the influence of the distortions on the optical system the dynamic and the thermoelastic analysis results are considered in the ray tracing shown in Figure 7. Thereby, different distortions are considered in each figure and each distortion is scaled in order to make the effects easily visible. It is always the same point in time shown and the lens colors represent the elastic deformation.
In Figure 7a only the results of the EMBS simulation were considered. The rigid body motion leads to a shift of the focal point closer to the lenses and due to a rigid body rotation of the last lens the focus is moved above the optical axis. The elastic deformations lead to a destruction of the focal point because of the different changes of the entrance and exit angles at each deformed surface for each ray.
In Figure 7b only the changes of the refraction index due to the thermal excitation is considered. One ray runs through the area where the refraction index is strongly affected and this is visible by the one ray that differs strongly from the ideal ray path. The path of this ray is shown in more detail in Figure 8. There, the influence of the thermally affected refraction index on the ray path gets   visible. The undisturbed ray (grey, dashed) runs in a straight line through the lens while the disturbed ray (red) proceeds in a curved path.
If thermal deformations are considered additionally to the thermally changed refraction index, the ray paths shown in Figure 7c result. The ray with the largest refraction index changes is still the most deflected one but in this case the thermal deformations change the angle of refraction such that the deviation of the ray from the ideal focus point decreases compared to Figure 7b.
In Figure 7d the rigid body motions and elastic deformations of the EMBS simulation, the thermally induced refraction index changes, and the thermal deformations are considered in the ray tracing. The resulting ray paths show the combination of all ray path distortions of the Figures 7a-7c. The focus point is destroyed, most of the ray intersections with the image plane are above the optical axis, and the strongly thermally influenced ray is still the most distorted one. The ray paths show that the dynamical and the thermal system distortions affect the ray paths and should therefore be considered in optical simulations if high-precision results are required.
Another possibility to show and analyze the influence of the distortions on the ray paths is presented in Figure 9. Thereby, a geometrical image simulation of a regular grid is performed. This means, a large number of rays are traced through the disturbed optical system with starting directions and positions such that the undisturbed system creates a regular grid on the image plane. Then, the   intersections of the rays with the image plane are plotted for a single point in time. In each geometrical image simulation other distortions are considered. In Figure 9a the undisturbed image is compared to the mechanically influenced system. Due to the rigid body motion the size of the image is reduced and the middle of the image is shifted. The distortion of the image visualizes the influence of the elastic deformations.
In Figure 9b, the influence of the thermal system excitation is visible. The decentral heat input leads to a higher temperature within the middle lens out of the lens center. This position of the high temperatures and so the position of the area with high refraction index changes is observable at the thermally influenced geometrical image. The upper right part of the image is strongly distorted, while the lower left part is nearly unchanged. There, the temperature of the lens stays close to room temperature and so the grid does not get distorted. If thermal deformations are considered additionally to the thermally induced refraction index changes, the distortion of the upper right corner increases. This visualizes, that the thermal deformations occur mainly in the strongly heated parts of the lens. The disturbances and their severity are dependent on the mechanical and thermal excitation magnitude. It is possible to scale the effects by either scaling the excitation or scaling the simulation results. It is important to keep in mind, that only linear effects are representable using the introduced method. Applying a high thermal excitation on the system, we get the results shown in Figure 9b. The thermal deformations were scaled additionally in order to make the thermal deformation effects easily visible. So, the influence of the thermal distortion on the refraction index has a larger effect on the ray paths than the thermal deformation. In the comparison of the lens deformations which lead to the images shown in Figures 9a  and 9b, the thermal deformation is several dimensions smaller than the mechanical deformation. This amplitude difference together with the fact that the thermal and the mechanical eigenfrequencies are in different frequency ranges show, that the independent calculation of the two deformations described in Sections 2.2 and 2.3 is valid. Figure 9c shows the combined thermal and mechanical distortion effects for a single point in time. 5 The shift of the image due to the rigid body motion, the grid distortion due to the elastic deformation, and the thermal influence on the upper right corner are visible. The combination of the distortions shows that all influences affect the imaging quality and, therefore, they need to be considered in the holistic simulation of high-precision optical systems.
Altogether, Figures 7 and 9 visualize the influence of time dependent dynamical and thermal distortions on optical systems. They show that it is important to consider these effects if high-precision optics are analyzed. Besides that, it results that the time dependent rigid body motions, elastic deformations, thermally induced refraction index changes, and the thermal deformations can be determined and considered in the holistic simulation of transient dynamical-thermoelastic-optical systems with the introduced method.

Conclusion
In order to consider transient dynamical and thermal distortions in the simulation of optical systems, a numerical method which combines dynamical EMBS simulation, thermoelastic FE analysis, and GRIN ray tracing has been presented and applied on a numerical example in this paper. During the procedure, transient rigid body motions and elastic deformations are calculated using EMBS simulation. Thermally induced transient refraction index changes and thermal deformations are determined using thermoelastic FE analysis. The resulting motions, deformations and refraction index changes are then transferred to relative surface coordinates and described using polynomial approximations. Then, the distortions are considered in a numeric GRIN ray tracing procedure. In an example of a triplet lens system the introduced method was applied. The resulting ray paths and geometrical image simulations show the importance of the consideration of the disturbances in the holistic simulation of transient high-precision optical systems. Mechanical disturbances like rigid body motion and elastic deformation as well as thermal disturbances like thermally induced refraction index changes and thermal deformations affect the ray paths and influence the imaging performance. The effects must be considered in the ray tracing if precise prediction of the operation performance is demanded. Besides illustrating how important it is to take these effects into account, the example shows that the presented method is capable of the calculation of transient dynamical and thermal effects of optical systems and it enables the consideration of these disturbances in the optical system simulation. So, the presented transient dynamical-thermoelastic-optical system simulation is an important expansion of classical ray tracing.