Issue 
J. Eur. Opt. SocietyRapid Publ.
Volume 19, Number 1, 2023
EOSAM 2021



Article Number  10  
Number of page(s)  9  
DOI  https://doi.org/10.1051/jeos/2023006  
Published online  27 February 2023 
Research Article
Transient optical simulation by coupling elastic multibody systems, finite elements, and ray tracing
Institute of Engineering and Computational Mechanics, University of Stuttgart, Pfaffenwaldring 9, 70569 Stuttgart, Germany
^{*} Corresponding author: luzia.hahn@itm.unistuttgart.de
Received:
1
February
2022
Accepted:
7
February
2023
Transient dynamical–thermoelastic–optical system simulation is an important expansion of classical ray tracing through rigid, resting lenses because the operating performance of highprecision optical systems can be influenced 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 finite element analysis and gradientindex 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.
Key words: GRIN ray tracing / dn/dT / Thermoelasticoptical system / Dynamicaloptical system / Elastic multibody system / Finite element methods / Inhomogeneous media / Transient analysis
© The Author(s), published by EDP Sciences, 2023
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 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 dynamicalthermoelasticoptical 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 dynamical–thermoelastic–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 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.
Fig. 1 Overall workflow of a holistic dynamical–thermal–optical system analysis. 
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.
2 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.
2.1 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(1)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(2)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 parameters and [7, 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(3)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 RungeKutta 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.
2.2 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–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 the linear equations of motion for the elastic body read(4)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 enables the approximation of the nodal displacements q by the reduced coordinates with(5)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 bestknown 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 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 of the reference frame K and a linear elastic deformation 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 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 of the undeformed surface is calculated. As shown in step two in Figure 2, the new reference surface minimalizes the nodal deformations and enables a faster calculation of the raysurface intersection point . 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 by a continuous function. Therefore, Zernike polynomials [17, 18] are useful but other twodimensional 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, timedependent rigid body motions and elastic deformations can be considered in the simulation of optical systems.
Fig. 2 Transforming the motion and deformation of a surface into a useful form for the ray tracing. 
2.3 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 timedependend nodal temperatures and the thermal deformations 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 timedependent 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 simulation. According to [23] the change in the refraction index with the temperature in optical glass can be described by(7)Thereby, D_{0}, D_{1}, D_{2}, E_{0}, E_{1}, and λ_{k} are material dependent constants, λ is the wave length of the light, and 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 thermalmechanical 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 twodimensional polynomials for the radial component and onedimensional 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.
Fig. 3 Temperature and refraction index distribution within a lens. (a) Temperature distribution within a lens; (b) approximated refraction index within a lens using cylindrical polynomials. 
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 timedependent mechanical rigid body motions, mechanical elastic deformations, thermal elastic deformations, and thermally induced refraction index changes is possible.
2.4 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}, NeweulM^{2} ^{2} is used for the EMBS simulation, TheelaFEA [21] for the thermoelastic FE analysis, and OMSim^{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.
2.5 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.
Fig. 4 Optical system with thermal (red arrow) and mechanical (springs) excitation. 
Fig. 5 Ray tracing and geometrical image simulation result for the undisturbed optical system. 
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.
3 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 timestep 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.
Fig. 6 Results of the dynamical and the thermoelastic system analysis. 
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.
Fig. 7 Results of the ray tracing through the optical system with different distortions considered (colored rays) compared to the ideal ray paths in the undisturbed system (grey, dashed). (a) Ray tracing with the consideration of the dynamic distortions; (b) ray tracing with the consideration Δn_{ϑ}; (c) ray tracing with the consideration of the thermoelastic distortortions; (d) ray tracing with the consideration of dynamic and thermoelastic distortions. 
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.
Fig. 8 Ray path within the middle lens without the consideration of Δn_{ϑ} (grey, dashed) and with its consideration (red). 
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 highprecision 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.
Fig. 9 Results of the geometrical image simulation with different distortions considered compared to the geometrical image simulation of the undisturbed system (grey). (a) Grey: no excitation blue: mechanical deformation; (b) grey: no excitation red: Δn_{ϑ} pink: Δn_{ϑ}, therm. deformation; (c) grey: no excitation violet: mechanical deformation, Δn_{ϑ}, thermal deformation. 
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 highprecision 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 highprecision 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.
4 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 highprecision 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.
Competing interests
The authors declare that they have no competing interests.
Funding
No funding has been applied or used.
Author’s contributions
The authors are the sole contributors to this work.
See https://www.itm.unistuttgart.de/en/research/dynamicalopticalsystems/ for a video sequence over a time period.
See https://www.itm.unistuttgart.de/en/research/dynamicalopticalsystems/ for a video sequence over a time period.
References
 Doyle K.B., Genberg V.L., Michels G.J., Bisson G.R. (2005), Optical modeling of finite element surface displacements using commercial software, vol. 5867, International Society for Optics and Photonics, SPIE, pp. 149–160. https://doi.org/10.1117/12.615336. [NASA ADS] [Google Scholar]
 Coronato P.A., Juergens R.C. (2003) Transferring FEA results to optics codes with zernikes: a review of techniques, vol. 5176, International Society for Optics and Photonics, SPIE, pp. 1–8. https://doi.org/10.1117/12.511199. [NASA ADS] [Google Scholar]
 Störkle J., Eberhard P. (2016) Using integrated multibody systems for dynamicaloptical simulations, in: Modeling, systems engineering, and project management for astronomy VII, vol. 9911, International Society for Optics and Photonics, pp. 542–556. https://doi.org/10.1117/12.2230692. [Google Scholar]
 Genberg V.L., Michels G.J., Doyle K.B. (2002) Making FEA results useful in optical analysis, vol. 4769, International Society for Optics and Photonics, SPIE, pp. 24–33. https://doi.org/10.1117/12.481187. [NASA ADS] [Google Scholar]
 Gatej A., Wasselowski J., Loosen P. (2012) Using adaptive weighted least squares approximation for coupling thermal and optical simulation, Appl. Opt. 51, 28, 6718–6725. [NASA ADS] [CrossRef] [Google Scholar]
 Doyle K.B., Genberg Victor, Michels G.J. (2002) Integrated optomechanical analysis, vol. 58, SPIE Press, Bellingham, Washington. [CrossRef] [Google Scholar]
 Wengert N. (2015) Gekoppelte dynamischoptische Simulation von Hochleistungsobjektiven (in German), in: Dissertation, Schriften aus dem Institut für Technische und Numerische Mechanik der Universität Stuttgart, vol. 40, Shaker Verlag, Aachen. [Google Scholar]
 Störkle J. (2018) Dynamic Simulation and Control of Optical Systems (in German), Dissertation, Schriften aus dem Institut für Technische und Numerische Mechanik der Universität Stuttgart, vol. 58, Shaker Verlag, Aachen. [Google Scholar]
 Meschede D. (2008) Optik, Licht und Laser (in German), 3rd ed., Vieweg+Teubner, Wiesbaden. [Google Scholar]
 Sharma A., Kumar D.V., Ghatak A.K. (1982) Tracing rays through gradedindex media: a new method, Appl. Opt. 21, 6, 984–987. [NASA ADS] [CrossRef] [Google Scholar]
 Hahn L., Störkle J., Eberhard P. (2019) Consideration of polarization during the ray tracing for mechanically stressed lenses in dynamicaloptical systems, Optik 193, 162923. [NASA ADS] [CrossRef] [Google Scholar]
 Hahn L., Eberhard P. (2021) Transient dynamicalthermaloptical system modeling and simulation, J. Eur. Opt. Soc. 17, 5. [CrossRef] [Google Scholar]
 Wengert N., Eberhard P. (2012) Using dynamic stress recovery to investigate stress effects in the dynamics of optical lenses, in: Proceedings of the 7th ICCSM, Zadar, Croatia. [Google Scholar]
 Nowakowski C., Fehr J., Fischer M., Eberhard P. (2012) Model order reduction in elastic multibody systems using the floating frame of reference formulation, Math. Model. 7, 40–48. [Google Scholar]
 Lehner M., Eberhard P. (2006) Modellreduktion in elastischen Mehrkörpersystemen (Model Reduction in Flexible Multibody Systems) 54, 4, 170–177. [Google Scholar]
 Schwertassek R., Wallrapp O. (1999) Dynamik flexibler mehrkörpersysteme (in German), Vieweg, Braunschweig. [CrossRef] [Google Scholar]
 Zernike F. (1934) Beugungstheorie des Schneidenverfahrens und seiner verbesserten Form, der Phasenkontrastmethode (in German), Physica 1, 7–12, 689–704. Elsevier B.V. [NASA ADS] [CrossRef] [Google Scholar]
 Doyle K.B., Genberg V.L., Michels G.J. (2012) Integrated otomechanical analysis, in: SPIE Press monograph, SPIE Press, Bellingham, Washington. [Google Scholar]
 Omar Z., Mitianoudis N., Stathaki T. (2010) Twodimensional Chebyshev polynomials for image fusion, in: 28th picture coding symposium, pp. 426–429. [CrossRef] [Google Scholar]
 Tangelder R.J., Beckmann L.H.J.F., Meijer J. (1993) Influence of temperature gradients on the performance of ZnSe lenses, in: Zuegge H. (ed) Lens and optical systems design, International Society for Optics and Photonics, pp. 193–201. [Google Scholar]
 Matter F., Ziegler P., Iroz I., Eberhard P. (2019) Simulation of thermoelastic problems with the finite element method, in: Proceedings in applied mathematics and mechanics, John Wiley & Sons, Inc. https://doi.org/10.1002/pamm.201900035. [Google Scholar]
 Nicholson D.W. (2003) Finite element analysis: thermomechanics of solids, CRC Press, Boca Raton. [CrossRef] [Google Scholar]
 Hoffmann H.J., Jochs W.W., Westenberger G. (1990) Dispersion formula for the thermooptic coefficient of optical glasses, in: Marker A.J. III (ed), Properties and characteristics of optical glass II, International Society for Optics and Photonics, pp. 219–230. [NASA ADS] [CrossRef] [Google Scholar]
 Bathe K.J. (1996) Finite element procedures, PrenticeHall, Upper Saddle River. [Google Scholar]
 Belytschko T., Liu W.K., Moran B. (2000) Nonlinear Finite Elements for Continua and Structures, John Wiley & Sons, Chichester. [Google Scholar]
All Figures
Fig. 1 Overall workflow of a holistic dynamical–thermal–optical system analysis. 

In the text 
Fig. 2 Transforming the motion and deformation of a surface into a useful form for the ray tracing. 

In the text 
Fig. 3 Temperature and refraction index distribution within a lens. (a) Temperature distribution within a lens; (b) approximated refraction index within a lens using cylindrical polynomials. 

In the text 
Fig. 4 Optical system with thermal (red arrow) and mechanical (springs) excitation. 

In the text 
Fig. 5 Ray tracing and geometrical image simulation result for the undisturbed optical system. 

In the text 
Fig. 6 Results of the dynamical and the thermoelastic system analysis. 

In the text 
Fig. 7 Results of the ray tracing through the optical system with different distortions considered (colored rays) compared to the ideal ray paths in the undisturbed system (grey, dashed). (a) Ray tracing with the consideration of the dynamic distortions; (b) ray tracing with the consideration Δn_{ϑ}; (c) ray tracing with the consideration of the thermoelastic distortortions; (d) ray tracing with the consideration of dynamic and thermoelastic distortions. 

In the text 
Fig. 8 Ray path within the middle lens without the consideration of Δn_{ϑ} (grey, dashed) and with its consideration (red). 

In the text 
Fig. 9 Results of the geometrical image simulation with different distortions considered compared to the geometrical image simulation of the undisturbed system (grey). (a) Grey: no excitation blue: mechanical deformation; (b) grey: no excitation red: Δn_{ϑ} pink: Δn_{ϑ}, therm. deformation; (c) grey: no excitation violet: mechanical deformation, Δn_{ϑ}, thermal deformation. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.