A RDL Modeling and Thermo-Mechanical Simulation Method of 2.5D/3D Advanced Package Considering the Layout Impact Based on Machine Learning

1. Introduction

Redistribution layer (RDL) is one of the most important components of the advanced packages like Chip on Wafer on Substrate (CoWoS) [1], Foveros [2], X-Cube [3], System on Integrated Chips (SoIC) [4,5], especially the multilayer RDL interposer based CoWoS-R [6], FO-PLP [7] and chiplet-based architectures [8,9]. To enable high-bandwidth signal processing between chiplets [10], several techniques such as dual damascene processing, semi-additive processing (SAP) and polymer damascene processing have been introduced. These techniques narrow the RDL trace width/pitch to as small as 2 microns [11,12] and even sub-micron [13]. To enhance the electricity transmission capacity, traces become thicker, with an aspect ratio of up to 4.2 [14], which is very different from traditional shell-like thin traces. Furthermore, wafer level integrations can stack up to 5 metal layers [1], with the area up to 1200mm² [15], or even as extensive as 2500mm² [1]. As a result, RDLs in advanced packages become extremely complex in shape and behavior, posing great challenges in accurately capturing its thermo-mechanical characteristics. Additionally, in advanced packages as substrate is extremely thinned to several microns or even sub-micron [16] which is similar to or even thinner than RDL, making the RDLs more prominent comparing to other package components [17,18]. For this reason, RDLs are playing crucial roles in package thermo-mechanical characteristics.

Integrated packages face severe thermo-mechanical risk due to the coefficient of thermal expansion (CTE) mismatch between materials, the accelerated thermal cycling test (TCT) [19] is an important approach for assessing their reliability. As a low-cost alternative of testing, finite element analysis (FEA) simulation under TCT conditions becomes a significant research area on design for reliability (DFR). Balancing between efficiency and accuracy of simulation is a critical issue in this field.

For traditional integrations, due to the thicker substrate and larger trace width/pitch, axisymmetric FEA model can meet the efficiency and accuracy requirements for DFR. By 2D axisymmetric model, Lee et al. [20], Che et al. [21], Machani et al. [22] simplified the RDLs into homogeneous rectangles and refined the outermost bump, then evaluated the fatigue life of the bump through transient simulation. By 3D axisymmetric model, Lee et al. [23], Che et al. [24] simplified the RDL into homogeneous films and refined the corner solder, thus assessed the reliability of stacked-chip packages. However, advanced packaging has different structural characteristics.

For advanced packages, RDL plays a more significant role in thermo-mechanical reliability. Therefore, more precise RDL modeling approaches are required for DFR of advanced packages. Dividing the RDL then employing material equivalization is a widely used strategy to improve the accuracy. In this strategy, division size and equivalization method are crucial for efficiency and accuracy. When efficiency is concerned, RDL will be divided into large size regions, and there are three typical material equivalization methods. The first method, employed by Wang et al. [25], involves calculating the isotropic materials based on the copper ratio within each region. The second method, used by McCaslin et al. [26], is determining the anisotropic equivalent material properties of regions consisting of simple shape traces through composite material mechanics. The third one, utilized by Valdevit et al. [27], is the morphology-based approach, where regions are categorized into 'lines', 'vias', or 'web' types, then specific mathematical models are developed to compute the material properties for each type of region. When accuracy is concerned, RDL will be divided into small size blocks, and there are also three typical material equivalization methods. The first method, employed by Lien et al. [28], involves calculating the isotropic equivalent material properties based on the metal percentage within each block. The second method, used by Lee et al. [29], focuses on calculating the anisotropic material properties and the equivalent reference temperature of blocks using composite material mechanics. The third one, introduced by Gibson [30], and Lee et al. [29], is a simulation-based approach that calculates the anisotropic equivalent material properties of blocks containing complex traces. To achieve higher accuracy, Yaddanapudi [31] divided the RDL into extremely small size pixels and directly defined each pixel as either metal or dielectric, a high-fidelity structural response of RDL was achieved by consuming large amount of solving resources. In general, existing work indicates that there is a trade-off between division sizes and the material equivalization method. The key to balancing efficiency and accurately lies in obtaining the anisotropic equivalent material properties quickly and accurately under larger division size.

Recently, machine learning (ML) has shown the capability of extracting interpretable models from scientific data automatically [32], and it has been increasingly employed in the material equivalization of composite materials in 3D and 2D structures. Regarding 3D, Dai et al. [33] used feature matrixes to digitally represent the lattices orientations of 3D structures, and then created a graph-Artificial Neural Network (ANN) to predict the structural performance. Regarding 2D, Liu et al. [34], Ye et al. [35] and Gong et al. [36] used high-resolution grayscale matrix to digitalize the planar strictures, then constructed convolutional neural networks (CNNs) to quickly predict the equivalent mechanical properties. Due to the accurate and fast predictive capabilities of models created based on ML, it has been introduced into the modeling and simulation of integrated package. For example, Selvanayagam et al. [37,38] partitioned an 8-RDLs interposer into 3×3 regions, then established surrogate models that link the copper ratio in each layer within each region to the warpage in TCT condition, then ultimately optimized the global package warpage.

In this work, we applied ML for material equivalization and developed a RDL modeling and simulation method for DFR of advanced packages. In this method, a RDL is divided into equal-size blocks. To describe the layout within each block, we further subdivided each block into pixels and represent each pixel using the metal percentage. Consequently, each block can be represented as a tensor composed of pixel metal percentages, and the entire RDL is digitalized as the collection of such block tensor. An ANN is constructed using either fully connected neural networks (FCNN) or convolutional neural networks (CNN), to model the relationship between the block tensor and the corresponding equivalent material properties. Subsequently, the ANN is trained by a small subset of tensors with the equivalent material properties obtained through finite element analysis (FEA) as labels, then the material surrogate model is constructed. Lastly, all blocks are transformed into FEA elements and assigned the material properties predicted by the model, then completing the modeling and simulations of the RDL.

To validate the efficiency and accuracy of the proposed method, line bending simulations were conducted for a 21.6×21.6mm² RDL. The detailed fine mesh model served as the benchmark, with the reaction force being regarded as the accuracy indicator. The results show that as the reduces of substrate thickness, the reaction force error of the traditional method that neglects layout impact increases sharply, indicating that critical layout impact must be considered in the simulation of advanced package. When the proposed method is used with a substrate thickness of 50μm, the reaction force error is 2.81% the layout impact can accurately be estimated utilizing only 200×200 elements, the number of elements is only 1/15 of the traditional approach.

As an application case, a simulation was introduced to investigate the thermo-mechanical response of a 2.5D integrated CPU chip at the maximum temperature state in TCT. The results indicate that for advanced packages the maximum stress more likely to occur in the RDL, which is different from the traditional integrations where the maximum stress occurs in the bumps. It also reveals that the stress of both the RDL and bumps are significantly impacted by the layout. Moreover, it is observed that the stress in RDL is particularly impacted by vias and bumps in adjacent layers.

The method precisely concerns these layout-related impacts with minimal resources and time, presents an opportunity to improve the efficiency of advanced package DFR.

2. Methods

The general workflow of the proposed modeling method is shown in Figure 1 and consists of the following steps:

Step 1: Two-level RDL digitalization, this step is consisted by global and local discretization. In the global level, the package level RDL pattern is divided into equal-size RDL blocks containing a part of traces, and there are Q_x blocks along the X direction and Q_y along Y. To capture the layout feature within block, in the local level, each RDL block is divided into pixels with uniform size. The subdivision results in PQ_x pixels along the X direction and PQ_y pixels along Y, the pixel value is defined as the metal percentage within pixel region. Consequently, the RDL block can be represented by a PQ_x×PQ_y tensor, and the package-level RDL can be digital represented by Q_x×Q_y tensors.

Step 2: ANN Training dataset preparing, this step prepares the ANN training dataset by randomly selected a small subset of RDL blocks. The dataset consists of input tensors and labels representing the anisotropic equivalent material properties of each block obtained through 3D FEA simulations [16], including Young's modules, Poisson's ratios, shear modules, thermal conductivities and CTEs in all directions.

Step 3: ANN-based surrogate model establishing, this step establishes an ANN-based material surrogate model. The input data of the surrogate model is the RDL block tensor and the output data include all equivalent properties. The surrogate model employs ANN and is trained by the dataset prepared in Step 2. Once established, the ANN eliminates the need for repeating Steps 2 and 3, it would be highly advantageous for the repeated design iteration process.

Step 4: Equivalent properties predicting, all equivalent properties for all the RDL blocks generated in Step 1 are predicted by the surrogate model established in Steps 2, 3.

Step 5: Global FEA model building, all RDL blocks are transformed into 3D hexahedral solid elements to construct the global RDL FEA model and assigned the material properties predicted by Step 4. The model consists of Q_x×Q_y elements encompassing the entire RDL, and can also serve as a part of the package-level FEA model.

As the method implementation, a series of programs are established, including: the ANN implemented by Pytorch, the main framework developed in Python and the FEA simulator powered by ANSYS 18.2. The RDL blocks creating processes are implemented by C++. In addition, the C# developed program is used as the graphical post-processor.

2.1. ANN Architectures

Activation function, rectified linear unit (ReLU) is selected as the activation function, shown by (1), (2):

$$\mathrm{f}\left(\mathrm{x}) =\max (0, \mathrm{x}\right)$$

(1)

$$\mathrm{f}\text{'}\left(\mathrm{x}\right)=\left\{\begin{array}{c}1,\mathrm{x}\ge 0\\ 0,\mathrm{x}<0\end{array}\right.$$

(2)

Loss Function, cross-entropy loss is selected as the loss function, shown by (3):

$$\mathrm{Loss}=-\sum _{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{w}}_{\mathrm{i}}{\mathrm{y}}_{\mathrm{i}}\log \widehat{{\mathrm{y}}_{\mathrm{i}}}$$

(3)

Where, n is the number of samples, y_i represents the true value of the i-th sample, $\widehat{{\mathrm{y}}_{\mathrm{i}}}$represents the predicted value of the i-th sample, w_i represents the weight of the i-th sample. The loss for a batch prediction, shown by (4):

$$\mathrm{Loss}=-\frac{1}{\mathrm{batch}}{\sum }_{\mathrm{j}=1}^{\mathrm{batch}}{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{w}}_{\mathrm{ij}}{\mathrm{y}}_{\mathrm{ij}}\log \widehat{{\mathrm{y}}_{\mathrm{ij}}}$$

(4)

Where, batch is the number of predicting procedures for one predicting batch. For the j-th predicting procedure of the batch, y_ij represents the true value of the i-th sample, $\widehat{{\mathrm{y}}_{\mathrm{ij}}}$represents the predicted value of the i-th sample and wi represents the weight of the i-th sample.

Optimizer, the stochastic gradient descent (SGD) is adopted to optimize neural network, and momentum is considered during the optimization process. Learning rate and momentum parameters are adjusted to fine-tune the optimization process, as indicated in equations (5) and (6)

$${\mathrm{v}}_{\mathrm{t}}=\mathrm{m} \times {\mathrm{v}}_{\mathrm{t} - 1} - \mathrm{lr} \times \mathrm{g}$$

(5)

$${\mathrm{w}}_{\mathrm{t}}={\mathrm{w}}_{\mathrm{t} - 1}+{\mathrm{v}}_{\mathrm{t}}$$

(6)

Where, t is the current iteration step number, relatively t–1 is the number of the previous step, lr represents the learning rate, g is the gradient, m means the momentum. v_t, w_t represents the velocity and weight in current iteration step respectively. v_t-1, w_t-1 represents the velocity and weight in the previous step respectively.

ANN networks, full-connected neural network (FCNN) and convolution neural network (CNN) are created. As shown in Figure 2, the FCNN includes an input layer, output layer and several hidden layers. The parameters of the input layer are represented by the vector shown by (7) and the size of the input vector shown by (8):

$$\mathrm{I}=[{\mathrm{I}}_1^{\left[0\right]}, {\mathrm{I}}_2^{\left[0\right]},\dots ,{\mathrm{I}}_{\mathrm{p}}^{\left[0\right]},\dots ,{\mathrm{I}}_{{\mathrm{Q}}_{\mathrm{in}}}^{\left[0\right]}]$$

(7)

$${\mathrm{Q}}_{\mathrm{in}}={\mathrm{PQ}}_{\mathrm{x}} \times {\mathrm{PQ}}_{\mathrm{y}}$$

(8)

where, ${\mathrm{I}}_{\mathrm{p}}^{\left[0\right]}$ is the metal percentage of the p-th pixel, Q_in is the vector size which is equal to the number of pixels of a single RDL block, PQ_x and PQ_y represent the number of RDL block’s pixels along the X and Y direction respectively. The output vector of the FCNN includes all 15 anisotropic equivalent material properties of the RDL block, shown by (9):

$$\mathrm{O} =\left[{\mathrm{E}}_{\mathrm{x}}, {\mathrm{E}}_{\mathrm{y}}, {\mathrm{E}}_{\mathrm{z}}, {\mathsf{\mu }}_{\mathrm{xy}}, {\mathsf{\mu }}_{\mathrm{xz}}, {\mathsf{\mu }}_{\mathrm{yz}}, {\mathrm{G}}_{\mathrm{xy}}, {\mathrm{G}}_{\mathrm{xz}}, {\mathrm{G}}_{\mathrm{yz}}, {\mathrm{k}}_{\mathrm{x}}, {\mathrm{k}}_{\mathrm{y}}, {\mathrm{k}}_{\mathrm{z}}, {\mathrm{A}}_{\mathrm{x}}, {\mathrm{A}}_{\mathrm{y}}, {\mathrm{A}}_{\mathrm{z}}\right]$$

(9)

The hidden layer can be represented by vector as (10):

$${\mathrm{A}}^{\left[\mathrm{n}\right]}=[{\mathrm{A}}_1^{\left[\mathrm{n}\right]}, {\mathrm{A}}_2^{\left[\mathrm{n}\right]},\dots ,{\mathrm{A}}_{\mathrm{i}}^{\left[\mathrm{n}\right]},\dots ,{\mathrm{A}}_{{\mathrm{H}}_{\mathrm{n}}}^{\left[\mathrm{n}\right]}]$$

(10)

where, ${\mathrm{A}}_{\mathrm{i}}^{\left[\mathrm{n}\right]}$is the value of the i-th node, Q_n is the number of nodes of the n-th hidden layer.

Shown as Figure 3, the CNN is consisted by convolution network and the full-connected network, the input parameters are the tensor with size of [PQ_x, PQ_y, 1], and the output vector is the same as (7). The convolution network contains convo lution layers and max pooling layers and the tensor is transf ormed through layers. The tensor size is [C_nx, C_ny, C_nz] after the n-th convo lution layer, and [P_nx, P_ny, P_nz] after the n-th max pooling layer. The output tensor of the convolution network can be compressed into the vector of size Q_Lin.

$${\mathrm{Q}}_{\mathrm{Lin}}={\mathrm{P}}_{\mathrm{nx}} \times {\mathrm{P}}_{\mathrm{ny}} \times {\mathrm{P}}_{\mathrm{nz}}$$

(11)

where, P_nx, P_ny, P_nz are the output tensor size in X, Y and Z directions respectively.

2.2. Training dataset augmentation

To augment the training dataset of blocks, the geometric symmetry and data transformations are employed. The initial block with anisotropic properties (E_x = E_x0, E_y = E_y0) is shown in Figure 4a, and it can be rotated by 90°, 180°, and 270°, as indicated in Figure 4b,d, respectively. When the angle of rotation is 90° or 270°, the anisotropic properties are transformed to E_x = E_y0, E_y = E_x0. The initial block and its rotated versions can also be flipped in the Y direction, as shown in Figure 4e–h, while keeping the same anisotropic properties. Eight unique blocks can be generated from the initial one, without requiring additional FEA solutions.

3. Result and discussions

3.1. Method validation

The two-layer model containing RDL and a silicon substrate is shown in Figure 5 is established for validation. The model consists of the substrate layer and the RDL (the M1 layer), the materials involved in the simulations are listed in Table 1, and all materials are assumed to be linear elastic. The substrate is made up of Si, its thickness (T_B) is set as a variable parameter, with a range of values from 20μm to 500μm. The RDL is consisted by Cu and PI, the thickness (T_R) is 10μm. The layout pattern of M1 layer is shown in Figure 6 and the domain direction of traces is the X direction.

The validation of the adaptability to the mechanical loads is based on the parallel trace bending load case (CaseP) and vertical trace bending load case (CaseV) shown in Figure 5. In the CaseP, nodes on X = 0 face are constrained in the X, Y and Z directions, nodes on the top line of X = Length face are set Z = -0.5 mm. In the CaseV, the constrain are set on Y = 0 face, and the displacement load are set on the top line of Y = Width face.

As the benchmark, a detailed RDL fine mesh model consisting of 640,365 nodes and 606,482 elements was established with a substrate thickness of 50μm. The solving process cost 176s. The stress distributions are depicted in Figure 7, stress concentrations can be observed at discontinue positions, indicating the critical impact from layout.

The sums of the Z-direction reaction forces of all nodes applied displacement loads in CaseP and CaseV are FR_DP: -7.174E-03N and FR_DV: -6.464E-03N, respectively. The reaction force errors serve as accuracy indicators, as shown in (12).

$$\left\{\begin{array}{c}{\mathrm{Error}}_{\mathrm{P}}=\left|({\mathrm{FR}}_{\mathrm{P}} - {\mathrm{FR}}_{\mathrm{DP}}) /{\mathrm{FR}}_{\mathrm{DP}}\right|\times 100\%\\ {\mathrm{Error}}_{\mathrm{V}}=\left|({\mathrm{FR}}_{\mathrm{V}} - {\mathrm{FR}}_{\mathrm{DV}}) /{\mathrm{FR}}_{\mathrm{DV}}\right|\times 100\%\end{array}\right.$$

(12)

Where, Error_P and Error_V are reaction force errors for CaseP and CaseV respectively. FR_P and FR_V are reaction force to be evaluated for CaseP and CaseV, respectively.

For the proposed method, considering the efficiency and accuracy set the global mesh divisions: Q_x = Q_y = Q = 200, the local pixel division: PQ_x = PQ_y = PQ = 15 and built 40,000 RDL blocks. Q_t blocks are randomly selected as the ANN training dataset, the default value of Q_t is 1000, and the testing dataset contains all the 40,000 blocks. A FCNN with two hidden layers is employed for the ML based method. As depicted in Figure 2, the first hidden layer contains 400 nodes (H1 = 400), while the second one has 300 nodes (H2 = 300). Set the learning rate: lr = 0.02, the momentum: m = 0.5, the training epoch: q = 250. The training and testing processes are run using a GPU of NVidia RTX 3070Ti, the average predicting time of all blocks are 3.652 s, and the loss variation is shown by Figure 8.

The calculated material properties shown in Figure 9 indicate the ML-based method predicted material properties are very close to those obtained by the simulation-based method, which demonstrates the accuracy of the ANN-based surrogate model. Additionally, notable difference between Young’s modulus along the trace direction (E_x) and perpendicular (E_y). Moreover, reduction in E_y is obvious within the gap area between traces, indicating the critical impact of RDL layout on the mechanical behavior of RDL.

The stress distribution shown in Figure 10 is well consistent with the result of the detailed fine mesh model presented in Figure 7. Additionally, the reaction force errors for CaseP and CaseV are 1.72% and 2.81%, respectively, validating the precision of this method. Furthermore, the proposed approach only requires 200×200 FEA elements, whereas the detailed fine mesh model consists of 15 times more elements. This demonstrates that the approach contributes significant resource savings in simulation.

Furthermore, to verify the accuracy of thermo-mechanical simulations, parallel trace warpage case (CasePT) and vertical trace warpage case (CaseVT) are created as shown in Figure 11. Except to the 150 ℃ temperature, nodes on X = 0 and Y = 0 faces are constrained in the X, Y and Z directions in CasePT and CaseVT cases, respectively. Figure 12 and Figure 13 show the Z direction displacement and von Mises stress distributions of the cases, they present almost identical deformation and stress distributions to the detailed fine mesh model.

3.2. Key factors influence

The influence of the substrate thickness, global mesh division, training dataset size, and ANN architecture to the prediction process is investigated.

The influence of the substrate thickness.

Three simulation methods were performed, including the traditional layout neglecting block-based method, the simulation-based method and the proposed ML-based methods, on seven substrate thickness values: 20, 50, 100, 200, 300, 400, and 500μm. Figure 14 summarizes the reaction force errors of each case, demonstrating the fact that as the substrate thickness decreases, the impact of RDL layout becomes more and more critical, neglecting the layout features will result in substantial errors.

The influence of the global mesh division.

Another two global division values of 100 and 300 are conducted to investigate the influence of the global division. Figure 15 summarizes the reaction force errors of different cases, illustrating that larger global mesh divisions can improve solution accuracy, albeit at the cost of increased solving time and resource consumption.

Figure 16 depicts stress distributions at Q = 100 and Q = 300, the stress distribution shown in Figure 16a is similar to those shown in Figure 7 and Figure 10, indicating that the proposed method can accurately capture the layout impact shown in Figure 6 using only 100×100 elements. Moreover, larger Q values can capture more accurate pattern features but may not justify the increased time and resource consumption. Regarding the block-based method, as shown in Figure 17, obvious stress inhomogeneity only becomes captured until Q ≥ 600. Moreover, even when Q = 800, Figure 18 demonstrates that the accuracy cannot match the proposed method. However, simulation costs become exorbitant.

The influence of the training dataset size.

ANNs was established using ten different dataset sizes ranging from 500 to 5000. Figure 19 summarizes the loss and training time and indicates that increasing dataset size generally results in improved prediction accuracy, although the effect becomes negligible beyond 2000 data points. Furthermore, the proposed dataset augmentation algorithm can significantly enhance the prediction accuracy as shown in Figure 20.

The influence of the FCNN architectures

Three different types of FCNNs with two hidden layers and varying numbers of hidden nodes are established. The training and prediction losses are summarized in Figure 21a and show that increasing the number of hidden layer nodes results in improved prediction accuracy. However, the effect becomes less significant as the number of nodes increases. Additionally, Figure 21b illustrates the accuracy results obtained from the FCNNs with three hidden layers and varying numbers of nodes, showing that while additional layers can improve accuracy, their benefit is limited and requires more training time.

Application of the CNN architectures.

Figure 22 summarizes the losses and training time for “1 CNN + 1 FCNN”, “1 CNN + 2 FCNNs” and “2 CNNs + 2 FCNNs”. The results demonstrate that while CNNs have shorter training time, they exhibit significantly larger losses across all three architectures, indicating undesired defeaturing during convolution and pooling operations negatively impacts prediction accuracy.

3.3. Large area 2.5D integrated CPU chip thermo-mechanical simulation

As an application, a 2.5D integrated CPU chip is constructed, as shown in Figure 23. All 3 RDLs are modeled by the proposed method, the division parameters are: Q = 200, PQ = 15. The training dataset size of Q_t = 1000 is utilized without applying data augmentation. Predictions are made using FCNNs with two hidden layers consisting of 400 and 300 nodes, respectively. In the package-level FEA model, the maximum element size is 0.12mm, it contains 1,072,658 nodes and 625,648 elements.

To simulate the maximum temperature state of the TCT, a temperature load of 150℃ was applied and the bottom surface of the BGA was constrained in X, Y and Z directions, as shown in Figure 24. The results obtained under this condition is shown in Figure 25, indicating following information:

Firstly, it is found from Figure 25a that the maximum stress occurs in the M1 layer rather than bumps, this indicates that for 2.5D advanced packages, as the substrate thickness decreases and the trace width/pitch reduces, the thermo-mechanical risk of RDL becomes more significant. The conclusion for traditional integrations that the maximum risk always occurs at bumps is no longer accurate.

Secondly, Figure 25b,c demonstrate the fact that the stress distribution of RDL is critical impacted by layout, neglecting the layout impact would bring in critical errors. And the proposed method can accurately capture the impact.

Furthermore, it can also be found from Figure 25(a) that the location and magnitude of the maximum stress of bumps and TSVs are influenced by the layouts and surrounding structures, and it is different from the conventional understanding that the outermost bumps are the most vulnerable. Therefore, the impact of the vias and bumps layouts should also be considered in the simulation.

Finally, the detailed view of the maximum stress region of M1 layer in Figure 25d illustrates that the stress on RDL is not only influenced by the RDL layout but also critical impacted by the vias and bumps within adjacent layers.

4. Conclusions

This paper proposes a machine learning (ML) based RDL modeling and simulation method for DFR of 2.5D/3D advanced packages, taking the critical layout impact into consideration. In the method, RDL is divided into blocks and subdivided into pixels of metal percentage, such that the RDL can be digitalized as metal percentage tensors. Consequently, an ANN-based surrogate model is built and trained by a small subset of tensors to predict the equivalent material properties of each block. Lastly, all blocks are transformed into elements to build a FEA model for simulations. For the validation of accuracy, line bending simulations are conducted for a 21.6×21.6mm² RDL using the proposed method, where the reaction force is compared with a benchmark result obtained by detailed fine mesh simulation. The results show that through the proposed method, the reaction force error can be as low as 2.81% and the layout impact can be accurately considered with 200×200 elements, just 1/15 of the traditional approach. The results also show that if the layout impact is neglected, the reaction force error will increase sharply as the substrate thins. Therefore, the impact must be considered in advanced package simulation. As an application case, a simulation was introduced to investigate the thermo-mechanical response of a 2.5D integrated CPU chip at the maximum temperature state in TCT. The results indicate that for advanced packages the maximum stress more likely to occur in the RDL, which is different from the traditional integrations where the maximum stress occurs in the bumps. It also reveals that the stress of both the RDL and bumps are significantly impacted by the layout. Moreover, it is observed that the stress in RDL is particularly impacted by vias and bumps in adjacent layers. The proposed method precisely concerns layout impacts with little resource, presents an opportunity to improve the efficiency of advanced package DFR.