Real-Time Thermomechanical Modeling of PV Cell Fabrication via a POD-Trained RBF Interpolation Network

This paper presents a numerical reduced order model framework to simulate the physics of the thermomechanical processes that occur during c-Si photovoltaic (PV) cell fabrication. A response surface based on a radial basis function (RBF) interpolation network trained by a Proper Orthogonal Decomposition (POD) of the solution fields is developed for fast and accurate approximations of thermal loading conditions on PV cells during the fabrication processes. The outcome is a stand-alone computational tool that provides, in real time, the quantitative and qualitative thermomechanical response as a function of user-controlled input parameters for fabrication processes with the precision of 3D finite element analysis (FEA). This tool provides an efficient and effective avenue for design and optimization as well as for failure prediction of PV cells.


Introduction
As Silicon wafer thickness is expected to decrease in the future [Metz, Fischer and Trube (2017)], it becomes critical to estimate and predict the thermomechanical response associated with each critical cell and module due to different thermal fabrication steps such as emitter diffusion, contact firing, and stringing/tabbing. Thermal shock, materials cooling with a different thermal coefficient of expansion, and other effects can result in crack initiation and propagations. To this end, we developed a response framework based on a radial-basis function (RBF) interpolation trained by a Proper Orthogonal Decomposition (POD) capable of accurately and virtually instantaneously predicting the physics of the thermomechanical processes involved in PV cell fabrication. This framework provides the capability of approximating thermomechanical stresses, strains, and deformations resulting from the thermal processes during c-Si PV solar cell and module fabrication at the same speed or faster than the very process that it is simulating, hence rendering a real-time modeling tool. This while taking advantage of the detail and accuracy of grid-converged Finite Element Analysis (FEA).
The key to achieving this goal is to generate, beforehand and off-line, an extensive set of solutions that compute such engineering design endpoints such as principal stresses, principal strains and deformation components over the PV cell surface. These set of solutions have been generated using FEA simulations within a predefined parameterized design space (various thickness of silicon, aluminum, and silver, number of busbars and maximum firing temperature). These solutions are then organized to form the basis snapshots of a POD matrix. The POD snapshot matrix is in turn used to construct an orthonormal basis of the mechanical behavior of the silicon wafers during the manufacturing process. Solutions derived from the POD basis are subsequently used in an interpolation network constructed using radial basis functions (RBFs) to predict a real-time response of the field solution for a given set of input parameters. The entire POD-RBF interpolation network basis data for interpolation is stored in a database that can be accessed instantly and remotely and therefore predict stresses, strains, and deformations virtually instantaneously. This particular efficiency advantage makes this technique applicable in the prediction of the potential formation and propagation of cracks during key thermal fabrication steps which plays a key enabling role in the PV cell design process. The POD-RBF interpolated network acts as a multidimensional interpolation that preserves the physics of the problem while providing accurate and fast predictive capabilities in the course of the design process without having to resort to multiple and computationally burdensome FEA models. This POD-RBF interpolated network is validated against the FEA field solution for several test cases. Although FEA is an accurate and established numerical analysis tool, one of its significant drawbacks lies in its need for a well-defined mesh that often requires user-interaction in order to yield accurate results. Therefore, most of the computational time required in FEA is attributed to the construction of an appropriate mesh capable of correctly capturing the problem physics. In addition, these appropriately constructed meshes for FEA preclude the ability to be automatically generated or adapted to different problem conditions and hence posing a major challenge when requiring that FEA acts as a tool for design and optimization. Therefore, the development of a fully automated, accurate, and stable numerical thermomechanical analysis tool is of paramount importance for an eventual design and optimization suite. To this end, a stress, strain and deformation numerical prediction tool fully-automated for design and optimization is developed based on a POD-RBF interpolation network. This tool is capable of yielding real-time predictions while preserving the accuracy of the underlying FEA without the need for user-interaction and its large computational demands. This efficiency advantage enables the tool to be used as a predictive tool during the PV manufacturing process and as a solver within a design and optimization framework.
The foundations of POD were introduced at the beginning of the 20 th century as a statistical tool [Pearson (1901)] with the objective of finding an optimal basis for least-square curve fitting. The technique has since been renamed and reintroduced in the literature for a large number of different applications. Additional names and variations of the POD method include: Karhunen-Loeve decomposition, principal component analysis or singular value decomposition [Karhunen (1947);Feeny and Kappagantu (1998)]. In addition, the POD method has been applied in a variety of problems including signal processing and control theory, human face recognition, data compression, fluid mechanics, parameter estimation and others. The POD can be employed to reproduce a reduced-order approximation of the solution field with very high fidelity as it is capable of capturing the dominant or principal components of the data by only exploring a few modes. This is possible due to the capability of the POD to yield the optimal basis for least-squares approximation by constructing a set of vectors in a rotated frame of reference, where the angles of rotation represent the POD basis [Bialecki, Kassab and Fic (2005); Fic, Bialecki and Kassab (2005)]. Radial basis function network can be used as a substitutive solution to the conventional neural network in various applications of signal processing and system control [Chen, Cowan and Grant (1991)]. Also, radial basis function-based neural network has been used for developing a nonlinear system model [Han, Guo and Qiao (2018)]. Also, recent work has been done on developing a power forecasting model of photovoltaic cell using a radial basis function neural network by Xu et al. [Xu, Chen, Zhou et al. (2019)]. The combination of the POD method with RBFs offers the ability to parametrize problems in terms of design parameters that can be varied to produce different columns of a data matrix that can later be interpolated and thereby mitigating the onerous task of having to solve the forward problem every time the parameters are varied. The POD-RBF framework is therefore capable of preserving the correlation between the forward problem and the design parameters to seek a solution of an ill-posed problem Kassab (2005, 2008)]. Furthermore, the POD-RBF framework is ideally suited for addressing design and optimization problems, as it not only provides robustness, accuracy, and stability, but provides order reduction, noise filtration, and regularization. Different application of the POD-RBF framework for design and optimization and parameter estimation problems can be found in the literature showing significant computational gains [Klimanek, Bialecki and Ostrowski (2010)] as well as adaptability [Rogers, Kassab, Divo et al. (2012)] and the capability of coupling with CFD solvers [Huayamave, Ceballos, Barriento et al. (2017)].

Material properties
When the silicon wafer enters the furnace belt as seen in Fig. 1(a), the firing process goes through four thermal steps. The first step is the initial temperature ramp up where the paste solvents in the wafer are volatilized, the second step is the burn out which removes all of the organic binder that was used in paste formation. The third step is the sintering which the silver metal forms a bond with the underlying silicon substrate to form metal contact. The final step is the wafer cool down phase. In the contact firing process as shown in Fig.  1(a), the initial solar cell configuration consists of several layers as shown in Fig. 1 Since the silicon nitride [Burkhardt and Marvel (1969); ; Shimada, Matsushita, Kuratani et al. (1984); Zhang and Grigoropoulos (1995)] and the silicon dioxide [Cahill (1990); Schafft, Seuhle and Mirel (1989); Tada, Kumpel, Lathrop et al. (2000)] layers have dimensions in the nanoscale and the other layers have dimensions in the microscale, the thermal effects of these two layers are neglected and therefore not included in the analysis configuration. In addition, temperature dependent material properties are necessary for the current analysis of thermal processes due to the wide range of temperatures at which the silicon wafers are exposed. The relevant thermal and thermoelastic properties for each of the layer i.e., silicon, aluminum, and silver are obtained from the literature survey.

Mechanical properties
Silicon is an anisotropic material and the complete description of its elasticity is thus given by a fourth order rank elasticity tensor, which is mathematically cumbersome to manipulate when a rotation in the orientation of the material is necessary. Fortunately, the silicon lattice has a property called cubic symmetry, in which all directions and planes rotated 90º are equivalent. This allows for an expression of the elastic properties of silicon in terms of orthotropic material constants. Such expression is given by Hopcroft et al. [Hopcroft, Nix and Kenny (2010)] as follows for a frame of reference for a standard silicon crystal in the (100) direction: With the following formulation for material elasticity, this is straightforward to input in commercial FEA software:

Thermomechanical FEM implementation
To accurately simulate the contact firing process as seen in Fig. 1(a), a single solar wafer is isolated, and a finite element analysis is performed. The wafer is composed of several layers as shown in Fig. 1(b). For our analysis, thermal effects of silicon nitride and silicon dioxide layers are omitted from the analysis since these two layers have nanoscale dimensions while other layers have microscale dimensions. Consequently, a composite model is created using a silicon layer, an aluminum layer, silver busbars layers and silver pads layers with the wafer frontal dimensions of 156 mm×156 mm as shown in Fig. 2.

Figure 4: Wafer boundary conditions
The analysis is performed utilizing a typical firing temperature profile shown in Fig. 5. The region of interest is the temperature profile between 40 seconds and 58 seconds when the wafer enters and exits the conveyer. Therefore, transient non-linear analysis is performed for 18 seconds.

Mesh domain and optimization
The layers are modeled as a composite using shell elements due to the length-to-thickness ratio. All the components are meshed using 4-node thermally coupled doubly curved shell elements for reduced integration, and mesh convergence analysis is implemented for the 3, 4, and 5 busbars geometry to obtain an accurate solution. After mesh convergence analysis, the final mesh of the 3 busbars model contained 19,148 shell elements, the final mesh of the 4 busbars model contained 23,622 shell elements, and the final mesh of the 5 busbars model contained 27,896 shell elements. Fig. 6 shows the final mesh for 3 and 5 busbars wafers. The field variables were collected from equally distributed locations using a script in the FEM solver (ABAQUS). This script was written to automatically query and store the values of the displacements (δx, δy, δz), principal stresses (σ1, σ2, σ3) and principal strains (ε1, ε2, ε3) at 2,809 (53×53) predefined and uniformly distributed locations throughout the domain interpolated from the non-uniform FEM mesh.

Establishing the POD-trained RBF interpolation network
The POD-RBF interpolation network is built in two stages:

Produce the POD expansion
First, the data snapshot vector under a given set of sampled design parameters is created. Specifically, the data snapshot vector is the collection of N sample values of u, which is the field variable under consideration. In the firing process, N ,which is obtained from the discretization of the FEM, is equal to 2,809 (53×53) equally distributed locations on the surface of the wafer where the field variable is sampled. The field variables under consideration are the discrete values of the three principal stresses, the three principal strains, and the three components of deformation, for a total of nine (9) field variables to track during the fabrication process. A collection of M snapshots denoted as u j (for j=1, 2 … M) is generated by altering the values of the design parameters according to those specified in Tab. 4. These design parameters are denoted as par=(st, at, bt, nb, ft), where st: principal stresses, at: aluminum thickness, bt: silver busbar/pad thickness, nb: number of busbars, and ft: maximum firing temperature. These parameters are kept in a [L×1] vector par where L is the number of design parameters upon which the field variables depend, i.e., L=5 for the contact firing process. Each u j is then stored in a column of a rectangular [N×M] snapshot (or data) matrix U.
The main objective of the POD is to establish a set of orthonormal vectors Φ j (for j=1, 2 … M) that resemble the snapshot matrix U in an optimal way. The matrix Φ is commonly referred to as the POD basis and is given by: Here V represents the eigenvectors of the covariance matrix C, defined as C=U T U. These eigenvectors can be easily determined through the nontrivial solution of the general [M×M] eigenvalue problem: Here Λ is a diagonal matrix that stores the eigenvalues λ of the covariance matrix C.
Because the covariance matrix C is the product of the transpose of the snapshot matrix U with itself, C is consequently symmetric and positive definite, and therefore all eigenvalues λ are real and positive. In addition, the eigenvalue decomposition routine can be configured to sort the eigenvalues λ in descending order; yielding a structure that is analogous to the energy modes of the system (from high to low). This energy decreases very rapidly with the increasing mode number and because higher modes hold very little energy (or data) of the system, they can be neglected without significantly influencing the accuracy of the representation while simultaneously filtering unwanted noise from the system. This process is known as the truncation of the POD basis, and it is accomplished by deciding which modes of energy of the system can be neglected based on the residual yielded in the subsequent truncated interpolation. This truncation also helps to smooth out the numerical noise in the model which is critical for an accurate and stable response surface interpolation. The resulting POD basis � , referred to as the truncated POD basis, consists of K<M vectors as: � = � (7) This truncation also corresponds to the truncation of the corresponding eigenvectors in the eigenvector matrix, leading to a truncated eigenvector matrix denoted as � , that stores the first K eigenvectors of the covariance matrix C. In this study, the value of K was chosen to be 10, which was found to produce a smooth interpolation without significant loss of accuracy as evidenced in the results. The truncated POD basis is also orthogonal and therefore satisfying � � = and yields optimal approximation properties in a reduced order. Once the truncated POD basis � is known, the snapshot matrix U can be reconstructed by: Here A represents the amplitudes associated with the snapshots u j . Now, recasting the orthogonality of the truncated POD basis � , these amplitudes A can be simply determined as = � .

Generate RBF interpolation framework from the POD expansion
In this stage RBFs are employed to interpolate the field variables as a function of the design parameters within the design space. RBFs then use the truncated POD expansion to formally establish an interpolation framework capable of reproducing the field as a function of any desired value of the design parameters. The procedure is developed herein in general assuming that the vector par may contains as many design parameters as employed in the POD expansion process (size L).
To illustrate this, consider the inverse Hardy Multiquadrics RBF as described in Hardy [Hardy (1971, 1990)]: Here c denotes the RBF shape factor which controls the steepness of the function, par corresponds to a vector containing arbitrary values of the design parameters, par i correspond to the set of design parameter values used to generate the snapshots u i for i = 1,2 … M. The function fi(par) therefore provides a continuous algebraic representation that depends uniquely upon the deviation from pre-established poles (par i ). In standard applications of RBF, par corresponds to a space coordinate (x,y,z) and thus � − � corresponds to the Euclidean (radial) distance between (x,y,z) and the poles (xi,yi,zi), and therefore the name: radial basis function. An explicit approximation of the dependence of a [M×1] column vector of the snapshot matrix U on parameters par can be regenerated as: At this stage a collocation process is formulated by requiring that the approximation ( ) retrieves the snapshot u i at the sampled parameter values par i for i = 1,2 … M: or, in matrix form: Using the orthogonality property of � , the matrix of coefficients B can be evaluated by simple inversion as: Having evaluated the coefficient matrix B, the POD-trained RBF interpolation given by Eq. (10) can now be used to reproduce the unknown field in the form of an [M×1] vector of field values at M discrete locations that corresponds to any arbitrary set of parameters par as: This can be thought of as a numerical Eigenfunction expansion of the solution or as a multidimensional response surface of the field as a function of design parameters. The POD expansion and the RBF interpolation were first coded on MathCAD (v15) to test the validity and accuracy of the approach for selected cases. Then a Fortran code was developed for the POD-RBF interpolation network to automatically read the FEM data from the ABAQUS script and output the data for post-processing.

Thermomechanical FEM results
All the FEM runs corresponding to the 243 combinations of design parameters were successfully performed and the field results of principal stresses, principal strains, and deformation components were used to build the POD-RBF interpolation network. The von Mises stress and out-of-plane deformation contour plots are displayed below for one of these FEM cases with a configuration of 3 busbars/pads, silicon substrate with thickness of 100 µm, aluminum substrate with thickness of 50 µm, busbars/pads with thickness of 50 µm, and a maximum firing temperature of 1073 K. The contour plots, as shown in Figs. 7 and 8, are generated at the bottom surface of the silicon substrate where the stresses are maximum due to the deflection (bowing) direction.

Experimental validation
To validate the computational analysis, experimental results were obtained during the contact firing stage of a solar cell manufacturing process provided by Suniva, Inc. Four groups of 15 silicon solar wafers (156 mm pseudosquares, cut from 200 mm rounds) were used for this experiment. They were etched in a KOH bath prior to subsequent printing and firing. The samples were printed with silver and aluminum pastes and fired through a belt furnace. Not all samples survived processing without breaking. The broken wafers were excluded from subsequent measurements. For this validation, four cases outside the parametrized space were considered and are summarized in Tab. 6. In addition, the substrate aluminum thickness and the silver busbars/pads thickness were kept constant at 50 µm.   The four cases were successfully performed and the errors between the out-of-plane deformations obtained by the experimental and the FEM results are shown in Tab. 7. The FEM von Misses Stress and deformations results for these cases are shown in Figs. 9 and 10. In addition, wafer and paste thicknesses were measured by taking the mass of the wafers during processing and measuring the feature sizes on the surface. Coupling the area and the material density allowed us to calculate the thickness of the various layers. The thickness was confirmed by independent measurements using a micrometer (wafer thickness) or an optical microscope (paste thickness) although these methods have larger errors and do not give the average over the entire sample. The bowing (deformation) of the wafers was measured by eye with a standard ruler relative to the flat surface (lab bench) to within an estimated accuracty of ±0.2 mm.

POD-RBF interpolation results
The result of the POD-RBF interpolation network is validated using a set of FEM thermoelastic solutions. After the nine (9) snapshot matrices U were formed for each of the nine (9) field variables (σ1, σ2, σ3, ε1, ε2, ε3, δx, δy, δz), the decomposition is performed and tested. The resulting nine (9) covariance matrices C each with dimension [243×243] are formed as C=U T U. After the eigenvalue decomposition of these covariance matrices, the resulting 243 eigenvalues λ range from a maximum value of about 10 9 to a minimum value of about 10 -5 . More importantly, the eigenvalues decrease very rapidly from the largest value to less than 10 3 after the first 8 eigenvalues, indicating that most of the system information (energy) is contained and can be extracted from the first few eigenvalues using a truncated POD basis. For silicon thickness (st) of 100 µm, aluminum thickness (at) of 10 µm, busbar and contact pad thicknesses (bt) of 30 µm, number of busbars (nb) of 3, and maximum firing temperature (ft) of 873 K, Figs. 11 to 13 below show contour plots of the first principal stress, first principal strain, and out-of-plane deformation for both the FEM solution (left) and the POD-RBF interpolation (right). The qualitative comparison of the two solutions shows virtually no error while a quantitative comparison through relative RMS errors reveals a difference of 0.05% for stress and 0.01% for strain and deformation. The developed real-time response framework based on the POD-RBF interpolation scheme has the capability of efficiently generating the principal stresses, principal strains, and deformation components for any arbitrary set of design parameters within the design space.
In fact, this can be tested even in the extreme case where a user will erroneously input a fractional value for an expected integer quantity such as the number of busbars. For instance, a case can be tested where the input parameters are such that st=135 µm, at=25 µm, bt=45 µm, nb=4.5, and ft=1000 K. Notice that the number of busbars input parameter is 4.5 representing a non-physical quantity. Despite this, the POD-RBF interpolation network is robust enough to make the necessary corrections to depict the first principal stress, first principal strain, and out-of-plane deformation component fields in such a way as to predict sort of a "ghost" busbars appearing on the surface of the PV cell as shown in Fig. 14. Of course, this is just a demonstration of the capabilities of the POD-RBF interpolation network for case that has no physical meaning where in practice the input parameters are simply validated to fall within the design envelope and represent physicallycorrect cases (such as an integer number of busbars).

Conclusion
A real-time response framework based on a POD-trained RBF interpolation network is developed to accurately predict the thermomechanical response that results from the variable conditions occurring during PV cell fabrication processes. The POD-RBF interpolation network acts as a multifaceted algebraic response surface that can produce instantaneous stress, strain, and deformation fields as a function of input parameters within a design envelope. In this paper, this approach was demonstrated for the contact firing process of PV cell fabrication where five input parameters were prescribed; namely the thickness of the silicon substrate, the thickness of the aluminum substrate, the thickness of the silver busbars, the number of busbars, and the final firing temperature. These design parameters were swept within a design space consisting of three values for each parameter leading to a total of 3 5 =243 total FEM cases to build the basis for the interpolation network. After the POD-RBF interpolation network was built, it was tested against FEM results outside the parameterized data showing excellent agreement with RMS errors of less than 1%. This tool not only provides real-time access to accurate thermomechanical solution fields but, given the appropriate models and constraints, it also provides an effective avenue for design and optimization as well as for failure prediction of PV cells.