Deformation Analysis of the Glass Preform in the Progress of Precision Glass Molding for Fabricating Chalcogenide Glass Diffractive Optics with the Finite Element Method

Precision glass molding (PGM) technology is a cost-efficient process for the production of micro/nanostructured glass components with complex surface geometries. The stress distribution, surface profile, and reduced refractive index of the molded lens are based on the lens being fully formed. The process of the deformation of the glass preform is rarely discussed, especially in the case of multi-machining parameters in the experiment. The finite element method (FEM) was adopted to analyze the glass preform deformation. Due to the phenomenon of incomplete deformation of the glass preforms in the experiments, two groups of finite element simulations with different boundary conditions were carried out with MSC.Marc software, to reveal the relationship between the deformation progress and the parameters settings. Based on the simulation results, a glass preform deformation model was established. The error between the model result and the simulation result was less than 0.16. The establishment method of the glass preform deformation model and the established model can be used as a reference in efficiently optimizing PGM processing parameters when the designed lens has two different base radii of curvature.


Introduction
Glass optical components are widely used in the camera, computer, laser projection, biological, and consumer electronics fields because of their optical performance, mechanical stability, thermal stability, and high chemical resistance [1][2][3][4]. The hybrid aspherical diffractive singlet achromat design can be used to reduce chromatic aberration in compact optical systems. The applications are almost in the areas of bifocal contact lenses, beam shapers, and barcode reader diffusers [5].
Precision glass molding (PGM) technology is a cost-efficient process for the production of micro/nanostructured glass components with complex surface geometries [6][7][8]. It is a technology that deforms glass preform into the desired shape by use of molds at a high temperature in an inert or vacuum environment [3,[9][10][11].
Among optical glass, chalcogenide glasses (ChGs) are a group of materials that contain one or more of the chalcogen elements (S, Se, and Te) and group IIIA-VA elements (Ge, Ga, As, and Sb). ChG has a wide optical window that can be extended from the visible to far-infrared spectrum [12]. Compared to Ge or ZnSe, ChGs are amorphous and can be fabricated to the lenses with PGM, which provides a significant advantage over diamond turning due to its high-volume manufacture and lower cost. The densities and dn/dTs of ChGs are lower than those of Ge, thus having advantages in weight and athermalization. In z 2 = c 2 ·r 2 2 1 + 1 − (1 + k 2 )c 2 2 r 2 2 + B 2 r 2 2 + B 4 r 4 2 + B 6 r 6 2 + B 8 r 8 2 + B 10 r 10 2 + 1 1 − n C 1 r 2 2 + C 2 r 4 2 . (2) Micromachines 2021, 12, x FOR PEER REVIEW 3 of 15 Chalcogenide glass diffractive-refractive hybrid lenses are molded on a Toshiba precision glass molding machine (Model No. GMP-415V, Toshiba Machine Co., Ltd., Numazu, Japan), as shown in Figure 1. The diameter and thickness of the designed lens are 22 mm and 4 mm, respectively. The aspherical surface and the hybrid aspherical diffractive surface are respectively described with Equations (1) and (2).
= ⋅ ( ) In Equation (1), c1 is the curvature at the center of the profile, k1 is the conic constant, r1 is the radius of the lens, and A2, A4, A6, A8, and A10 are the even order coefficients. In Equation (2), c2 is the curvature at the center of the profile, k2 is the conic constant, r2 is the radius of the lens, n is the coefficient of refraction of the lens, and B2, B4, B6, B8, and B10 are the even order coefficients of the aspheric base. The material of the lens is Schott grade IRG26. The coefficients of the designed lens are listed in Table 1.  In Equation (1), c 1 is the curvature at the center of the profile, k 1 is the conic constant, r 1 is the radius of the lens, and A 2 , A 4 , A 6 , A 8 , and A 10 are the even order coefficients. In Equation (2), c 2 is the curvature at the center of the profile, k 2 is the conic constant, r 2 is the radius of the lens, n is the coefficient of refraction of the lens, and B 2 , B 4 , B 6 , B 8 , and B 10 are the even order coefficients of the aspheric base. The material of the lens is Schott grade IRG26. The coefficients of the designed lens are listed in Table 1.

Aspherical Surface
Hybrid Aspherical Diffractive Surface The glass preform is a biconvex spherical lens with a curve radius of 20 mm, a thickness of 7 mm, and a diameter of 20 mm. The diagram of the designed lens and the glass preform is shown in Figure 2. The glass preform is a biconvex spherical lens with a curve radius of 20 mm, a thickness of 7 mm, and a diameter of 20 mm. The diagram of the designed lens and the glass preform is shown in Figure 2. Rapidly solidified aluminum RSA-905 was selected as the mold core material. Three PGM experiments were carried out with the processing parameters listed in Table 2. The photographs of the molded lens are shown in Figure 3. The difference among the molded diffractive surfaces was obvious, while the difference among the molded aspherical surfaces was small.  Rapidly solidified aluminum RSA-905 was selected as the mold core material. Three PGM experiments were carried out with the processing parameters listed in Table 2. The photographs of the molded lens are shown in Figure 3. The difference among the molded diffractive surfaces was obvious, while the difference among the molded aspherical surfaces was small. The glass preform is a biconvex spherical lens with a curve radius of 20 mm, a thickness of 7 mm, and a diameter of 20 mm. The diagram of the designed lens and the glass preform is shown in Figure 2. Rapidly solidified aluminum RSA-905 was selected as the mold core material. Three PGM experiments were carried out with the processing parameters listed in Table 2. The photographs of the molded lens are shown in Figure 3. The difference among the molded diffractive surfaces was obvious, while the difference among the molded aspherical surfaces was small.  A compact flash card can be inserted into the card slot of the PGM machine to store the data detected during automatic operation. The sampled data is saved in a file at each cycle. Sampling was performed from the start to the end of a cycle. The sampling period was set as 9000 s. One thousand data items were collected within the sampling period.  A compact flash card can be inserted into the card slot of the PGM machine to store the data detected during automatic operation. The sampled data is saved in a file at each cycle. Sampling was performed from the start to the end of a cycle. The sampling period was set as 9000 s. One thousand data items were collected within the sampling period. That is, the sampling cycle is sampling period/1000 items, which is nine seconds. The waveforms of the force and displacement of the lower mold since the forming stage began are shown in Figure 4. A compact flash card can be inserted into the card slot of the PGM machine to store the data detected during automatic operation. The sampled data is saved in a file at each cycle. Sampling was performed from the start to the end of a cycle. The sampling period was set as 9000 s. One thousand data items were collected within the sampling period. That is, the sampling cycle is sampling period/1000 items, which is nine seconds. The waveforms of the force and displacement of the lower mold since the forming stage began are shown in Figure 4. The deformation of the glass preform is related to the relative position of the upper and lower molds. The upper mold was installed on the fixed upper axis of the molding system, while the lower molds were installed on the mobile lower axis. Therefore, there is a need to describe the displacement of the lower axis with time. The displacement data are processed with regression analysis with IBM SPSS Statistics software (25, SPSS Inc., Chicago, IL, USA). When the time is t, the displacement is yi for the i th experiment. Based on the results, the line types of the curve consists of linear, conic, and cubic curves. The three kinds of curves are expressed as Equations (3)-(5). The deformation of the glass preform is related to the relative position of the upper and lower molds. The upper mold was installed on the fixed upper axis of the molding system, while the lower molds were installed on the mobile lower axis. Therefore, there is a need to describe the displacement of the lower axis with time. The displacement data are processed with regression analysis with IBM SPSS Statistics software (25, SPSS Inc., Chicago, IL, USA). When the time is t, the displacement is y i for the ith experiment. Based on the results, the line types of the curve consists of linear, conic, and cubic curves. The three kinds of curves are expressed as Equations (3)-(5).
The fitting parameters are listed in Table 3 when R-squared is bigger than 0.98. Whichever line type was chosen, the slope of the displacement of No. 3 was either 200 times of the No. 1, or 10 times of No. 2, before the force reached the setting values. Therefore, the actual press speed may be a function depending on the set values of the pressing force, speed rate, press slope, and molding temperature. Then, there are two phenomena that need to be explained. The first one is why the pressing force was still rising, when the deformation rate decreased from 81 s to 200 s. The second one is why the relationships of the deformation velocities among the three simulations were not the same as the relationships among the processing parameters.

Viscoelasticity of Glass
When applying specified stress on glass at the molding temperature, strain (ε(t)) goes on increasing with time (t), which is called creep. In creep, the deformation of glass is a superposition of instantaneous elastic deformation, delayed elastic deformation, and viscous flow, which can be described by viscoelastic mechanical models. When applying a constant stress σ 0 at t = 0, the strain-time relationship can be obtained by Equations (6) and (7) according to Kelvin model [27].
E and η are elastic modulus and viscosities of the viscoelastic material, respectively. τ is defined as the term "retardation time" to describe the creep rate in the viscoelastic deformation.

Deformation Analysis with FEM and Discussion
The deformation of the glass includes two parts, which are the deformation of the upper half part and the lower half part of the glass preform. The respective deformation speeds of the upper half and the lower part are not known from the monitoring values in the experiment. FEM could be helpful for analyzing the deformation. The boundary condition could be displacement control and the element surface pressure control. The simulation model is firstly established, and then two groups of simulations with the two kinds of boundary conditions are carried out.

Simulation Model
To gain more information about the simulation result of the molded lens, threedimensional simulation was carried out. If the geometries of the molds were the same as the actual molds, the computation workload would be too large. Thus, the profile of the mold core was the same as the design, but the height of the mold core was smaller than the actual mold core. As the simulation in this paper only included the molding stage of the PGM process, which does not involve heat transfer, the simplified geometry model would be acceptable.
The three-dimensional geometric model is established in the finite element analysis software MSC.Marc (2016.0.0, MSC Software Corp., Newport Beach, CA, USA), as shown in Figure 5. There were three contact bodies including the upper and lower molds and the glass preform. The geometries of the upper and lower molds were artificially and subjectively divided into 20,335 and 20,638 elements, which include pentahedrons and hexahedrons. The geometry of the glass preform was divided into 66,206 elements, which were tetrahedrons. The use of these finite elements can meet the requirements of geometric model construction.
PGM process, which does not involve heat transfer, the simplified geometry model would be acceptable.
The three-dimensional geometric model is established in the finite element analysis software MSC.Marc (2016.0.0, MSC Software Corp., Newport Beach, CA, USA), as shown in Figure 5. There were three contact bodies including the upper and lower molds and the glass preform. The geometries of the upper and lower molds were artificially and subjectively divided into 20,335 and 20,638 elements, which include pentahedrons and hexahedrons. The geometry of the glass preform was divided into 66,206 elements, which were tetrahedrons. The use of these finite elements can meet the requirements of geometric model construction. The PGM machine instruction book shows that the pressing force that is detected by the load cell is displayed in real-time. However, when the actual pressing forces that were applied to the glass were the same, the theoretical deformations should be the same, which were different from the actual deformations shown in Figure 4. Thus, the actual pressing force should be a function dependent on other relative setting processing parameters. Therefore, the recorded forces should not be the boundary conditions. The pressures applied to the mold core in the simulations were separately designed for qualitative analysis.
As the values of the boundary conditions were not the processing parameters of the experiments, the viscoelastic parameters were used for calculating the simulations to describe the glass deformation over time. The use of the viscoelastic parameters of different glasses at the temperatures which make the glasses have a similar viscosity was accepted in this paper. The viscoelastic parameters of glass D-ZK3L at 550 °C, which can be a molding temperature, were used for the finite element simulations. The viscoelastic properties of glass D-ZK3L, other thermal properties of glass IRG206, and mechanical properties of the mold were listed in Table 4 [28][29][30]. The diagram of the established finite element model with boundary conditions is shown in Figure 6. The boundary conditions were applied to the end faces of the mold cores. The PGM machine instruction book shows that the pressing force that is detected by the load cell is displayed in real-time. However, when the actual pressing forces that were applied to the glass were the same, the theoretical deformations should be the same, which were different from the actual deformations shown in Figure 4. Thus, the actual pressing force should be a function dependent on other relative setting processing parameters. Therefore, the recorded forces should not be the boundary conditions. The pressures applied to the mold core in the simulations were separately designed for qualitative analysis. As the values of the boundary conditions were not the processing parameters of the experiments, the viscoelastic parameters were used for calculating the simulations to describe the glass deformation over time. The use of the viscoelastic parameters of different glasses at the temperatures which make the glasses have a similar viscosity was accepted in this paper. The viscoelastic parameters of glass D-ZK3L at 550 • C, which can be a molding temperature, were used for the finite element simulations. The viscoelastic properties of glass D-ZK3L, other thermal properties of glass IRG206, and mechanical properties of the mold were listed in Table 4 [28][29][30]. The diagram of the established finite element model with boundary conditions is shown in Figure 6. The boundary conditions were applied to the end faces of the mold cores.    According to the simulation results, the displacement in y direction of any node could be obtained. The displacement data of the three nodes were extracted: the center node of the top surface of the glass preform (Node_T in Figure 5), the center node of the bottom surface of the glass preform (Node_B in Figure 5), and the edge intermediate node of the glass preform (Node_M in Figure 5). Therefore, the total deformation of the glass preform can be calculated by subtracting the displacement of Node_B from the displacement of Node_T. The deformation of the upper part of the glass preform can be calculated by subtracting the displacement of Node_M from the displacement of Node_T. Therefore, the deformation of the lower part of the glass preform can be calculated by subtracting the total deformation from the deformation of the upper part of the glass preform.

Simulation with the Boundary Condition of Displacement Control
To analyze the different deformation speeds of the upper and lower parts of the glass when the total deformation is linear, the displacement of the lower mold is the boundary condition. A group of simulations were carried out. The deformations of the two half parts of the glass preform were analyzed.
A table that showed the data on displacement dependent with time was established in the software and was quoted in the process of the calculation in each simulation. The displacements of the mold were set as 0.013, 0.026, and 0.052 mm/s in the three simulations. These simulations were named 1X, 2X, and 3X. The deformation amounts of the upper half part and the lower half part of the glass and the shape of glass in the 2X experiment are shown in Figure 7a. The deformation amounts of the three simulations are plotted as three group curves as shown in Figure 7b. The three group curves also were processed through regression analysis with IBM SPSS Statistics software. The fitting parameters are listed in Table 5 when R-squared is bigger than 0.95. The line types of the curves representing the deformation of the upper and lower half part of the glass preform were highly linear. The slopes of the curves represent that the deformations of the upper and lower half part of the glass preform in the three simulations nearly had the same multiple relationships of the three boundary conditions. When the pressure applied to the glass was enough, the deformation velocity of the upper half part and the lower half part should be constant. The sum of the two constants should be the slope of the displacement boundary condition. However, the specific value should be related to the viscosity of the glass and the curvature of the molded glass surface.  The three group curves also were processed through regression analysis with IBM SPSS Statistics software. The fitting parameters are listed in Table 5 when R-squared is bigger than 0.95. The line types of the curves representing the deformation of the upper and lower half part of the glass preform were highly linear. The slopes of the curves represent that the deformations of the upper and lower half part of the glass preform in the three simulations nearly had the same multiple relationships of the three boundary conditions. When the pressure applied to the glass was enough, the deformation velocity of the upper half part and the lower half part should be constant. The sum of the two constants should be the slope of the displacement boundary condition. However, the specific value should be related to the viscosity of the glass and the curvature of the molded glass surface. The displacement of the lower axis is highly linear in the first 81 s of the forming stage in the PGM experiments, according to Section 4.2. Therefore, the deformation of the upper half part and the lower half part of the glass was linear. For the No. 1 experiment, the deformation velocity of the upper part of the glass was assumed as V α1 , and the deformation velocity of the lower part of the glass was assumed as V α2 .
The PGM machine instruction book pointed that the pressing force that was detected by the load cell was displayed in real-time. However, when the actual press forces that were applied to the glass were the same, the deformation should be the same. Thus, the actual press force is a function dependent on the time (t), and the actual press slope of the No. 1 experiment is assumed to be S 1 . When a force is applied to the glass preform with the spherical surface with a curve of 20 mm, the varying contact surface will lead to varying pressure.
In the progress of precision glass molding, the diffractive surface is upward and the aspherical surface is downward. The upper contact surface area, A(t), of the glass preform and the mold core can be calculated by the geometry relationship as shown in Figure 8a and can be expressed by Equation (8). The rate of change of contact area is assumed to be a(t), which can be expressed by Equation (9) when V α1 is smaller than 0.02 mm/s. The actual pressure σ α1 (t) applied to the upper surface of the glass can be expressed as Equation (10). The projection radius, Rp(t), of the lower contact surface shown in Figure 8b can be expressed by Equation (11). The normal force, Fn(t), applied to the lower contact surface can be expressed by Equation (12). The normal pressure applied to the lower surface is σα2(t), expressed as Equation (13).
The strain, εα1(t), of the upper surface can be concluded as Equation (14), according The projection radius, R p (t), of the lower contact surface shown in Figure 8b can be expressed by Equation (11). The normal force, F n (t), applied to the lower contact surface can be expressed by Equation (12). The normal pressure applied to the lower surface is σ α2 (t), expressed as Equation (13).
F n (t) = 2 9 2 V α2 t 0 The strain, ε α1 (t), of the upper surface can be concluded as Equation (14), according to Kelvin model. The strain of the lower surface is much smaller than the upper surface from Figure 7b.
where α is a calculated constant. To know the property of the curve, like in Equation (6), when τ is assumed to be 81 s, the curves when the parameters are σ 0 /E, 0.95 × σ 0 /E, 0.6 × σ 0 /E are plotted as in Figure 9. It was seen that the strain was highly linear from 0 s to 81 s, while the strain speed obviously decreased after 81 s. Therefore, the first phenomenon in the problem description part can be explained. The varying press force and the varying contact area led to a constant pressure before the 200 s mark in the three experiments. The strain-time curve for the viscoelasticity under constant pressure has similar properties of the curves in Figure 9.

Simulation with the Boundary Condition of Pressure Control
To obtain the deformation degree when the pressing force was different, a group of simulations were carried out. As the area of the bottom surface of the lower mold core was constant, the boundary condition in this section was the element surface distribution pressure. The pressure values that were applied to the element surface varied dependent on the time, t. The pressures in the No.1, No.2, and No.3 simulations were 0.1 t MPa, 0.2 t MPa, and 0.3 t MPa, respectively. The personal computer parameters and the simulation operation times were listed in Table 6. Since it was a three-dimensional simulation, it took nearly two hours to calculate each simulation. The deformations of the glass preforms in the three simulations are plotted as curves, as shown in Figure 10. The deformations of the lower half part of the glass preforms were nearly the same, but there were many big differences in the deformation of the upper half part of the glass preforms. The deformations of the upper half part of the glass preforms were nearly linear. The deformation rates were about 0.025, 0.04, and 0.051 mm/s in the three simulations, according to IBM SPSS Statistics software. The relationships of the deformation rates were not the same as those of the applied pressures in the simulations. Table 6. Personal computer parameters and the simulation operation time.

Properties
Informations Central Processing Unit Intel(R) Core(TM) i7-5820K CPU @ 3.30GHz Graphic Processing Unit NVIDIA GeForce GTX 960 Random Access Memory 32.0 GB

Simulation with the Boundary Condition of Pressure Control
To obtain the deformation degree when the pressing force was different, a group of simulations were carried out. As the area of the bottom surface of the lower mold core was constant, the boundary condition in this section was the element surface distribution pressure. The pressure values that were applied to the element surface varied dependent on the time, t. The pressures in the No. 1, No. 2, and No. 3 simulations were 0.1 t MPa, 0.2 t MPa, and 0.3 t MPa, respectively. The personal computer parameters and the simulation operation times were listed in Table 6. Since it was a three-dimensional simulation, it took nearly two hours to calculate each simulation. The deformations of the glass preforms in the three simulations are plotted as curves, as shown in Figure 10. The deformations of the lower half part of the glass preforms were nearly the same, but there were many big differences in the deformation of the upper half part of the glass preforms. The deformations of the upper half part of the glass preforms were nearly linear. The deformation rates were about 0.025, 0.04, and 0.051 mm/s in the three simulations, according to IBM SPSS Statistics software. The relationships of the deformation rates were not the same as those of the applied pressures in the simulations.  The actual press force should be a function with the setting press speed. According to Section 4.2, the actual press slope of the No.1 experiment was assumed to be S1. Therefore, the actual press slopes of the No.2 and No.3 experiments were assumed to be S2 and S3. Since the upper deformation was similarly linear, the strain calculation was similar to that in Section 4.2. Therefore, the rate of the strain Vα 1(t) of the upper half part of the glass in the first simulation can be calculated as Equations (15) and (16). Then, the ratio of the deformation velocities of the upper part of the glass in No.1 experiment Vα 1(t), No.2 experiment Vβ 1(t), and No.3 experiment Vγ 1 (t) can be calculated by Equations (17) and (18). The actual simulation results are shown as Equation (19). The error between the model result and the simulation result was about 0.116 for the second simulation, while it was about 0.151 for the third simulation. The error may be caused by the approximation of the contact area in Equation (9). However, Equation (17) plays a certain role in setting the manufacturing parameters.    The actual press force should be a function with the setting press speed. According to Section 4.2, the actual press slope of the No. 1 experiment was assumed to be S 1 . Therefore, the actual press slopes of the No. 2 and No. 3 experiments were assumed to be S 2 and S 3 . Since the upper deformation was similarly linear, the strain calculation was similar to that in Section 4.2. Therefore, the rate of the strain V α1 (t) of the upper half part of the glass in the first simulation can be calculated as Equations (15) and (16). Then, the ratio of the deformation velocities of the upper part of the glass in No. 1 experiment V α1 (t), No. 2 experiment V β1 (t), and No. 3 experiment V γ1 (t) can be calculated by Equations (17) and (18). The actual simulation results are shown as Equation (19). The error between the model result and the simulation result was about 0.116 for the second simulation, while it was about 0.151 for the third simulation. The error may be caused by the approximation of the contact area in Equation (9). However, Equation (17) plays a certain role in setting the manufacturing parameters.
V α1_simulation , V β1_simulation , and V γ1_simulation are the deformation velocities of the upper part of the glass in No. 1, No. 2, and No. 3 simulations, respectively. Y(t) is a function that changes over time.
The ratio of the deformation rates of the glass preforms in the experiments was 0.007:0.012759:0.017852. Thus, the ratio of the rates of the applied force in the experiments was about 0.007 2 :0.012759 2 :0.017852 2 , immediately 1:2.92:5.72. The ratio does not have a linear relationship with the pressing speed. Therefore, the actual press force in the PGM experiments needs to be further studied.

Conclusions
At present, there is little discussion on the deformation of glass preforms under multi-machining parameters of PGM technology, but more discussion on the properties of the molded lens. This paper investigated the glass preform deformation under multimachining parameters when fabricating ChG diffractive lens with PGM technology through experiments, FEM simulations, and theoretical analysis. The primary findings are as follows: 1.
In different experiments, the ratio of the deformations of the upper half part and the lower half part of the glass preform are approximately a constant when the pressing forces are big enough, but the pressing speeds are set differently. When the pressing force is set linearly over time and the primary pressing force is too small to achieve the setting pressing speed, after a period of time, the deformation speed of the glass preform starts to decrease significantly.

2.
FEM can be used to quickly establish the glass preform deformation model in the time when the glass is in a viscoelastic state. The error of the obtained glass preform deformation model is less than 0.16.

3.
The established model can be used as a reference to efficiently optimize PGM processing parameters when the designed lens has two different base radii of curvature. Data Availability Statement: Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.