Application of Thermo-Viscoelastic Laminated Plate Theory to Predict Warpage of Printed Circuit Boards

1. Introduction

Along with rapid assembly for the purpose of creating thinner printed circuit boards, the side effect of warping during the reflow process is inevitable. As a result, the assembly process encounters serious challenges, such as difficulty in implementing a highly integrated assembly, as well as degradation in the reliability of connectivity. Aside from attempts to simulate these issues employing the FE analysis method, it is also necessary to conduct an estimation of the warpage at the early stage of design, for which development of more simplified estimation tools is strongly desired. From a material behavior point of view, if any glass transition points exist within the temperature range of the reflow process, the relaxation effect (due to viscoelastic characteristics of the material of the boards) appears as a deformation, which acts as a barrier to achieving a distinct estimation of the amount of warpage. Generally speaking, the viscoelastic behavior of resin materials exhibits very sensitive temperature dependency. This makes it difficult to accurately capture the characteristics of resin materials in actual experiments, and accordingly, to build numerical models based on such an experimental measurement. There have been many analysis cases published using sophisticated FE approaches (Shrotriya et al., 2005 and Valdevitet al., 2008), but these approaches have not extended beyond applications as a handy design tool.

To easily estimate the thermal deformation in the laminated structure, some approaches employ the multilayered beam theory (Lim, 2008 and Yuju, 2003). The method proposed in this paper enhances these approaches to incorporate the layered plate theory, and includes the effect from the temperature-dependent viscoelasticity and the temperature-dependent coefficient of thermal expansion of resin materials. The program, which is equipped with the developed method, can give an arbitrary temperature history to a multilayered plate consisting of an arbitrary number of layers. As well, the practical approach for measurement of viscoelasticity characteristics has been specifically developed by the author’s group, with the aim of performing measurements to obtain highly accurate data (Kobayashi, 2008). The method proposed herein was verified using the achievements from this development. Verification results are also shown in this paper. This method instantly calculates the amount of warpage and stress in a multilayered plate by giving the values of the thickness and material constants of the plate. Advantages of the proposed method may cover a wide range of real world applications, such as design optimization problems.

2. Viscoelasticity theory and its incremental form solution

In order to incorporate viscoelasticity into the multilayered plate theory, the one-dimensional linear viscoelasticity theory can be expanded into the plane stress field, and then treated as an incremental form of solution. The generalized Maxwell model is applied to exhibit linear viscoelastic behavior. This model is composed of a parallel series of multiple Maxwell models, each of which is assembled with a serial connection of a dashpot and a spring defined by:

$$E_r\left(t^{\prime }\right)=E_e+{\displaystyle \sum _{n=1}^N{E}_n\exp \left(-t^{\prime }/{\tau }_n\right)}$$E1

(1)

where E_r(t’) is called the relaxation modulus of longitudinal elasticity defining the stress relaxation keeping the strain constant. E_n, τ_n, and N are the material constants denoting the coefficient to the Prony series, the relaxation time, and the number of terms in the Maxwell model, respectively; t' is the reduced time. Where the reduction rule of time-temperature is applicable, the time-temperature conversion factor, a_T(T), is obtainable from the following formula:

$$t^{\prime }={\displaystyle {\int }_{\text{}0}^t\frac{1}{a_T\left(T\right)}dt}$$E2

(2)

where t is the real time and T is the temperature.

Using Eq. (1) and Eq. (2), the one-dimensional behavior of the viscoelastic material with temperature dependency can be represented. Considering the plane stress state, where uniform in-plane deformation is assumed to take place, the stress-strain (σ-ε) relation can be represented with the relaxation form as follows:

$$\sigma \left(t\right)=\frac{1}{1-\nu }{\displaystyle {\int }_0^t{E}_r\left(t^{\prime }-{\tau }^{\prime }\right)}\frac{\partial \epsilon \left(\tau \right)}{\partial \tau }d\tau$$E3

(3)

where ν is the Poisson’s ratio. This study is conducted in accordance with the plate theory of shin shells assuming the thicknesses of a circuit board are not relatively large. Based on this pre-condition, the effect from the out-of-plane shear deformation of the plate is not taken into consideration. Therefore, the bending behavior of the plate is mainly governed by only the Young’s modulus. Note that Poisson’s ratio is incorporated into Eq. (3) for the purpose of expanding the use of the beam theory into the plate theory. Accordingly, the Poisson’s ratio in this study is treated as constant value. Further, the Poisson’s ratio of viscoelastic materials in its strict sense is allowed to have its own time-dependency independent of the Young’s modulus, such that the Poisson’s ratio can be experimentally measured and the measured Poisson’s ratio becomes available for the associated analysis.

Equation (3) should be converted into an incremental form so that it can be treated numerically. Taking t_m as arbitrary time, and assuming that the strain varies with a constant gradient of Δε(t_m)/Δt_m during a time increment of Δt_m= t_m-t_m-1, the relationship between the stress increment, Δσ(t_m), and the strain increment, Δε(t_m), can be represented using Eq. (1) and Eq. (3) by following the expression below (Eq. (4)). Equation (4) takes a form with respect to the strain increment; in the multilayered plate theory described in the section below, necessary equations are derived with respect to the continuity of the strain in each layer, where σ_n is the stress in the n-th term of the Maxwell model.

$$\Delta \epsilon \left(t_m\right)=\frac{1}{A}\left\{\Delta \sigma \left(t_m\right)+{\displaystyle \sum _{n=1}^N{a}_n\left(\Delta {t^{\prime }}_m\right){\sigma }_n\left(t_m-\Delta t_m\right)}\right\}$$E4

(4)

where

$$\{\begin{array}{c}\Delta {t^{\prime }}_m={t^{\prime }}_m-{t^{\prime }}_{m-1}\\ \begin{array}{l}{a}_n\left(\Delta {t^{\prime }}_m\right)=1-\exp \left(-\Delta {t^{\prime }}_m/{\tau }_n\right)\\ b_n\left(\Delta {t^{\prime }}_m\right)=a_n\left(\Delta {t^{\prime }}_m\right)\cdot {\tau }_n/\Delta {t^{\prime }}_m\end{array}\\ A=\frac{1}{1-\nu }\left\{E_e+{\displaystyle \sum _{n=1}^N{b}_n\left(\Delta {t^{\prime }}_m\right)\cdot E_n}\right\}\\ {\sigma }_n\left(t_m-\Delta t_m\right)\\ =\frac{1}{1-\nu }{\displaystyle {\int }_0^{t_m-\Delta t_m}{E}_n\exp \left\{\frac{-{t^{\prime }}_m+\Delta {t^{\prime }}_m+{\tau }^{\prime }}{{\tau }_n}\right\}\frac{\partial \epsilon }{\partial \tau }d\tau }\end{array}$$E5

(5)

3. Multilayered plate theory including viscoelasticity

3.1. Assumptions

For modeling a printed circuit board with a multilayered plate as illustrated in Fig. 1, the following assumptions were made:

• The in-plane property of the plate is homogeneous and isotropic.

• The thickness of the plate is sufficiently thin to generate a no stress component in the direction normal to the surface.

• The plate is not subjected to any constraint.

• The plate will warp under the uniformly distributed temperature.

The bending moment to deform the plate convexly is defined as positive; the deflection along this direction is also defined as positive. H_i and ν_i in the figure denote the plate thickness and the Poisson’s ratio, respectively. The curvature induced under temperature variation becomes identical in all in-plane directions of the xy plane, since isotropicity and non-constraint are assumed. Therefore, the warpage of the multilayered plate can be calculated only on the deformation within the xz plane, as shown in Fig. 2.

Figure 1.

Multilayered plate

Figure 2.

Deformation in xz plane

In multilayered plate theory, global warpage is calculated based on the assumption that the strain in each layer is independent and the respective interface is continuous; the strain generated in each of the layers is composed of three components:

1. Thermal strain

2. Strain due to in-plane force in the plate

3. Strain due to bending moment on the plate

Each strain component is classified as follows:

1. Thermal strain

When the temperature-dependent thermal expansion coefficient of the i-th layer at an arbitrary temperature, T, is expressed as α_i(T), the thermal strain is expressed as the following equation:

$${{\epsilon }^{\prime }}_i={\alpha }_i(T)\cdot (T-T_0)-{\alpha }_i(T_1)\cdot (T_1-T_0)$$E6

(6)

where T₁ represents the initial temperature in the analysis. The temperature-dependent thermal expansion coefficient is the average thermal expansion coefficient based on a reference temperature, T₀. The second term in Eq. (6) is necessary in order to prevent strain at the initial temperature from being generated. In order to express it as an incremental form, Eq. (6) is differentiated by T.

$$\frac{d{{\epsilon }^{\prime }}_i}{dT}=\frac{d{\alpha }_i(T)}{dT}\cdot (T-T_0)+{\alpha }_i(T)$$E7

(7)

Therefore, the thermal-strain increment in a given temperature change, ΔT (t_m), can be obtained from the following equation:

$$\Delta {{\epsilon }^{\prime }}_i(t_m)=\left\{\frac{d{\alpha }_i(T)}{dT}\cdot (T-T_0)+{\alpha }_i(T)\right\}\Delta T(t_m)$$E8

(8)

1. Strain due to in-plane force in the plate

On each layer of the multilayered plate illustrated in Fig. 2, generation of in-plane force and bending moment are considered. Taking the in-plane force in the i-th layer as P_i(t_m), the incremental force can be expressed with Eq. (9). H_i and B are the thickness and width of each layer, respectively; σ’’_i(t_m) indicates the stress due to the in-plane force.

$$\Delta P_i\left(t_m\right)=H_i\cdot B\cdot \Delta {{\sigma }^{″}}_i\left(t_m\right)$$E9

(9)

Introducing the effect from viscoelasticity using Eq. (4), as previously derived, Eq. (8) can be written in the incremental form as follows;

$$\Delta {{\epsilon }^{″}}_i\left(t_m\right)=\frac{1}{H_i\cdot B\cdot A_i}\left\{\Delta P_i\left(t_m\right)+{\displaystyle \sum _{n=1}^{N}a\left(\Delta {t^{\prime }}_m\right)\cdot P_{i,n}\left(t_m-\Delta t_m\right)}\right\}$$E10

(10)

where

$$P_{i,n}\left(t_m-\Delta t_m\right)=H_i\cdot B\cdot {{\sigma }^{″}}_{i,n}\left(t_m-\Delta t_m\right)$$E11

(11)

1. Strain due to bending moment

The incremental strain due to the bending moment in each layer can be similarly expressed by Eq. (13);$${\overline{z}}_i$$ is the distance measured from the neutral plane of the i-th layer.

$$\Delta M_i\left(t_m\right)=\frac{H^3{}_i\cdot B}{12{\overline{z}}_i}\cdot \Delta {{\sigma }^{‴}}_i\left(t_m\right)$$E12

(12)

$$\Delta {{\epsilon }^{‴}}_i\left(t_m\right)=\frac{12{\overline{z}}_i}{H^3{}_i\cdot B\cdot A_i}\left\{\Delta M_i\left(t_m\right)+{\displaystyle \sum _{n=1}^{N}a\left(\Delta {t^{\prime }}_m\right)\cdot M_{i,n}\left(t_m-\Delta t_m\right)}\right\}$$E13

(13)

where

$$M_{i,n}\left(t_m-\Delta t_m\right)=\frac{H^3{}_i\cdot B}{12{\overline{z}}_i}\cdot {{\sigma }^{‴}}_{i,n}\left(t_m-\Delta t_m\right)$$E14

(14)

The strain increment, Δε’’’_i(t_m), can be represented as Eq. (15).

$$\Delta {{\epsilon }^{‴}}_i\left(t_m\right)={\overline{z}}_i\cdot \Delta C_i\left(t_m\right)$$E15

(15)

where ΔC_i(t_m) is the increment of the curvature,

$$C_i(t_m)=\frac{1}{R_i(t_m)}$$E16

(16)

3.2. Strain continuity and equilibrium equation

The global in-plane strain generated at each layer is expressed as the sum of the above-mentioned strain components a) through c). As this global strain must be kept continuous across the interface on each layer, the following expression can be written:

$$\begin{array}{l}{\alpha }_i\Delta T\left(t_m\right)+\frac{1}{H_i\cdot B\cdot A_i}\left\{\Delta P_i\left(t_m\right)+{\displaystyle \sum _{n=1}^{N}a\left(\Delta {t^{\prime }}_m\right)\cdot P_{i,n}\left(t_m-\Delta t_m\right)}\right\}-\frac{H_i}{2}\cdot \Delta C_i\left(t_m\right)\\ ={\alpha }_{i+1}\Delta T\left(t_m\right)+\frac{1}{H_{i+1}\cdot B\cdot A_{i+1}}\left\{\Delta P_{i+1}\left(t_m\right)+{\displaystyle \sum _{n=1}^{N}a\left(\Delta {t^{\prime }}_m\right)\cdot P_{i+1,n}\left(t_m-\Delta t_m\right)}\right\}\\ +\frac{H_{i+1}}{2}\cdot \Delta C_{i+1}\left(t_m\right)i=1_{}t{o}_{}I-1\end{array}$$E17

(17)

The applied load must also satisfy the following equilibrium equations:

$${\displaystyle \sum _{i=1}^I\Delta P_i\left(t_m\right)=0}$$E18

(18)

$$\begin{array}{l}{\displaystyle \sum _{i=1}^I\Delta M_i\left(t_m\right)+\Delta P_1\left(t_m\right)\cdot \left({\overline{z}}^{\prime }-\frac{H_1}{2}\right)+\Delta P_2\left(t_m\right)\cdot \left({\overline{z}}^{\prime }-H_1-\frac{H_2}{2}\right)}\\ +\cdots +\Delta P_I\left(t_m\right)\cdot \left({\overline{z}}^{\prime }-H_1-H_2-\cdots -\frac{H_I}{2}\right)=0\end{array}$$E19

(19)

Assuming the thickness of each layer is negligibly small compared with its resultant curvature, the respective curvature on each layer is identical as Eq. (20). Therefore ΔP_i (t_m), ΔM_i(t_m) and ΔC_i(t_m) are obtainable by simultaneously resolving Eq. (17), Eq. (18), and Eq. (19).

$$\Delta C_1(t_m)\approx \Delta C_2(t_m)\approx \cdots \approx \Delta C_I(t_m)=\Delta C(t_m)$$E20

(20)

3.3. Derivation of stress and deflection

By way of the above-mentioned procedure, in-plane forces and curvatures are obtained. Using these results, stress and deflection can be derived following the steps below. Firstly, the stress generated in each layer is calculated as the sum of the in-plane stress and the bending stress, as follows:

$$\sigma {}_i(t_m)=\frac{P{}_i(t_m)}{H{}_{i}B}+\frac{12z}{H_i^{3}B}{M}_i(t_m)$$E21

(21)

As this model is assumed to be isotropic and unconstrained, the curvature must be identical in all directions within the xy plane. As per Fig. 3, L is taken as the distance from the center to the corner of each layer; the cross-section of each layer is assumed to have the shape shown in Fig. 4, and θ is taken as the slope at the tip of each layer. Accordingly, the maximum tip deflection, δ(t_m), is obtained as follows:

$$\delta \left(t_m\right)=\left(1-\cos \theta \right)/C\left(t_m\right)$$E22

(22)

Assuming θ is relatively small and in-plane elongation is also small, resulting in L = L’, the following expression can be obtained:

$$\delta \left(t_m\right)=\frac{L^2\cdot C\left(t_m\right)}{2}$$E23

(23)

Figure 3.

Definition of L

Figure 4.

Deformation of the cross section

4. Measurement of viscoelastic material properties

Aside from the fact that thermal stress analysis taking viscoelasticity into account has not been adequately carried out to date, one must recognize that it is difficult, in practice, to obtain sufficient accuracy in such experimental measurements. Accordingly, any data measured in such experiments is not accurate enough to relate to the viscoelastic numerical model. Viscoelastic materials exhibit very sensitive temperature dependency, particularly in the vicinity of the glass transition point. Thus, intricate temperature control is required throughout the entire duration of the measurement operation. In addition, the time-domain constants obtained usually span a wide range of digits, in the range of 20 to 30. This leads to the additional task of determining desirable factors, which may prove difficult unless advanced optimization techniques are applied.

This paper also presents a test case for obtaining the characteristics of epoxy resin material. A device for measuring dynamic viscoelasticity, RSA III (TA Instruments), was used. Dynamic viscoelasticity characteristics were measured for angular velocities 3.16, 10, 31.6, and 100 rad/sec with an ascending temperature rate of 2 °C/min, in the temperature range of -40 to 60 °C. As shown in Fig. 5 and Fig. 6, the storage modules, E’, and the loss modulus, E’’, of the epoxy resin material were measured. This measurement device is equipped with a temperature-controlled oven with a solid structure and a large volume flow rate, providing excellent performance in temperature control.

Figure 7 shows the master curve, representing the relationship between the storage/loss modulus and reduced angular velocity, for the reference temperature, T_R. To create the master curve, the WLF formula as shown below was applied for the temperature-time reduction.

$$\log \left(a\left(T\right)\right)=-\frac{C_1\left(T-T_R\right)}{C_2+\left(T-T_R\right)}$$E24

(24)

$$T_R=T_g+50\begin{array}{cc}& C_1=8.86\begin{array}{cc}& C_2=101.6\end{array}\end{array}{T}_g=-{16}^{\circ }\text{C}$$E25

In this study, a simple optimization technique was employed in which generally recommended values for C₁ and C₂ were used, and T_g was estimated using the quasi-Newton method (Kobayashi, 2008). As a result, a single smooth master curve could be used to cover the wide range of angular velocities. The coefficients for the Maxwell model were obtained from the master curve by the optimization approach so that the relaxation curve from Eq. (1) could be calculated numerically. The result is shown in Fig. 8. The relationship between time-domain constants and frequency-domain constants is expressed by Eq. (25). Using this formula, the numerical model that gives a close approximation of the actual measurement data can be obtained, as seen in Fig. 8.

$$E^{\prime }=E_e+{\displaystyle \sum _{n=1}^N\frac{{\tau }_n^2{{\omega }^{\prime }}^2}{1+{\tau }_n^2{{\omega }^{\prime }}^2}{E}_n}{E}^{″}={\displaystyle \sum _{n=1}^N\frac{{\tau }_n{\omega }^{\prime }}{1+{\tau }_n^2{{\omega }^{\prime }}^2}{E}_n}$$E26

(25)

These procedures for identification have been organized in Excel; please refer to the Appendix.

Figure 5.

Storage modulus of epoxy resin

Figure 6.

Loss modulus of epoxy resin

Figure 7.

Master curve for Storage/Loss modulus

Figure 8.

Master curve for relaxation modulus

5. Numerical verification by FEM

To verify the calculation method proposed in this study, a three-layered plate—40 mm in length and 20 mm in width (Fig. 9), consisting of epoxy resin with 0.5 mm thickness (viscoelastic, T_g = 105 °C), FR-4 substrate with 0.5 mm thickness (viscoelastic, T_g = 120 °C), and aluminum with 0.1 mm thickness (elastic)—was analyzed. The initial temperature of the laminated plate was set at 180 °C, and it was cooled to 25 °C within 2,000 sec, as shown in the temperature history (Fig. 10). Because the laminated plate was gradually cooled, its inside temperature was assumed to be uniform.

Figure 11 shows the master curves of the relaxation moduli of the epoxy and FR-4 substrates. The temperature-dependent viscoelasticity was assumed to be in accordance with the Arrhenius formula, as expressed by Eq. (26).

$$\log \left(a\left(T\right)\right)=-\frac{\Delta H}{2.303R}\left(\frac{1}{T}-\frac{1}{T_0}\right)$$E27

(26)

Figure 9.

Three-layered plate model

The values of the coefficients for Epoxy resin shown in Eq. (26) were as follows: ΔH₁ = 2.4420 × 10⁵ (before T₀); ΔH₂ = 3.0480 × 10⁵ (after T₀); and T₀ = 97.7 °C. Those for the FR-4 substrate were as follows: ΔH₁ = 8.8300 × 10⁴ (before T₀); ΔH₂ = 4.3220 × 10⁵ (after T₀); and T₀ = 113.2 °C.

Figure 10.

Applied temperature history

Figure 11.

Relaxation modulus

Please note that the material constants of the epoxy resin used for this verification analysis differ from those expressed in the measurement example in the preceding section. Figure 12 shows the thermal expansion coefficients of the epoxy and FR-4 substrates. The aluminum that composed Layer 3 was an elastic body with the Young’s modulus E = 70,000 MPa and the thermal expansion coefficient α ₃ = 23.2 × 10^-6/°C.

For comparison with the multilayered plate theory developed in this study, an FE analysis of the same model was performed using the shell elements of Abaqus ver. 6.8. Figure 13 shows the deformation of the plate obtained by the FE analysis. Figure 14 shows the relationship between the deformation and temperature. Figure 15 shows the relationship between stress and temperature for each layer. As shown in Figs. 14 and 15, the deformation and stress of the laminated plate obtained by the multilayered plate theory were confirmed to be identical to those obtained by the FE analysis.

Figure 12.

Thermal Expansion curves

Figure 13.

Deformation of plate (FEM)

Figure 14.

History of plate deflection

Figure 15.

History of stress

6. Conclusions

We were able to establish a simplified and easy-to-use tool for estimating the warpage in printed circuit boards based on the multilayered plate theory combined with the effects from the temperature-dependent thermal expansion coefficient and the temperature-dependent viscoelastic characteristics of the resin material. The results derived from this method are confirmed to be in agreement with the FE analysis results.

Appendix (Identification of the generalized Maxwell model and development of a curve fit program)

When the generalized Maxwell model for time domain is identified based on the master curve in frequency domain using Eq. (25), the following three points should be noted.

E_e, E_n, τ_n of positive definite

Since the generalized Maxwell model is regarded as a mechanical model, it is preferable that all the values of these coefficients are always positive. However, the input rule for them is different among general purpose FEMs. For example, some codes allow negative value input, while others strictly prohibit negative value input. Hence, there is no unified rule among all the codes. Since it is empirically observed that master curves may oscillate due to the affect of terms with negative values, it is considered reasonable to control the input data to make them positive definite.

Number of terms in the generalized Maxwell model

It is common practice to give the abscissas of a master curve, i.e., frequencies using a logarithmic scale covering a range of 20 to 30 digits. To make this master curve approximate a smooth curve, it is said that the number of terms (number of two-element Maxwell models) should be selected so they are equal to or above the number of digits of the frequencies. In order to confirm this, a simple calculation was carried out using a single two-element Maxwell model. The calculation was performed under the following conditions.

Elastic model: E = 100 Pa

Viscoelastic model: E₁ = 100 Pa, τ₁ = 1 sec

The relaxation modulus calculated using this model is shown in Fig. 16. In the figure, the curve with the solid line in the relaxation module is noted to decay one digit over time. It shows that a single Maxwell model is capable of representing relaxation behavior over a time domain of about one digit. Accordingly, when each dashpot is provided with the sequence of τ_n such that the next one has an order of time greater by one digit than the former, the relaxation behavior over the full range of a time domain can be expressed without a break. If the abscissa of a master curve, for example, is represented with a time domain of 20 digits, the number of terms in the generalized Maxwell model may be selected to include 20 or more.

Figure 16.

Relaxation behavior of single two-element Maxwell model

Smoothness of relaxation spectra

The next task is to organize the model so that the contribution from each term is approximately smoothed. In accordance with the knowledge derived by Emri et al. (1993), keeping the smoothness of discrete relaxation spectra is effective in securing the desirable accuracy of approximation results. The Kronecker’s delta in the expression is denoted by δ.

$$H_n(\tau )={\displaystyle \sum _{n=1}^N{E}_n{\tau }_n\delta ({\tau }_n-\tau )}$$E28

(27)

An example of these relaxation spectra is shown in Fig. 17. An attempt was made in this example such that the envelope for these discrete spectra is approximated to be piecewise quadratic so that the smoothness can be maintained subject to the curvature change along this envelope being not too large. Through the testing of such provisions, followed by an approximate calculation, it becomes possible to perform a curve fit operation for a master curve even though data is missing. For viscoelastic materials with sharp temperature dependency, it becomes very difficult for the temperature control in the measurement device to catch up with the actual material response, and as a result, critical defects are bound to occur (Fig. 18(b)); therefore, smoothing manipulation for those relaxation spectra is a highly effective measure.

In the curve fit program developed by the author’s company, the generalized Maxwell model is identified based on the master curve shown in Fig. 7. This program is designed to completely fulfil the constraint conditions discussed in the preceding section. A sample output from this program is shown in Fig. 19. The user is only required to enter ”Input data,” “Number of Prony terms,” and “Poisson’s Ratio” in the specified input field, and then press the “Optimization” button. The program automatically performs an approximate calculation. The optimization operation uses the quasi-Newton method. For the quasi-Newton method, it is necessary to set up the initial condition in the vicinity of the optimized value. However, this program incorporates an algorithm that can automatically estimate, from the test results, an initial condition that easily converges.

The master curve for epoxy resin material shown in Fig. 8 covers a range of about 12 digits in terms of angular frequency, so the number of terms in the Maxwell model is set to 12. As previously mentioned, positive values for all of E_e, E_n, and τ_n are maintained during the calculation. The coefficients of the identified Maxwell model are automatically written into the respective format of the input file for Abaqus, Marc, or LS-DYNA.

Figure 17.

Smoothing manipulation for relaxation spectra

Figure 18.

Examples of measured master curve