Machine Learning Inversion Method for Elastoplastic Constitutive Parameters of Encapsulation Materials

TL;DR

This paper introduces a machine learning method to accurately measure material properties of encapsulation materials used in 3D electronics.

Contribution

A novel neural network-based inversion method is proposed for high-accuracy elastoplastic parameter measurement using nanoindentation.

Findings

The relative error of material parameters is less than 3% with a 95% confidence interval.

The inversion convergence error for indentation response parameters is less than 0.1%.

Finite element simulations verified the method's impact on process stress in TCV products.

Abstract

Accurate measurement of material mechanics parameters is crucial for evaluating process quality and product reliability and is a major challenge in the development of 3D heterogeneous integration technology. Aiming to perform high-accuracy measurements of the elastoplastic nonlinear constitutive parameters of microelectronic materials using the nanoindentation testing technique, we take advantage of a neural network to construct a forward characterization model to characterize these response characteristic parameters for different materials, design an improved algorithm for obtaining a reverse iterative solution of the forward characterization model, and develop a material mechanics parameter measurement method to solve overdetermined equations using the least-squares method. This method was further improved by addressing the issues of algorithm stability and solution uniqueness, achieving high-precision and fast reverse solutions for elastoplastic constitutive parameters. The relative error of the material parameters is less than 3% (95% confidence interval), the maximum error is less than 8%, and the inversion convergence error of the key indentation response characteristic parameters is less than 0.1%. The difference between the measured material parameters and the theoretical model in the influence on the process stress of TCV (through ceramic via) products is verified through finite element simulation.

Figures

Full text

1. Introduction

The mechanical properties of materials have a significant impact on the reliability of microelectronic integrated products. Three-dimensional heterogeneous integration (3D HI), a key technology for improving the performance of microelectronic products, needs to address the complex issue of compatibility between different malarial systems. The characteristics of microelectronic materials are closely coupled with process methods; therefore, accurate measurement of material mechanical parameters is crucial for evaluating process quality and product reliability and is a major challenge in the development of 3D HI technology. Traditional material mechanics testing methods require the preparation of macroscopic (cm level) test samples that meet the requirements of mechanical testing machines. However, this method is not feasible for the microelectronic materials formed during manufacturing processes, such as μm scale wiring and interconnection hole materials (usually gold or copper) formed during the electroplating process, in situ grown dielectric films [1], etc. The nanoindentation testing technique utilizes high-resolution sensing devices to continuously monitor payload (P) and displacement (h) during the indentation process and solves the mechanics parameters of the material based on the acquired P–h curve [2]. As nanoindentation testing only requires μm scale samples and supports near non-destructive in situ testing, it has gradually become a fast and simple mechanical performance testing method in the field of microelectronic packaging in recent years [3,4], widely used in the research of TSV (through silicon via) [5,6,7], TGV (through glass via) [8], MEMSs (micro-electro-mechanical systems) [9,10], nano metal bonding [11], etc.

Through the adoption of the widely accepted Oliver Pharr model [12] and the Cheng–Cheng method [13], nanoindentation testing can now be used to reliably acquire the elastic modulus and hardness values of materials [14]. To accurately design a method for determining the reliability of microelectronic products, it is necessary to obtain the key elastoplastic nonlinear constitutive parameters of materials, such as yield strength (σ_y_) and hardening index (n). However, there is a complex high-order nonlinear coupling relationship between these mechanical parameters and nanoindentation test data [15], which is difficult to calculate directly using analytical methods. Early research mainly relied on dimensional analysis [16,17,18], in which functions used to fit material parameters are constructed based on a certain amount of testing/simulation data. Since these methods are limited by data sources and the complexity of dimensionless functions, achieving high calculation accuracy is difficult, and the calculation accuracy varies significantly in different material parameter ranges. Another approach is to combine dimensional analysis with finite element method (FEM) for iterative calculation. This approach can improve calculation accuracy [11,19] but has the drawbacks of high computational cost and long calculation time, making it unsuitable for practical applications. In the bid to expand the application scope of the nanoindentation testing technique, explorations into high-precision, high-confidence, and fast nonlinear material parameter inversion methods are urgently needed.

Deep learning is a data-driven method for extracting patterns and features and has achieved great success in solving complex nonlinear problems such as image processing and natural language processing. Deep learning has unparalleled advantages in complex nonlinear modeling compared with traditional methods and is expected to provide a better route for the inversion of material nonlinear constitutive parameters in nanoindentation testing [20]. Marimuthu et al. [21] conducted regression modeling on P–h curve sampling data and material parameters and introduced a physical constraint equation of unloading initial slope (S) to construct an inversion method based on a physics-informed artificial neural network (PI-ANN). Park et al. [22] determined the minimum input feature parameters (h_max_, S, W_p_/W_t_, and C) and training data scale through correlation analysis. The application of transfer learning has also greatly reduced the scale of high-precision training dataset [23,24]. Jiao et al. [25] systematically studied the solution uniqueness problem in the reverse solution of elastoplastic constitutive parameters from the perspective of machine learning and proposed the use of pile-ups generated during the indentation process as the input feature of machine learning. Previous research has fully demonstrated the potential of deep learning in the inversion of nanoindentation data, but the high sensitivity of nanoindentation curves to material parameters and the solution uniqueness of inversion [25,26] make obtaining accurate and stable solutions using deep learning difficult, which thus necessitates in-depth research.

This study proposes a novel method for measuring material mechanics parameters based on the principle of solving overdetermined equations using the least-squares method. First, a forward characterization model of the response characteristic parameters of the nanoindentation process under different material mechanics parameters is constructed based on a neural network. Then, the least-squares inverse iteration is performed based on the forward characterization model to obtain the solution. Thus, the problem of measuring material parameters is transformed into the problem of finding the least-squares optimal solution of nonlinear overdetermined equations. This method is further improved by addressing the issues of algorithm stability and solution uniqueness. The method proposed in this paper avoids a major problem in the reverse solution of material parameters: the solution being difficult to obtain because the analytical relationship between material parameters and indentation experimental data is hard to establish. This paves the way for achieving rapid, accurate, and stable extraction of material mechanical parameters from nanoindentation experimental data. FEA simulation data demonstrates that the proposed algorithm can yield reverse solution results with high accuracy (both the strengthening index n and yield strength σ_y_ are <5%) and high stability and can avoid the multi-solution problem (multiple calculations yielding the same solution, and a convergence mean square error (MSE) < 1 × 10^−8^), which is superior to the results reported in References [11,19,23].

2. Methodology

2.1. Constitutive Model

This study focuses on a typical elastoplastic constitutive model featuring linear elasticity governed by Hooke’s Law and power law hardening behavior. This model corresponds to the elastoplastic constitutive equation presented in Equation (1), with its associated stress–strain curve illustrated in Figure 1. Given that most materials—particularly metals—adhere well to this model, the research is dedicated to inverting the power law constitutive parameters of materials using nanoindentation test data [27].

$$[eqn]$$

where E is Young’s modulus (Pa);

• R is the strength (Pa);
• σ_y_ is the yield strength (Pa);
• ε is the total strain;
• ε_p_ is the plastic strain;
• n is the strain strengthening index.

2.2. Reverse Analysis Framework

The overall framework of the material mechanics parameter reverse calculation method based on the principle of solving overdetermined equations using the least-squares method proposed in this article is shown in Figure 2. It mainly includes a neural network model constructed based on the basic material parameters and the characteristic parameters of the pressure–displacement curve, and a neural network reverse solution algorithm based on the improved LM (Levenberg–Marquardt) iteration algorithm.

In this study, finite element simulation was used to conduct a virtual experiment of nanoindentation, obtaining a large amount of high-precision data with a low cost for the construction of the neural network model. The simulation tool used in the virtual experiment was the static mechanical analysis module of the commercial software Ansys Workbench. The tetrahedral Berkovich indenter was substituted with a 2D axisymmetric model with the same projected area to accelerate the simulation speed [28]. The grid density in the area near the indenter and the pressed material was set to 50 nm, and the number of cells was 18,396. The mesh size in the farther regions was set to 150 nm (Figure 3). The mesh quality of the finite element simulations performed via Ansys in this study has undergone mesh independence analysis. This adequately ensures that the extracted simulation result data (e.g., S, C, and P_m_) are not affected by the model’s mesh division, with the actual discrepancy less than 0.1%. The minimum number of calculation steps during the loading and unloading stages was set to 1000. Considering that Poisson’s ratio (ν) has a small impact on the calculation results, it was set to a fixed value of 0.32. Five variables were set in the virtual experiment, and the range of values and factor levels were determined based on the possible elastoplastic material parameters of electroplated copper (Table 1). A total of 2016 rounds of a full-factor orthogonal DOE simulation experiment were conducted, each round yielding one piece of P–h data. Automatic creation of the simulation model and submission of the simulation results were achieved by writing Python (The version is Iron Python 2.7.11) scripts.

The indentation P–h curve data corresponding to each group of material parameters was obtained through a virtual simulation experiment (Figure 4). The curve feature parameters were defined and extracted to train the neural network model, including the following:

• (a)P_m_ is the maximum pressure, mN.
• (b)h_m_ is the maximum displacement, μm.
• (c)h_f_ is the residual displacement, μm.
• (d)C is the loading curvature, calculated using the equation $[eqn]$ .
• (e)S is the stiffness coefficient, defined and calculated using Equation (2) based on the unloading curve, where parameters B and m are obtained through fitting according to the unloading power function. The stiffness coefficient can also be obtained directly using a nanoindentation instrument; thus, m and B can be calculated directly using Equations (3) and (4):

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

• (f)W_t_ is the total indentation work, which is the area under the P–h curve when the indentation reaches its maximum displacement. It can be calculated using Equation (5):

$$[eqn]$$

• (g)W_u_ is the unloading work, which is the area under the unloading curve. It can be calculated using Equation (6):

$$[eqn]$$

• (h)W_l_ is the absorption work, $[eqn]$ .

2.3. Neural Network Structure and Training Algorithm

(1) The neural network used in this study is a fully connected BP neural network, which has the structure shown in Figure 5. It has an input dimension of 6 and an output dimension of 4. The network consists of two hidden layers: the first hidden layer is configured with 30 neurons and adopts the hyperbolic tangent (tanh) function as its activation function; the second hidden layer contains 20 neurons with the sigmoid function as its activation function. The output layer employs a linear activation function.

The first hidden layer of the neural network uses the tanh function (7), and the second hidden layer of the neural network uses the sigmoid function (8):

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ is the weighted sum (with added bias) of all outputs of the previous layer that is input into the neuron function of the l-th layer.

(2) Training algorithm with Bayesian regularization

To improve the generalization prediction ability of the model constructed based on the neural network proposed in this paper and to prevent overfitting, we adopted Bayesian regularization to train the algorithm. In addition, the L2 regularization term is constructed by adding the sum of squares of the model parameters based on the MSE loss function to limit the size of the model parameters, thereby limiting the number or magnitude of the model parameters and improving the generalization performance of the model [29,30]. Thus, the loss function E_ω_ used to train the neural network is as follows:

$$[eqn]$$

where λ is the regularization hyperparameter, which is set to 0.5 in this study.

For data partitioning, the training, validation, and test sets accounted for 70%, 15%, and 15% of the total dataset, respectively. During training, input features were linearly normalized to the range of [−1, 1] (the specific range was adjusted based on the sensitivity differences in each parameter; see Table 2 for details), while output targets were normalized to [0.1, 0.9]. The model adopted mean squared error (MSE) as the performance evaluation metric, with the optimization objective of minimizing the error. The training algorithm used was Bayesian regularization (trainbr), with the relevant hyperparameter settings as follows: maximum number of training epochs = 1000, minimum gradient threshold = 1 × 10^−7^, initial damping coefficient (Mu) = 0.005, damping coefficient decay ratio = 0.1, growth ratio = 10, and maximum damping coefficient = 1 × 10^10^.

2.4. Reverse Solution Algorithm

(1)Definition of the reverse solution model

The above neural network (as shown in Figure 5) can be written as a function mapping the representation of independent vs. dependent variables: $[eqn]$ , where

$$[eqn]$$ $$[eqn]$$

The equation is further divided into a system consisting of four equations. This system of nonlinear equations mapping the input and output variables is the aforementioned neural network model constructed in this study.

$$[eqn]$$

In the indentation experiment, n and σ_y_ are unknown variables to be measured, serving as variables to be solved in the system of equations. According to the relationship between the number of equations and the number of unknown variables, when the number of equations is greater than the number of unknowns, the system of equations is overdetermined and has no exact solution. Therefore, when other input and output variables are known, the problem is transformed into finding the least squares solution of the above neural network equation system.

This study proposes an improved iterative optimization algorithm designed based on the LM algorithm for a reverse solution of the nonlinear equation system composed of the aforementioned neural network. The optimization objective function is as follows:

$$[eqn]$$

where y is the predicted output of the model; t is the output of sample data; K is the number of output dimensions, i.e., the number of equations; and N is the number of sample data groups. It is feasible to obtain the solution of the equation system by combining N groups of data. However, in normal experiments, there is no guarantee that the material parameters corresponding to each group of test data are the same. If the data groups are combined, there will be a problem of imbalanced solution conditions and failure of convergence. Therefore, N is selected as 1, and N greater than 1 is only used as a reference candidate. Combining the data groups is feasible when the material parameters of each set of test data are the same.

Therefore, the n and σ_y_ values (obtained after reverse iteration) that yield the smallest value from the optimization objective function are the material parameters to be obtained through inversion based on the nanoindentation P–h curve.

(2)Improvement of the inversion objective function

To solve the system of equations constructed based on a neural network, we designed an improved MSE objective function, shown in Equation (14), in which the weight factors can be adjusted in multi-objective optimization:

$$[eqn]$$

where λ is the error weight coefficient in the range from 0 to 1 (but cannot be 0 at the same time), used to control the contribution of each output error to the iterative solution.

(3)Improvement of the inversion iteration algorithm

For the neural network with an input of X, the iterative solution is based on the following formula:

$$[eqn]$$

where h**lm is the iteration step size of each input variable. According to the LM algorithm, the calculation formula for the iteration step size is as follows [31,32]:

$$[eqn]$$

where μ is the damping parameter in the LM algorithm; H represents the Hessian matrix $[eqn]$ ; and g represents the gradient vector of the loss function $[eqn]$ , where f is the error between the model output and the sample output defined in the loss function.

To achieve a flexible definition of the variables to be solved through inversion in input X, we improved the LM algorithm in the following way: We dot multiplied a logical mask row vector M on the basis of the original Jacobian matrix, where each element of M has a value of 0 or 1, and the length is equal to the dimension number of input parameter X. By assigning a value of 1 to an element M, iterative updating of the input parameter at the corresponding element position is controlled; if the value of 0 is assigned, the iteration step size is always equal to 0, which is equivalent to keeping the initial value of the corresponding element position constant. Thus, the reverse iterative solution of the remaining unknown variables can be achieved when some input quantities and all output quantities are known.

$$[eqn]$$

The convergence criterion for the training process was defined such that the iteration terminates if any of the following conditions is satisfied: reaching the maximum number of training epochs, the loss function falling to the preset target, the gradient norm dropping below the minimum threshold, or the damping coefficient exceeding the upper limit. This setup is designed to balance training efficiency and the stability of model convergence.

3. Application Verification and Discussion

3.1. Model Evaluation

The coefficient of determination R^2^ was used to evaluate the training results of the neural network model used for a reverse solution, which is defined by Equation (18). According to common practice in neural network training, the input and output variables were normalized to the [−1, 1] interval before starting model training. The training ended when the convergence accuracy MSE reached around 1 × 10^−9^. Figure 6 shows that the R^2^ of every output variable of the model relative to the sample data reached 1, indicating that the model can characterize the correlation between the selected inputs and outputs very well.

$$[eqn]$$

When the proposed model is used for reverse solution testing on the to-be-solved input variables n and σ_y_, problems of unstable convergence at some points and multiple solutions can be found. Due to the phenomenon of non-convex optimization in multi-objective and multivariate nonlinear iteration, complex problems such as difficult convergence and multiple extrema are common, which can affect the stability and uniqueness of the reverse solution of model.

The problem of solution stability is manifested in the following phenomenon: significant difference in convergence accuracy exists between the results of different iterations of the to-be-solved variables if different initialization values are used. As a result, failure to converge or premature termination of the iteration can occur, making it impossible to stably obtain the inversion convergence solution every time.

The problem of solution uniqueness is manifested in the following phenomenon: the solution always converges to relatively high precision for the models obtained from different training batches or when the same model is assigned different initial conditions, but each time, the solution is significantly different from the previous one.

Taking the iterative optimization process of the binary variable as an example, the variation in unstable iterative convergence gradients is shown in Figure 7a, the iterative convergence process of non-unique solutions is shown in Figure 7b, and the ideal case is shown in Figure 7c. The solution can always converge to the same extremum point regardless of the initialization values, and the solution is stable and unique.

Based on the above analysis, we devised measures to address the problems of unstable and non-unique solutions in the reverse solution process: evaluate and analyze the sensitivity of input and output parameters using global sensitivity [33], and on this basis, adjust the model training process.

Global sensitivity, represented by Sen_i_, is defined as the degree to which an output parameter is affected by a certain one-dimensional input X**i under the variation in disturbances to other inputs xj≠i in the entire sample analysis space. For each pair of input–output samples ( $[eqn]$ , $[eqn]$ ), the derivative at that sample point is calculated and then summed up. The direct summation of derivatives at each point reflects the directional polarity of the influence, which is defined as vector sensitivity; the accumulation of the absolute value of the derivative at each point reflects the total degree of influence, which is defined as scalar sensitivity. The vector sensitivity and scalar sensitivity can be calculation using the following equation:

$$[eqn]$$

When the input parameters are normalized to the same range, sensitivity analysis reveals that different input parameters have significantly different influences on the outputs regarding sensitivity. The parameters with lower sensitivity (mainly E_r_, n, σ_y_, and S) have a weak relationship with the model output. This indicates that when the normalization range is uniform, there is still a significant scale difference in the influence of the parameters on the outputs. During the model reverse solution process, the influence of initial values can easily cause the model to become trapped in local optima, resulting in convergence failure and unstable solution. Based on the relative relationship of sensitivity under the uniform normalization range for all parameters, we can adjust the hyperparameters by setting narrow normalization range for parameters with lower sensitivity, as shown in the table below.

By comparing Figure 8 and Figure 9, we can see that parameter sensitivity can be enhanced through adjustment and compensation, achieving better overall balance.

A comparison of the sensitivity values before and after adjusting the normalization ranges is presented in Table 3. The effectiveness of the compensation technique proposed in this paper can be verified from the table as follows:

The sensitivity of key constitutive parameters is significantly enhanced: In the balanced sensitivity mode, the influence of material constitutive parameters (yield strength σ_y_ and hardening exponent n) on most output variables is substantially strengthened. For instance:

(a) The sensitivity of σy to the maximum load Pm rose from 0.0232 to 0.1475 (an increase of approximately 535%).

(b) The sensitivity of the hardening exponent n to P_m_ also increased from 0.0495 to 0.1105. This demonstrates that the proposed technique effectively improves the inversion algorithm’s capability to identify the core parameters to be solved.

The sensitivity distribution is more balanced, alleviating the parameter coupling issue: In the unbalanced mode, some parameters exhibit extremely low sensitivity (e.g., the sensitivity of stiffness S to residual depth h_f_ is merely 0.0035), which tends to make the inversion process susceptible to measurement noise or prone to converging to local solutions. After compensation, the sensitivity distribution of all parameters becomes more uniform (e.g., the sensitivity of S to h_f_ is increased to 0.0149), which helps enhance the numerical stability of the inversion system and the reliability of its solutions.

The compensation technique optimizes the exploration efficiency of the parameter space: Balanced sensitivity ensures that the weight of all input parameters’ influence on the outputs is of the same order of magnitude, preventing certain parameters from being “ignored” by the algorithm due to dimensional or numerical range discrepancies. Thus, the inversion process is guaranteed to conduct more efficient and comprehensive searches across the entire parameter space.

In addition, a sensitivity analysis shows that all four selected output variables are correlated with the input variables, and the higher the sensitivity, the stronger the correlation. However, the variables n and σ_y_ have relatively low sensitivities compared with other input parameters. When the outputs are partially used (the number of equations in equation system (12) >2 but <4), it is impossible to iteratively converge to an exact solution, and there is a problem of multiple solutions. Therefore, all four output variables are used for solving unknown variables.

Through the above two adjustment steps, stability and uniqueness of the solution are achieved. As can be seen from the comparison of the relative error distribution for the reverse solution shown in Figure 10, the 95% confidence percentage error is improved from 5% to 3%, and the maximum error is controlled within 8%.

3.2. Physical Rationality and Accuracy Verification of Solution

To verify that the numerical iterative solution conforms to physical constraints, a system of equations is constructed for the elastoplastic constitutive curve’s intersection point under the conditions of n = 0 and n ≠ 0 (Figure 1), and the characteristic strain calculation formula at the intersection point of the curves is obtained (20). The solutions that conform to physical constraints should satisfy the intersection point formula:

$$[eqn]$$

where the characteristic stress σ_r_ at the intersection point is σ_y_2 when n2 = 0.

Based on the prior physical constraints of the intersection point of the stress–strain curve, the characteristic parameters E_r_, P_m_, h_m_, W_l_, and C can be extracted during the loading stage of the material P–h test curve under the condition of n ≠ 0. Then, the characteristic stress σ^R^r at n = 0 (i.e., ideal plastic material) can be solved through model inversion. Furthermore, based on the engineering experience value of ε_r_ = 0.033 (0.025–0.033) and the material parameters n and σ_y_ obtained through inversion, the characteristic stress σ^R^r of the intersection point can be obtained using the intersection point formula. When the characteristic stress values calculated using the two methods are close (difference < 10%), the solution conforms to the physical constraints.

In this study, we measured five sets of fully annealed copper materials at different positions and then conducted inversion to solve material parameters, with the aim of checking to what degree the solution conformed to the physical constraints. The results in Table 4 show that the numerical solution is reasonable. (The five groups of nanoindentation tests in this study were conducted using a Hysitron TI 750 nanoindenter, with a constant maximum load mode (5 mN). The load was held for 10 s once the maximum load was attained, followed by the unloading phase) (see Figure 11).

The material parameters derived from inversion (hardening exponent n and yield strength σ_y_) were systematically compared with the data reported in the published literature for the same material (TSV-Cu) [19,34,35], and the results are presented in Table 5.

For most samples in this study (Samples 1, 2, 4, and 5), the inverted yield strength σ_y_ spans from 35.70 to 46.14 MPa. This range is in excellent agreement with the value of 42.83 MPa reported by Li et al. [19] at the same annealing temperature (300 °C). It is also reasonably consistent with the strength range (47.91–51.83 MPa) of as-deposited electroplated copper cited in other studies [34,35]. This strong consistency provides robust external validation of the high accuracy and reliability of the machine learning inversion algorithm proposed in this study. Our results successfully validate the well-established principle that electroplated copper possesses a higher yield strength than fully annealed bulk copper, stemming from its fine-grained microstructure.

The hardening exponent n of the samples in this study (ranging from 0.5467 to 0.5665) is generally slightly higher than the range reported in the literature (0.49–0.52). This discrepancy may be attributed to the unique microstructures formed by the specific electroplating process employed (e.g., finer grain size), which results in enhanced work hardening capacity. It is worth emphasizing that the inversion method adopted in this study has undergone a sensitivity analysis and normalization range optimization (Section 3.1), leading to a significantly improved capability in identifying parameters such as the hardening exponent n and yield strength σ_y_. Thus, it can capture the mechanical responses that may arise from subtle differences in microstructure.

The yield strength σ_y_ of Sample 3 (86.45 MPa) is notably higher. This is not an experimental or inversion error: the inverted results have passed physical rationality verification (Table 4), and the P–h curve reconstructed using these parameters shows a high degree of agreement with the experimental curve (Figure 12, with a convergence error of <0.1%). This finding precisely demonstrates that the proposed method possesses exceptional detection sensitivity and robustness in resolving the microstructural inhomogeneities of materials (e.g., regions with abnormally fine grains or higher impurity concentrations).

Through a systematic comparison with the literature data, we not only confirmed the reliability of the core results of this study but also highlighted the superiority of the proposed inversion method in resolving subtle differences in material mechanical properties. Our work supplements the existing database and provides a powerful tool for accurately characterizing the process–property relationship.

The trained model was used to perform inversion on the elastoplastic parameters of the selected fully annealed pure copper material. The measured values (Exp), reverse solution (RS) values, and simulation values obtained based on inversion parameters (virtual experimental value FEA) of the nanoindentation P–h curve were compared. The results are shown in Figure 12 and Table 6. The indentation P–h curve obtained in the test almost coincides with the curve plotted based on the parameter values obtained through inversion (or FFA), and the characteristic parameter error of the indentation curve is less than 0.1%.

3.3. Computational Efficiency

A comparison of the NN-LM method and traditional finite element iterative optimization is presented in Table 7.

The traditional FEM relies on repeatedly invoking finite element software for forward calculations, with a single forward simulation taking approximately 10 min. For parameter inversion, dozens to hundreds of iterations are typically required to achieve convergence, resulting in a total inversion time of hours or even longer.

The NN-LM method proposed in this study replaces finite element forward simulation with a well-trained neural network, reducing the forward calculation time from the minute to millisecond scale. Although the Levenberg–Marquardt (LM) algorithm requires an average of 200 iterations to converge to the target error (relative error < 1 × 10^−5^), the high efficiency of neural network forward propagation enables the entire inversion process for a single sample to be completed within 10 s, representing a speedup of more than two orders of magnitude compared with the traditional method.

3.4. Comparison of 3D Packaging Stress

We used a finite element model of local weak positions in a high-density 3D SiP package formed using TCV to analyze and compare the differences between the method proposed in this paper and the simplified constitutive model based on engineering experience parameters regarding the influence on the stress of 3D packaging structures. The results are shown in Figure 13. The thermal expansion coefficient of the ceramic substrate was 7.1 ppm, the elastic modulus was 310 GPa, and Poisson’s ratio was 0.25. The grid was automatically partitioned into 28,570 elements and 53,843 nodes. WorkBench2023 was used to compare the stress conditions of copper-plated materials under the conditions of no structural constraints and constant temperature load (300 °C) using multilinear, bilinear, and power law constitutive models.

The stress distributions of the materials obtained through model simulation are shown in Figure 14. The values of maximum principal stress at the same copper-plated hole obtained using different finite element models are extracted for comparison (Table 8). For the same grid node position, the stress obtained using the bilinear model is the highest, followed by that of the multilinear model, and the stress obtained using the power law model is the lowest (about 88%).

4. Conclusions

A neural network model for high-precision characterization of material power law constitutive model parameters was constructed based on characterization processing of the P–h curves plotted based on the data obtained in in situ near non-destructive testing of the parameters of microelectronic packaging materials using the nanoindentation testing technique. Then, a novel inversion algorithm was designed based on the neural network model for fast, accurate, and stable extraction of material mechanics parameters. The transformation process between the neural network model and overdetermined equation system inversion solution was worked out, and an inversion objective function with weighted factors and an improved LM algorithm with mask vectors were constructed to solve unknown inputs when some inputs and outputs are partially known. In addition, the least squares iteration of complex nonlinear equations between multiple variables was analyzed, along with the solution stability and uniqueness of inversion. By analyzing the differences between the influence of inputs on outputs regarding sensitivity, the normalization range of the modeling data was optimized to improve the stability of the reverse solution. The physical rationality of the reverse numerical solution was verified based on the proposed prior physical constraints of the elastoplastic constitutive curve’s intersection point. A comparative discussion of the material parameters derived from inversion against the data reported in the published literature not only verifies the reliability of the proposed method and the rationality of its results but also situates this study within the broader context of material characterization research, highlighting its significant value in the accurate elucidation of process–property relationships. A comparison with the results obtained through finite element simulation indicated that the proposed method can achieve an accuracy of <3% (95% confidence interval) in measuring material parameters, with a maximum error of less than 8%. A simulation of the stress of an actual packaging structure verified that the measured material parameters and the differences between material models have a significant impact on the stress, which provides useful guidance for optimization of the packaging process.