A Novel Moisture Diffusion Modeling Approach Using Finite Element Analysis

In this study, a novel wetness and moisture concentration analysis approach is presented. A finite element method is utilized for the solution technique mainly using thermal and surface effect elements. Numerical results obtained from the current approach are compared against other existing finite element-based solutions and the newly introduced peridynamics theory. For numerical analysis, a reflow soldering stage is simulated for a multi-material system with time-dependent saturated moisture concentrations. Different solubility activation energies and temperature conditions are considered. Numerical results demonstrate that the developed methodology can make accurate predictions under different conditions and it is more general than some other existing models which are limited to certain conditions.


Introduction
Polymer-based materials such as underfills, molding compounds, etc. are common materials used in microelectronic and optoelectronic components. However, these polymeric components may suffer from a large amount of volume expansion due to moisture absorption. Since electronic packages are multi-material systems, such volume expansions can result in hygroscopic stresses and swelling. Moreover, popcorn cracking can occur at high temperatures during reflow soldering since moisture will transform into vapor with high pressure during this process [1]. Therefore, in order to achieve durable electronic packages, moisture diffusion analysis is required.
The finite element method (FEM) is commonly used for this purpose. Since electronic packages contain multi-material interfaces, moisture concentration is not continuous along these interfaces, which introduces discontinuities. To overcome this problem, a wetness approach [2] was introduced since the wetness field is continuous at the interfaces. If the saturated concentration does not depend on time and temperature, it is possible to make a thermal-moisture analogy and a regular thermal analysis can be performed by making a necessary calibration of its parameters with corresponding moisture diffusion parameters. A review of thermal-moisture analogy schemes can be found in Yoon et al. [3]. According to this review, it was concluded that the direct analogy is limited to single-material systems and the normalized analogy can be used for multi-material systems if thermal loading conditions are isothermal, spatially as well as temporally. As saturated concentration is only a function of relative humidity, regardless of temperature, Jang et al. [4] proposed an advanced thermal-moisture analogy scheme for anisothermal moisture diffusion problems. On the other hand, if the saturated concentration is time-dependent, then thermal-moisture analogy is no longer available. For such cases, piecewise normalization and internal source approaches were developed [5,6]. The piecewise normalization approach is computationally expensive since it requires multiple load steps to ensure accuracy. On the other hand, although the internal source approach is easier to implement, its convergence is reliant on the number of iterations in each time step. More recently, peridynamic theory (PD) [7,8] was utilized to determine moisture and wetness distribution in electronic packages for time-dependent saturated moisture concentration condition. This technique is suitable for treating discontinuities and does not require an iterative solution. Moreover, Chen et al. [9,10] developed a unified activity-based diffusion theory which does not have any interface discontinuity issues since water activity is continuous at the interface. This theory can also be used to unify existing normalized concentration-based theories. In addition, ANSYS, a commercially available finite element software, provides coupled field (CF) and diffusion (DF) elements for moisture diffusion analysis. These elements are suitable for either time-or temperature-dependent saturated moisture concentrations. However, the diffusion element is limited to a uniform temperature field. An in-depth investigation of ANSYS diffusion elements was provided by Liu and Park [11] and Sutton [12].
In this study, a novel wetness and moisture concentration modeling approach is presented for saturated moisture concentration, which is a function of time. This study is an extension of the previously published conference paper [13]. The method uses conventional thermal and surface effect elements and it is implemented in an ANSYS framework. The approach is computationally efficient and straightforward in its computer implementation. Various verification and demonstration studies are considered for multi-material systems having time-dependent saturated moisture concentrations, with different solubility activation energies and temperature conditions. Results from the current approach match well with PD and ANSYS CF element results for time-dependent diffusivity. However, significant differences are observed in the case of temperature-dependent diffusivity. Besides, for different solubility activation energies, there is also a difference between the present results and those obtained using the ANSYS DF element.

Classical Wetness Formulation
The first Fick's law can be used to represent moisture diffusion as where J represents diffusion flux, D is diffusivity, and C(t) is the moisture concentration. The conservation of mass solute during the diffusion process can be satisfied by in which G(t) is the rate of diffusing substance generation per unit volume. Combining Equations (1) and (2), the second Fick's law for homogeneous domains can be obtained as For multi-material systems (nonhomogeneous domains), moisture concentration is discontinuous at multi-material interfaces. Therefore, a normalized wetness parameter w was introduced in Reference [2] as where C sat (t) represents the saturated moisture concentration, which is the maximum concentration that a material can absorb. Note that the wetness parameter is continuous at multi-material interfaces and the condition for chemical potentials to be equivalent is also satisfied at the interface. Therefore, the concentration ratio at the interface of materials 1 and 2 remains constant at any instant as As a result, the wetness field becomes continuous at the interface, i.e., w 1 = w 2 . After substituting the wetness parameter given in Equation (4) into Equation (3), the second Fick's law equation can be rewritten as If the saturated concentration, C sat , is constant, Equation (6) will reduce to Equation (7) can be directly solved by using the thermal-moisture analogy. However, if the saturated concentration is time-dependent, then the thermal-moisture analogy cannot be utilized. For this case, the commercial finite software ANSYS provides CF and DF elements [14], which are based on the modified form of the first Fick's law as and where v represents the transport velocity vector, k is the Boltzmann constant, Q is the particle heat of transport, and T(x, y, z, t) represents temperature. Combining Equations (8) and (2) and using the wetness parameter definition given in Equation (4), the conservation of mass solute equation can be rewritten as If the saturated moisture concentration rate is neglected, Equation (10) can be shortened as A similar derivation can also be conducted for Equation (9) and the conservation of mass solute equation can be derived as If the saturated concentration is temperature-dependent, i.e. the C sat (T) and chain rule are utilized, Equation (12) can take the form of Note that in Equation (13) the diffusivity, D, can be a function of time or temperature.

Peridynamic Wetness Formulation
The moisture diffusion equation given Equation (2) can be expressed in peridynamic form as [7] ∂ ∂t where f w ′ , w, x ′ , x, t is the wetness response function which represents the exchange of moisture between material points x and x ′ that are connected through hygro-bonds, as shown in Figure 1. In PD, material points interact with each other in a nonlocal manner and a material point can interact with other material points within its neighborhood named as the horizon, H (see Figure 1), with a radius of δ.
(1-D) 3 6 (2-D) 4 6 (3-D) 1 ,,,, The wetness response function can be defined as where the PD bond parameter, d(t), in Equation (15) can be expressed in terms of the classical diffusivity, D(t), and the saturated concentration, C sat (t), as In Equation (16), A is the cross-sectional area of a one-dimensional structure and h is the thickness of a two-dimensional structure.
The PD moisture diffusion equation given in Equation (14) can be solved by using the meshless approach. Therefore, the integral term can be replaced by a finite summation and Equation (14) can be rewritten as with where w (k) and w (j) are wetness values at material points x (k) and x (j) , respectively, and N represents the number of material points inside the horizon of the material point x (k) . During the solution process for each time step, the time dependency of C sat can be approximated with the backward Euler method between the consecutive time steps as and an iterative solution is not required to obtain the solution.

ANSYS Coupled Field and Diffusion Element Formulations
As mentioned earlier, the commercial finite element software, ANSYS, provides two different types of elements for moisture diffusion analysis. The first type is diffusion (DF) elements, which are available in different forms including Plane 238 (eight-node two-dimensional element), Solid 239 (20-node three-dimensional element), and Solid 240 (10-node three-dimensional element). These elements have only moisture concentration as the degree of freedom (DOF) per node. However, if the saturated moisture concentration, C sat , is defined with the command "MP, CSAT" in ANSYS, the DOF becomes the wetness parameter instead of the concentration. The formulation of these elements is based on Equation (11) and they allow the saturated moisture concentration to be a function of time or temperature. However, DF elements are only available for a uniform temperature. If the temperature field is not uniform, ANSYS offers coupled field (CF) elements including Plane 223 (eight-node two-dimensional element), Solid 226 (20-node three-dimensional element), and Solid 227 (10-node three-dimensional element). These elements have two DOFs at each node as wetness, w, and temperature, T. The formulation of these elements is based on Equation (13). The finite element approximation of field variables for the element e can be expressed as and where S is the vector of element shape functions, while T e and w e represent the vectors of unknown nodal temperature and wetness values of the element. Utilizing the virtual work principle for the element e with the arbitrary virtual quantities of w e and T e , Equation (13) can be rewritten as where . w e and . T e represent the nodal time derivative of wetness and temperature, respectively. Moreover, C d is the element diffusion damping matrix, C dt is the element thermal-diffusion damping matrix, K d is the element diffusion conductivity matrix, K dN 1 is the nonlinear part of the element diffusion conductivity matrix associated with thermomigration, K dN 2 is the nonlinear part of the element diffusion conductivity matrix due to ∂C sat /∂T, K dtN is the nonlinear part of the element transport conductivity matrix, K dt 1 is the element transport conductivity matrix, K dt 2 is the element thermal-diffusion conductivity matrix due to ∂C sat /∂T, K dtN 2 is the nonlinear part of the element transport conductivity matrix due to ∂C sat /∂T, and R is the nodal diffusion flow rate vector. The explicit definitions of these matrices are given in Reference [14] as where V e is the volume of the element. Utilizing Equation (22), the coupled finite element thermal-diffusion matrix equations can be written as where C t , K t , and Q represent the element specific heat matrix, the element thermal conductivity matrix, and the element heat generation load matrix, respectively. The thermal conductivity matrix, K t , is a combination of the element mass transport conductivity matrix, K tm , the element diffusion conductivity matrix, K tb , and the element convection surface conductivity matrix, K tc . If thermomigration is negligible, Equation (33) can be simplified as

ANSYS Thermal and Surface Effect Element Formulations
In addition to elements provided by ANSYS described in the previous section for moisture diffusion modeling, an alternative approach is presented in this study by utilizing conventional thermal and surface effect elements. Utilizing the virtual work principle in conjunction with the divergence theorem, Fick's second law equation given in Equation (6) can be written as where S e and n represent the surface of the element and unit normal vector to this surface, respectively. In Equation (35), the time derivative of the saturated concentration, i.e., dC sat /dt, can be expressed by using backward Euler method as After substituting the wetness field expression, w, given in Equation (20) into Equation (35) for arbitrary δw e , Equation (35) can be written in a matrix form as where Note that this equation has the same form as the heat flow equation, which can be written as where ρ, c, D(t), h f , and q represent the density, the specific heat, the thermal conductivity, the film coefficient, and the heat generation rate per unit volume, respectively. Since the wetness and heat flow equations given in Equations (37) and (42) have the same form, it is possible to relate the thermal parameters, ρ, c, D t , h f , and q to diffusion parameters C sat (t), DC sat (t), ∂C sat /∂t, and G. Therefore, for one-dimensional analysis, the wetness field equation given in Equation (37) can be made equivalent to the thermal heat flow equation given in Equation (42) by using a combination of thermal link (LINK33) and thermal surface effect (SURF151) elements between two nodes (see Figure 2) and adjusting parameters as demonstrated in Tables 1 and 2. Note that the thermal link element is necessary to construct C t e and K tb e matrices and the thermal surface effect element is necessary to construct the K tc e matrix. In Table 1, the subscript n in D n+1 and C sat,n+1 represents the current time step number. In Table 2, the rate of C sat can be approximated with the backward Euler method as ∆C sat /∆t, in which ∆ implies the difference of values between the time steps. Moreover, the nodal heat generation rate per unit volume, q, can be used in place of the diffusing substance generation rate per unit volume, G, in ANSYS by defining the nodal loads with the "BF, Node, HGEN, Value" command.

Numerical Results
In this section, numerical results for the analysis of the desorption process in a bar during reflow soldering are presented by considering conventional thermal link and surface effect elements, diffusion (DF) and coupled field (CF) elements, and peridynamics. The two-material bar shown in Figure 3 is initially fully saturated at 85 • C/100% RH with isolated lateral surfaces. The material properties for diffusion and thermal analyses are given in Tables 3 and 4 based on the data provided in Reference [15].   Vapor pressure activation energy, E VP (J/mol) 4.08737 × 10 4 4.08737 × 10 4 Table 4. Material properties for materials 1 and 2 for thermal analysis.

Material 1 Material 2
Density, ρ (kg/m 3 ) 3 × 10 3 3 × 10 3 Specific Heat, c (J/kgK) 1.5 × 10 3 1.5 × 10 3 Thermal Conductivity, k (W/mK) 0.2 0.6 The initial and boundary conditions of the wetness field are specified as where 2L and h represent the length and thickness of the bar, respectively. The solubility activation energy of materials 1 and 2 are either equal or unequal at the interface. The time-dependent saturated concentration, C sat (t), and diffusivity, D(t), can be expressed as and where R is the universal gas constant (R = 8.3145 J/Kmol), D 0 is the diffusivity factor, S 0 is the solubility factor, P 0 is the pressure factor, E D is the activation energy of the diffusivity, E S is the activation energy of the solubility, E VP is the vapor pressure activation energy, and RH indicates the relative humidity. For CF elements, the saturated concentration, C sat (T) can be specified as temperature-dependent. On the other hand, the diffusivity can be specified as either time-dependent, D(t), or temperature-dependent, D(T). Moreover, the temperature field can be either uniform or nonuniform. However, the DF element is not suitable for nonuniform temperature conditions since the temperature is not known a priori. In addition, the solutions with diffusion and coupled field elements are obtained by neglecting transport velocity and thermomigration, respectively.
The finite element (FE) discretization of a bar with conventional thermal elements (present approach) and with two-dimensional DF or CF elements are shown in Figures 4 and 5, respectively. The bar has a length of 2L = 2 mm, thickness of h = L/100, and cross-sectional area of A = h 2 . The discretization size is specified as ∆x = h. The total time is specified as t = 80 s with a time step size of ∆t = 2 s. The moisture concentration results are plotted at four different time instances, t = 20, 40, 60, and 80 s. The results obtained from the current approach are compared against peridynamics as well as ANSYS DF and CF element results for uniform and nonuniform temperature conditions with equal and unequal solubility activation energies of materials 1 and 2. In the case of unequal solubility activation energy, material 2 has E S 2 = 4.50 × 10 4 J/mol.  . The finite element discretization of a bar for "diffusion" or "coupled field" elements.

Bar Subjected to Uniform Temperature
In the first case, the two-material bar configuration is subjected to the uniform temperature field of As shown in Figures 6 and 7, the predictions based on the present approach, CF elements with time-dependent diffusivity, D(t), and temperature-dependent saturated concentration, C sat (T), and peridynamics are in excellent agreement for both equal and unequal values of E S under uniform temperature conditions. However, the CF elements with temperature-dependent diffusivity, D(T), and saturated concentration, C sat (T), and the DF elements show significant deviation from the other solutions. When diffusivity is imposed as a time-dependent property, D(t), its value is known before the next time step as part of time integration, and does not require approximation. However, if temperature-dependent diffusivity, D(T), is considered, this requires approximation for the next time step. Therefore, this will yield results significantly different from the correct solution. Moreover, the DF elements do not predict the correct concentration because they disregard the effect of the time derivative of the saturated concentration.

Bar Subjected to Nonuniform Temperature
In the second case, a nonuniform temperature condition is considered. The nonuniform temperature is obtained by applying the temperature boundary conditions in the form of The DF element is not applicable due to the nonuniformity of the temperature field. Moreover, the CF element with time-dependent diffusivity, D(t), is also not applicable because diffusivity, D(t), given in Equation (50), highly depends on temperature which is the unknown DOF in this problem. The concentration results along the bar are plotted in Figures 8 and 9. As can be seen from these figures, the predictions from the present approach and peridynamics are in very good agreement, whereas there is a significant difference between the results from the CF elements with temperature dependence, D(T).

Conclusions
This study presents a new finite element-based modeling approach for moisture diffusion analysis in the presence of time-dependent saturated moisture concentration by utilizing conventional ANSYS thermal and surface effect elements. The capability of the present approach is demonstrated by comparing results from the present approach with those obtained using ANSYS diffusion and coupled field elements as well as peridynamic theory. For comparison purposes, a fully saturated one-dimensional two-material bar configuration is considered with either the same or different solubility activation energies in the presence of a time-dependent saturated moisture concentration. Both uniform and nonuniform temperature conditions are analyzed. For all cases, the current approach is capable of predicting correct concentration values. This is not true for ANSYS diffusion and coupled field elements. If the diffusivity is a function of temperature, coupled field elements could not provide an accurate solution. On the other hand, although the diffusion element is shown to be suitable for a uniform temperature, it can yield incorrect results, especially for unequal values of solubility activation energies in a multi-material system. As a summary, the present approach is superior to the diffusion and coupled field elements of ANSYS, especially under nonuniform temperature conditions.