Prediction of Adhesively Bonded Polyurethane-to-Steel Double Butt Joint Failure Using Uncertainty Analysis

Adhesively bonded joint has gained increasing popularity due toweight reduction and relief of stress concentration.However, damage in the mode of adhesive failure and cohesive failure, as well as the adherend failure, could still occur, and it has been realized that the reliability of the adhesively bonded joint depends on numerous complex and even nonlinearly interacting factors. Consequently, the prediction of load-bearing capacity and damage localization for an adhesively bonded joint can be difficult due to the not well-known effects of the variations or uncertainties in the imperfected adhesive-adherent bonding surfaces or the environmental conditions as well as the material and geometrical nonlinearities. Although abundant uncertainties are present, the standard analysis tool in industries is still deterministic finite element analysis (FEA). ,e routine practice of such analysis is applying a cohesive zone model (CZM) implemented with a proper traction-separation law to predict the onset and gradual degradation of adhesionwhichmay underpredict or overpredict the load-bearing capacity of an adhesively bonded joint. In this study, deterministic FEA with a CZM is first applied to predict the load-displacement curve of an adhesively bonded polyurethane-to-steel double butt joint. A comparison of the prediction with the experiment reveals the inability of a deterministic approach for accurately predicting the load-bearing capability of the joint and the failure propagation route.,en, uncertainty analysis using polynomial chaos expansion (PCE) is applied to examine if it can enhance the prediction of the joint failure. ,e results show the predictions generated by the uncertainty analysis correlate better than the deterministic analysis with the test data, hence demonstrating the potential of uncertainty analysis in improving the prediction of the failure mode and load-bearing capability of an adhesively bonded polyurethane-to-steel double butt joint.


Introduction
Adhesively bonded joints have been widely adopted in numerous industrial applications, e.g., automotive industry, architecture, and aerospace. Compared to the conventional mechanically fastened and welded joint, adhesively bonded joint becomes increasingly popular owing to the greater strength-to-weight ratio, reduction in the stress concentration, more flexibility in the assembly process, and higher resistance to fatigue [1]. For an ideal adhesively bonded joint design, the load-bearing capacity prediction and possible failure mode analysis are important. e failure of an adhesively bonded joint can be generally categorized into three modes [2]: (I) adhesive failure, which occurs in the interface between the adhesive and the adherent, (II) cohesive failure, which means the adhesive itself damages during services, and (III) adherent failure. e initiation and propagation of these three failure modes can be associated with numerous factors. To name but a few, the pretreatment of the adherent surface, the physical and environmental conditions during the adhesive applying process, such as temperature and humidity, etc. More importantly, owing to the indispensable imperfection involved in the adhesive-applying process, especially for bonding dissimilar adherent (e.g., soft porous material bonded with the metal), the actual bonding area between the adhesive and the adherent is only partial of the nominal bonding area even though a spacer is used [3]. In addition, geometrical [4] and material nonlinearity [5] and some complex behaviors such as creeping [6] could be introduced when the joint is subjected to time-variable, longterm history external excitations [7]. erefore, with the abundance of variations and uncertainties in these multiple and possibly nonlinearly interacting factors, an accurate prediction of the failure mode and the load-bearing capability of an adhesively bonded joint can be difficult [2].
Despite being associated with uncertainties, a deterministic approach is still the standard practice in industries. Depending on the situation whether the problem can be simplified as a plane stress problem, the deterministic approach uses either 2D or 3D finite element analysis (FEA) coupled with an adhesive failure criterion [8] or a failure initiation-propagation-separation simulating scheme [5] to predict the strength and failure mode of adhesive joints. e method of using the adhesive failure criterion is based on the theory of continuum mechanics to investigate the stress or strain distribution along the midlayer of the adhesive and perform failure prediction using maximum normal or shear stress and strain [9] or critical longitudinal strain [10] as indicators. Although considerable success has been reported, the limitation of such methods, however, becomes prominent induced by the ductile adhesive problem being tackled and the result's mesh dependency due to stress singularity [11]. e method of CZM (cohesive zone model) has been shown to be more predictive as it accounts for both the crack initiation indicated by continuum mechanics theory and the crack propagation triggered by strain energy release [12]. e failure prediction generated by the CZM has been found to be less mesh-sensitive, and this technique is recommended to be the best practice for the time being. However, as discussed by Ramalho et. al. [2], the success of the CZM can be compromised by inaccuracy in the description of the traction-separation law. Also, the prediction can be further deteriorated by improper placement of cohesive elements along the potential failure propagation path. e potential failure propagation path is not always well known as a prior especially when (i) the strength of the adherend is not significantly stronger than the adhesive and (ii) the adherend has a porous cellular microscopic structure, which leads to partial bonding. One example is the adhesively bonded polyurethane foam-to-steel joint, which has been widely applied in the reefer trailer and thermal-insulating architecture industries. Hence, for such cases, an uncertainty analysis with incorporation of stochasticity in the bonding area is worthy of attention.
Recently, a few research studies have been aimed at applying uncertainty analyses to the adhesive strength evaluation for both single-lap and double-lap CFRP-to-steel joints and the pull-out force prediction of a FRP-to-concrete bonded joint [13,14]. However, the problems are all simplified as analytical models, and it has been pointed out by Ramalho et al. [2] that the analytical model inevitably involves simplifications and assumptions limiting its applicability. One difficulty in applying uncertainty analysis coupled with numerical models is the large computational cost because sampling-based uncertainty analysis often requires tens of thousands runs to achieve convergent statistical estimations. As the degrees of freedom of a stochastic problem (the number of random variables or random processes) increase, such uncertainty analysis becomes soon implausible. Tipireddy and Kumar pioneered in applying nonsampling-based uncertainty analysis approach, e.g., polynomial chaos expansions (PCE), to estimate failuresafety frontier of an adhesive anchor [15]. Although still analytically treated the problem, it shed some light for the possibility of coupling numerical models with uncertainty analysis for the adhesive strength estimation by reducing the computational resources to an acceptable level. e objective of the present work is, therefore, to explore the potential of an uncertainty analysis coupled with numerical models in improving the prediction of load-bearing capacity of an adhesively bonded joint. For this purpose, experiments are carried out to determine the mechanical properties of the adherent and the adhesive and the load-bearing capacity of the joint. en, a load-bearing capacity prediction generated by a deterministic model implemented with a mix-mode CZM is compared to the measurement results for demonstration of the limitation of the deterministic approach. By treating the debonding area between the adherent and the adhesive as a random variable, which has been found to be a dominant factor affecting the adhesion strength by Sorrentino et al. and Lai et al. [3,16], an uncertainty analysis incorporated with polynomial chaos expansion is applied to quantify the variation of the loadbearing capacity prediction. Finally, the potential of the present approach in enhancing the adhesive strength prediction is discussed.

Experiments
As depicted in Figure 1, an adhesively bonded polyurethaneto-steel double butt joint, which is a popular design for reefer trailers, was fabricated as the test specimen. It consists of two outer steel covering plates bonded with two inner boards by an adhesive with a thickness of 0.3 mm. e inner board is a composite of two 0.6 mm thin layers of steel sheets and a thick layer of polyurethane thermal-insulating foam with a thickness of 78.8 mm.
e width and the length of the specimen are 180 mm and 2250 mm. e bonding between the steel layer and the polyurethane foam layer is also achieved by adhesives with a thickness of 0.3 mm. e bonding process consisted of polishing, debris cleaning, and adhesive application, as well as pressurization for solidification in sequence.
In order to simulate the failure mode and predict the load-bearing capacity of the joint, the mechanical properties of all the components, i.e., the steel, the polyurethane foam, and the adhesives, need to be experimentally determined, which are demonstrated in Sections 2.1∼2.3 followed by a three-point bending to measure the bending moment bearing capacity of the joint.

Uniaxial Tension of Steel.
e adherent is a type of prepainted organically coated steel sheet. It consists of three parts, namely, organic coating layer, metallic coating layers, and the steel core. is type of steel sheet is particularly popular in the reefer trailer industry as it features high quality of visual appearance and high resistance to abrasion and weathering.
ree dog-bone-shaped specimens made of steel were fabricated using precision laser cutting and uniaxial tension were performed with displacement-controlled at a rate of 1 mm/min using an electromechanical testing machine. e specimen dimension and the measured engineering σ-ε curves are displayed in Figure 2 showing the measured σ-ε curves of the three specimens are fairly consistent, and Young's modulus is extracted as 212 ± 5 GPa using simple linear curve fitting. Poisson's ratio is taken from the literature as 0.3 [17]. e yielding and the ultimate strength are 365 ± 1 and 526 ± 1 MPa, respectively.

Uniaxial Compression and Tension of Polyurethane and the Hyperelastic Material Model Characterization.
e polyurethane foam studied in this work is provided by BASF ® company with the core density being 40 kg/m3. e polyurethane foam is a common practice in the thermal insulation application because of the low heat transmission rate (typically around 0.021 W/m·k).
Polyurethane foam may undergo large volume change subjected to external loading, and the stress-strain relationship is nonlinear. Hyperfoam model is a hyperelastic material model, which has been shown to be suitable for polyurethane foam characterization. Without considering the thermal effect, the strain energy of the hyperfoam model is written as [18] where μ i , α i , and β i are the parameters to be identified, λ 1 , λ 2 , and λ 3 are the three principle stretches, and J is the volume strain with J � λ 1 λ 2 λ 3 ≠ 1, owing to polyurethane foam being compressible. For the uniaxial mode (including uniaxial tension and uniaxial compression), the following relations hold [18]: in which ε U is the uniaxial strain so that the strain energy can be simplified as and the axial stress can be derived as e hyperfoam parameters, μ i , α i , and β i , can be extracted by nonlinear least square fit with correlating the experimental uniaxial test data and equations (2)∼(4). e uniaxial compression and tension test setup are shown in Figures 3 and 4, respectively. e experiments were carried out using a universal mechanical testing machine with the strain rate being 2 mm/min for both compression and tension tests. For each type of test, three specimens were prepared, and the measured strain-stress curves are shown in Figure5. As depicted in Figure 5(a), the compression of the polyurethane foam shows a clear nonlinear character with three distinct phases: elasticity (strain range 0%-8%), plateau (strain range 8%-50%), and densification (strain range 50%-90%). is nonlinear phenomenon is induced by the microscopic cellular structure of the polyurethane foam. Within the small strain range, the foam behaves elastically due to cell wall bending. Next, a plateau of deformation occurs due to cell wall buckling. Finally, the cell walls crush together causing the densification phase and leads to the rapid increase of stress. ough these three phases generally exist for all polyurethane foams, the exact strain range of the sequential phases could significantly differ with the foaming formula, foam density, environmental condition, and pressurization during the foaming process. One of the main reasons is these aforementioned conditions could influence the closure rate of the cells of the foam. Furthermore, if the polyurethane foam is adhesively bonded to a metal adherent, the actual adhesion area could be largely uncertain owing to the foam cells not being fully closed. Justification of this point can be well provided by an illustration of closed and open polyurethane foam presented in Figure 6.
Using the uniaxial test data with the aid of nonlinear least square fit, the parameters of the hyperfoam material model are identified and listed in Table 1. To verify the identified parameters, μ i , α i , and β i are substituted into the hyperfoam model, and the compression and tension are simulated with Altair HyperWorks. e models for simulating compression and tension of the polyurethane foam specimen are presented in Figures 5(a) and 5(b). For both models, the polyurethane foam is discretized with 8-node linear hexahedral elements, and the material model is defined as the hyperfoam model with the aid of identified μ i , α i , and β i . e loading and the constraint are consistent with the experimental setup, and then, the engineering strain and stress are calculated by the implicit static solver. It is noticeable that a thin FRP (fiberreinforced polymer) sheet covers the end of the tensile polyurethane foam for both the experiment and the simulation.
is is for protecting the foam for avoidance of premature damage induced by the clamping force. Figures 5(c) and 5(d) compare the simulated and the measured stress-strain curves, and good agreement has been achieved.

Characterization of the Adhesive Material.
e adhesive material is supplied by H.B. Fuller ® company with a product name of Kӧrapur 840. is two-component semistructural adhesive is used for bonding the outer steel plate with the inner boards (refer to Figure 1). Of the inner boards, bonding between the steel layer and the polyurethane foam was also achieved by the adhesive, and thanks to the provision of the mechanical properties by the supplier, Table 2 lists the data for modelling the behavior of the adhesive. Figure 7 shows the experimental setup of the three-point bending test. e specimen was simply supported by two steel I-beams, and bolts were used to constrain the I-beams to the ground. Four steel blocks are used to constraint the lateral movement of the specimen. e effective span of the specimen is 1500 mm, and its parallelism with the ground is ensured by a laser level meter. e force applied to the specimen was exerted by a hydraulic actuator placed at the middle of the specimen and controlled by a MTS electrodynamic test system. e crack initiation-growth-failure process was recorded using a camera.

Numerical
Model. e deterministic prediction of the load-bearing capacity of the adhesively bonded joint is generated by FE analysis using HyperWorks® package. Figure 8 shows the mesh of the joint. e adherends, namely, steel cover plate and the steel sheet, are discretized with 4node linear plane-strain elements due to their thicknesses being significantly smaller than other dimensions. e polyurethane foam is discretized with 8-node linear hexahedral elements. e mechanical behavior for the steel cover plate and steel sheet is assumed to be elastically isotropic for simplicity, while a hyperfoam constitutive model is used for describing the polyurethane foam with the aid of the least square data fit results provided in Section 2.2. e adhesive layers are modelled with single-layer 4-node cohesive elements, and they share nodes with the surrounding elements of the adherends. It has been pointed out that the mesh around the cohesion area needs to be fine enough especially near the overlapping tips as, theoretically, stress singularity occurs here. erefore, mesh biasing effects are needed to denser the adhesion region based on [19] l e � 32Nt 2 9πEG c , (5) in which l e is the optimal element size, N is the element number within the length of the cohesive zone, t is the maximum traction, and E and G c are Young's moduli and critical energy release rate. e interaction between the joint and the two I-beams is modelled with frictionless finite sliding surface-to-surface contact, while the lower surface of Advances in Materials Science and Engineering the I-beams is fully restraint. e lateral movement of the joint is constrained by the L-type stoppers. e actuator is discretized with linear hexahedral elements, and its pressurization exerting on the joint is achieved through defining frictionless finite sliding surface-to-surface contact pair.
A CZM is adopted to describe the adhesive tractionseparation behavior. For a pure tension (shear) mode, the adhesion initially behaves elastically until the cohesive strength is reached (denoted in Figure 9 by t 0 n and t 0 s for tension and shear strength, respectively, in which the   Advances in Materials Science and Engineering subscripts n and s represent the normal and the shear modes). From this point on, damage initiates, and the cohesive stiffness degrades up to complete failure which is denoted by δ f n and δ f s [20]. Since the double butt joint studied in this work underwent three-point bending, shear stress and normal stress coexist, thereby necessitating adoption of the mix-mode CZM. e constitute equations governing the initial elastic phase relating the stress t n(s) and the strain ε n(s) in the mix mode can be written as A reasonable approximation is provided by Campilho et al. [21] as E nn � E, E ss � G, and E ns � 0 for thin layers of the adhesive. E and G are the tensile and shear modulus of the adhesive. en, transit from the elastic phase to the stiffness degradation phase occurs when the following quadratic strain criterion meets [21]: where is the Macaulay operator indicating compression does not induce any adhesive damage. ε 0 n is the initial tension damage strain, and ε 0 s is the initial shear damage strain. e damage propagation is described as softening in the adhesion stiffness as in which t m is the stress under the mix mode predicted by elastic characteristics without damage. Δ 0 and Δ max are the effective displacement at damage initiation and the maximum effective displacement during loading history, respectively. Δ f is the failure displacement under the mix mode. By assuming a classic triangular (bilinear) CZM [20][21][22], Δ f can be related to the total energy release during the adhesive separation G as [20] G � 1 in which G n and G s are the energy release for pure tension and shear mode, respectively. m is the proportion of tension in the mode mixture, and it is computed based on decomposition of the current state of deformation (also termed as nonaccumulative measure of energy). α is set as unity in the present study as recommended by Campilho et al. [21]. In order to better illustrate how we determine the energy release rate of the adhesive, fracture tests were performed. A DCB (double cantilever beam) test (depicted in Figure 10(a)) was conducted to estimate the mode-I energy release G n . According to the ASTM 5528 standard based on MBT (modified beam theory) [23], it has been pointed out by Khoshravan and Asgari Mehrabadi [24] that the original MBT [23] needs to be refined because G n can be overestimated due to imperfections and tolerances in the beam fabrication. erefore, a refined MBT formula with a correcting factor |Δ| is adopted here to estimate G n as [24] G n � 3Pδ 2b(a +|Δ|) , (10) where P is the applied load, δ is the crack's opening displacement, b is the specimen width, a is the crack length as shown in Figure 10(a), while the correcting factor |Δ| is related to P and δ. Let C � δ/P and C 1/3 be linearly dependent on a as

Advances in Materials Science and Engineering
A and B are two unknowns, and they can be determined by least square fit because C and a have already been measured. |Δ| equals to −B/A so that G n can be determined using equation (10).
An ENF (end-notched flexure) test was conducted to estimate the mode-II energy release G s [25], as shown in Figure 10(b). G s can be calculated as follows: where P, C (�δ/P), b, L, and a s are the load applied to the specimen, loading point compliance, specimen width, half of the bending span, and critical crack length, respectively.

Comparison of Deterministic Prediction with Experimental
Results. e force-displacement curve is calculated by incrementally increasing the vertical displacement of the actuator and exporting the reaction force of postprocessing the simulation data. e calculation result is compared to the measured data in Figure 11, showing the deterministic model overpredicts the strength and the stiffness. is overprediction can be illustrated by Figure 12(a) in which the calculated damage indicator is zero. Figure 12(b) shows the value of cohesive opening, which is defined as ε � 〈ε n 〉 2 + ε 2 s . Because ε 0 n � ε 0 s � 0.3, the threshold of cohesive opening yields ε � 0.09. It is shown in Figure 12(b) that the maximum cohesive opening is 0.069 so that the adhesive in not predicted by the deterministic model to be failure.

Advances in Materials Science and Engineering
However, the initiation of adhesive failure is clearly experimentally evidenced in Figure 12(c) that the steel sheet starts losing bonding with the polyurethane foam and propagating to a complete steel-foam interface debonding failure (Figure 12(d)). Hence, the deterministic model is unable to accurately predict the adhesive damage occurrence so that the strength of the joint is consequently overpredicted.

Polynomial Chaos Expansion.
Polynomial chaos expansion (PCE) is a nonsampling method for solving a stochastic problem. Unlike the Monte Carlo method which relies on repeatedly sampling from a distribution and the convergency rate which is only 1/ � n √ (n is the number of samples), PCE has been shown to be convergent in a L 2 sense (quadratically) for any second-order stochastic process, thereby making it suitable to be coupled with FE codes [26]. Let X t; ξ { } be a stochastic process, and it can be expanded as [26] where ξ is a vector of independent variables with known probability density function, and ξ is assumed to be normally distributed in this study. ϕ i is a family of orthogonal polynomials (e.g., Hermite polynomials) with the following relationship: where δ ij is Kronecker delta, 〈., .〉 is the ensemble average, and Ω is the compact support of the probability space   defined by ξ. Since ξ and ϕ i have been predefined by the user, the expansion coefficient, i.e., x i (t) in equation (13), is the only unknown to be determined. ere are two methods for calculating the expansion coefficient, namely, intrusive method and nonintrusive method [27]. e intrusive method involves substitution of equation (13) into the governing equation of the treated stochastic problem and uses Galerkin projection to project the resultant equations onto the random space spanned by the polynomial basis. is is achieved with the aid of orthogonality illustrated by equation (14), and a set of coupled equations can be consequently derived, and the solutions are the expansion coefficient. Although being mathematically rigorous, this method however becomes prohibitive when a large-scale finite element problem is treated because the Galerkin projection onto a commercial FE code is not generally accessible. e nonintrusive method, also named as the stochastic collocation method, is more convenient to use. It can use a commercial deterministic FE code as a black box and statistically estimate the expansion coefficient by running the FE code with randomly generated samples. For example, let Y be a function of two i.i.d (independently identically distributed) normal random variables ξ 1 and ξ 2 so that Y can be expanded using second-order truncated Hermite PCE as ere are six unknowns to be determined, i.e., x 0 ∼x 5 ; therefore, at least six collocation points need to be sampled from a standard normal distribution. However, it has been recommended by Huang et al. [26] that the collocation points need to be at least twice the number of unknowns, and the roots of the higher-order Hermite polynomials are the ideal candidates for the collocation points. Substitution of the collocation points into equation (15) yields Y T � NX T with N being an n × 6 matrix with each row given by where n and i represent the number and the index of the collocation points. Since the collocation points outnumber the unknowns, Y T � NX T is an overdetermined set of equations, and the approximate solutions can be generated by the least square or regression method.

Numerical Predictions Generated by Uncertainty Analysis.
As discussed in Section 2.2, the actual bonding area between the polyurethane foam and the steel sheet must only be partial of the nominal bonding area due to the cell of the foam being not fully closed. To account for this uncertainty, the cellular closure rate of the polyurethane foam is measured, thanks to the assistance provided by BASF Co., Ltd. e results show the closure rate is highly variable with a mean value of 10% of the total foam volume and a standard deviation of 5%. Meanwhile, the exact debonding position is also difficult to be localized due to the scale of the cell being too small (roughly 100 μm measured by a SEM [28]) so that the hollow cell would be too many to model with FE technique. erefore, the proportion of the debonding area to the nominal bonding area between the polyurethane foam and the steel sheet is modelled as a normally distributed random variable by random removal of the cohesive elements in the uncertainty analysis.
It can be found in the statistical estimation convergence study presented in Appendix that the 10th order of polynomial chaos expansions would suffice. According to the discussion in Section 4.1, at least 20 collocation points are needed to calculate the expansion coefficient. Figure 13 shows the median displacement-force curve calculated by uncertainty analysis correlates much better with the measurement than the deterministic prediction. It seems random removal of cohesive elements softens the stiffness of the joint making the uncertainty analysis better reflect inherent debonding between the steel sheet and the polyurethane foam induced by the polyurethane's hollow cellular structure. However, for the deterministic model, it assumes that bonding between the adherents is perfectly established. erefore, the overprediction of the joint strength is unavoidable. Figure 14 displays the results of one representative realization of the uncertainty analysis. Compared to the deterministic prediction shown in Figure 12, the damage indicator and the damage opening value of the uncertainty analysis (Figures 14(a) and 14(b)) clearly suggest the adhesive failure will occur along the edge of the actuatorspecimen contact interface. e photograph of the adhesive failure initiation and propagation recorded in the measurement (Figures 14(c) and 14(d)) proves the correctness of the failure prediction. Hence, the enhancement of predictivity of the uncertainty analysis has been demonstrated by accurately predicting both the load-bearing capacity and the failure localization.

Conclusion
While coupled with abundant uncertainties arising from adhesive application process, adherent surface roughness, temperature, and humidity, the numerical analysis of an adhesively bonded joint is still largely treated as a deterministic problem.
In this study, first by applying a deterministic approach to predict the bending moment bearing capacity of an adhesively bonded polyurethane-to-steel double butt joint, we show the deterministic approach overpredicts the strength of the joint. e predicted stiffness is prominently higher than the measured value. is discrepancy can be attributed to the polyurethane foam being not fully closed which leads to the actual bonding area between the polyurethane foam and the other steel adherent which is only partial of the nominal bonding area.
is partial bonding would inevitably deteriorate the quality of adhesion causing a decrease in the actual joint strength. e disadvantage of the conventional deterministic model is also illustrated as it does not successfully capture the adhesive failure. To explore an effective modelling technique, we proposed a methodology for incorporating stochasticity in the bonding interface into modelling. An uncertainty analysis with the bonding area between the polyurethane foam and the steel adherent treated as a random variable using polynomial chaos expansion is conducted to test its potential on improving the strength prediction of the joint. Compared to the deterministic results, the joint's strength and stiffness prediction generated by the uncertainty analysis have been significantly enhanced.
e force-displacement curve shows better correlation with measurement, while the damage localization and failure mode are ideally predicted by the uncertainty analysis in the meantime.
e present study suggests that it would be beneficial to account for stochasticity if the adherent has a porous cellular microscopic structure, e.g., polyurethane foam, because the bonding could be imperfect. From the prospective of design guidelines, since the deterministic approach overpredicts the strength and stiffness of the joint, the uncertainty analysis could also serve as a better performance evaluation tool. e significance and the applicability of the proposed approach to a larger-scale adhesively bonded structure will be further investigated in a future study.

Convergence Study of the Order of Polynomial Chaos Expansion
In order to determine a proper order of the PCE, an order convergence study is presented here. Since the load capacity is the most important in an adhesively bonded joint, the probability density of the maximum bearable load of the polyurethane-to-steel butt joint predicted by different orders of PCE is shown in Figure 15. It shows that the results   Advances in Materials Science and Engineering generated by the 6 order and 8 order are quite different from the 10 order and 12 order. However, the 10 order result is relatively close to the 12 order, and 10 order of PCE is much computationally efficient (save 20% CPU time), thereby verifying 10 order being a reasonable approximation.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare there are no conflicts of interest regarding the publication of this research.