Optimized Polynomial Virtual Fields Method for Constitutive Parameters Identification of Orthotropic Bimaterials

Heterogeneous materials are widely applied in many fields. Owing to the spatial variation of its constitutive parameters, the mechanical characterization of heterogeneous materials is very important. +e virtual fields method has been used to identify the constitutive parameters of materials. However, there is a limitation: constitutive parameters of one material have to be a priori; then, constitutive parameters of the other one can be identified. Aiming at this limitation, this article presents a method to identify the constitutive parameters of heterogeneous orthotropic bimaterials under the condition that constitutive parameters of both materials are all unknown from a single test. A constitutive parameter identification method of orthotropic bimaterials based on optimized virtual field and digital image correlation is proposed. +e feasibility of this method is verified by simulating the deformation fields of a two-layer material under three-point bending load. +e results of numerical experiments with FEM simulations show that the weighted relative error of the constitutive parameter is less than 1%. +e results suggest that the variation coefficient-to-noise ratio can perform a priori evaluation of a confidence interval on the identified stiffness components. +e results of numerical experiments with DIC simulations show that the weighted relative error is 1.44%, which is due to the noise in the strain data calculated by DIC method.


Introduction
Heterogeneous materials, such as functionally gradient materials, wood [1,2], polymeric foams [3], biological materials [4], and composites [5], have been widely applied in engineering. Constitutive parameters identification of heterogeneous materials is not only the important subject in experimental mechanics but also has attracted extensive notice in the fields of solid mechanics, structure health monitoring, and medical diagnosis. e traditional method of identifying constitutive parameters can directly determine all constitutive parameters of constitutive models through standard experiments (e.g., uniaxial tensile test). However, the investing of experimental involved in such standard test methods may increase when using anisotropic models and inhomogeneous materials, as well as bimaterial [6].
With the development of the photomechanics, full-field measurement methods (e.g., digital image correlation and moiré interferometry) combined with inverse methods have been developed so that the multiple constitutive parameters can be identified by a heterogeneous strain field [7]. e inverse methods belong to the updating methods (e.g., finite element updating method) and nonupdating methods (e.g., virtual fields method). e finite element updating method (FEUM) extract constitutive parameters through minimizing the cost function between the response of the finite element model and the real behavior iteratively; its simulation must be close to experimental conditions and its initial guess choice will affect the identification result [8]. e virtual fields method (VFM), which was put forward by Perrion and Grédiac [9,10], is a typical nonupdating method used in a linear material model. Compared with the updating methods, the advantage of VFM is that it requires less specimen boundary conditions and geometric configurations [11,12].
A volume of research has focused on full-field measurement and VFM to identify the constitutive parameters of materials.
e VFM within peridynamic framework was proposed to identify material properties in linear elasticity [13]. e special virtual fields method was adopted to resolve the difficulty in identifying constitutive properties of incompressible and nonhomogeneous solids [14]. e VFM was combined with an elaborately designed test configuration in order to realize the identification of the anisotropic yield constitutive parameters [15]. Several inverse methods based on the VFM have been developed, including the Fourier series-based VFM [16,17], the sensitivity-based VFM [18], the eigenfunction VFM [19,20], and several optimizations to improve the accuracy of identification results [21,22]. e heterogeneity of strain field caused by multiple constitutive parameters of the complex constitutive model has been solved by the traditional VFM. For heterogeneous materials, the heterogeneity of strain field caused by the spatial variation of constitutive parameters requires further improvement of the virtual field method. According to the available literature studies, several methods and applications have been studied to identify constitutive parameters identification of heterogeneous material using VFM [6,23,24]. For instance, a bimaterial constitutive parameter identification method based on the combination of VFM and Moiré interference was proposed [25], and the VFM was also extended to identify the parameters [26]. In addition, several methods have been researched to lower the effects of noise in strain fields [26][27][28][29]. e optimized fields method, such as polynomial optimized virtual field [30] and piecewise optimized virtual field [31], was developed to improve convenience and accuracy. Although the above researches make identification of constitutive parameters more rapid and convenient, there is still a limitation for bimaterials: the constitutive parameters of one material must be a priori in order to determine the constitutive parameters of the other material. e purpose of this work is to propose an optimized virtual field that can not only extract the constitutive parameters of the heterogeneous orthotropic bimaterials without knowing any material constitutive parameters but also require only one single test.
In this study, to extract constitutive parameters of the bimaterials under the condition that constitutive parameters of the both materials are unknown, special optimized polynomial virtual fields combined with threepoint bending test was proposed. In Section 2, the basic principle of optimized polynomial VFM and its improvement applied in bimaterials without knowing any material constitutive parameters are introduced. In Section 3, the FEM simulated experiments and the simulated DIC experiments are conducted to verify the feasibility of the abovementioned method. Finally, Section 4 is the conclusions.

Basic Principle of Optimized Polynomial Virtual Fields
Method.
e integral form of the mechanical equilibrium equation of the VFM based on virtual work for a continuous solid can be expressed as [31] where σ is the Cauchy stress tensor, u * is the virtual displacement vector, ε * denotes the strain tensor corresponding to u * , T represents the external force associated with the stress tensor σ by the boundary, S is a vector which is the volume force applied over the volume V of the specimen, and a denotes the distribution of acceleration. Such a distribution will cause an additional volumetric force distribution equal to −ρa with D'Alembert's principle, the virtual displacement field being kinematically admissible (KA).
e constraints of KA virtual fields u * (1) , u * (2) , u * (3) , and u * (4) are as follows: e exact values and white noise constitute the measured strain component. Assume that the noise components are uncorrelated to each other, and presume that the noise between points is also uncorrelated. erefore, the principle of virtual work is as equation (5) (see Appendix A2 for the detailed equation derivations).
where N 1 , N 2 , and N 6 denote the processes of scalar zeromean stationary Gaussian generalized by the corresponding strain components ε 1 , ε 2 , and ε 6 , respectively, and c denotes the measured amplitude of the random variable strain.
If the noise is ignored, approximate parameters expressed as Q app , these components are defined by equation (6). So, directly identify Q 11 , Q 22 , Q 12 , and Q 66 using these four special virtual displacement fields as equation (7) (see Appendix A3 for the detailed equation derivations).

Advances in Materials Science and Engineering
Similar results are obtained for Q 22 , Q 12 , and Q 66 . Denoting V(Q),the vector containing the variances of Q 11 , Q 22 , Q 12 , and Q 66 , one can write as equation (8) (see Appendix A4 for the detailed equation derivations).
where A ij and B ij are the coefficients of the monomials, m and n are the polynomial degrees that define that the maximum number of monomials has used, and L and w are the typical dimensions of the x 1 and x 2 directions, respectively.
Using the Lagrange multiplier approach, the Lagrangian function L (i) can be constructed for each constitutive parameter sought, and its constraint is equation (13), and the objective function is (η (i) ) 2 . erefore, the expression of the Lagrange function is where λ (i) is the vector containing Lagrange multipliers. e KA condition can produce several linear equations, and the number of this equations is depending on the number of supports. e special virtual field condition also leads to several linear equations, the number of which depends on the type of constitutive model of specimen material. e two types of conditions can produce the following linear system as equation (14) (see Appendix A5 for the detailed equation derivations).
Attention that (η (i) ) 2 coefficients also correspond on the unknown Q ij parameters. us, the optimization problem is solved for a given set of Q ij parameters, and the process is iterative: the first guess of the Q ij parameters is used to find the first group of four special virtual fields. en, these four special virtual fields are used to find a new set of Q ij parameters, and repeat abovementioned steps until the optimal Q ij parameter is found. In practice, no matter what the parameter Q ij is initially selected, this iterative process will converge quickly. Usually, only two loops are sufficient to converge.

Optimized Polynomial Virtual Fields Applied to the Bimaterials without Knowing Any Material Constitutive
Parameters. As shown in Figure 1, for the plane-stress specimen with two different materials Part A and Part B, and both of them are set with the elastic orthotropic material. But the constitutive parameters of A and B are unknown. As mentioned above, the sum of exact values and white noise is the measured strain. Assume that the noise components are uncorrelated to each other, and presume that the noise between points is also uncorrelated. erefore, the governing equation of the virtual fields method for the total specimen is as equation (15) (see Appendix B1 for the detailed equation derivations).
e constraints of the eight special virtual fields denoted that u * (1) , u * (2) , u * (3) , u * (4) , u * (5) , u * (6) , u * (7) , and u * (8) must satisfy the following equation: Advances in Materials Science and Engineering S a , and Q 66b using the abovementioned four special virtual displacement fields. If the noise is ignored, approximate parameters expressed as Q app are identified. e components of which are defined as equation (17) (see Appendix B2 for the detailed equation derivations). 6 Advances in Materials Science and Engineering To minimize the effect of noise, taking Q 11a as an example, the variance of Q 11a is as follows: Advances in Materials Science and Engineering As described above, the virtual fields are expanded by polynomial basis functions based on equation (11). e unknown coefficients A ij and B ij are included in vector Y as follows: e integrals of equation (15) are calculated using vectors B11, B22, B12, and B66 through discrete sum approximation so that e others that the definition from the virtual fields (vectors H11, H22, H12, and H66) are related to the terms of equation (18).
e first KA virtual field conditions of three-point bending are u * 2 � 0, when x1 � 0, and x2 � 0, which implies e second KA virtual field conditions of three-point bending re u * 1 � 0, when x1 � L, which means that e optimized virtual fields are totally defined. e following equations can identify the constitutive parameters: 8 Advances in Materials Science and Engineering

Numerical Verification and Discussion
To evaluate the accuracy of the abovementioned optimized polynomial VFM, numerical experiments with FEM simulations and DIC simulations are employed in this section.

Numerical Verification of Optimized Polynomial Virtual
Fields with FEM Simulations. In FEM simulation, in order to facilitate the simulation of bimaterials configuration, a twodimensional rectangular beam specimen was adopted, and applied a three-point bending load to get inhomogeneous in-plane deformation fields. Based on the principle of the symmetry of three-point bending, only the left of the specimen was modelled, as shown in Figure 1. e length, width, and thickness of the bimaterial beam specimen were set as L � 30 mm, w � 20 mm, and t � 2.3 mm, respectively, and the width of the two material regions was the same: w a � w b � 10 mm. e force was set as F � 2544N. e four-node plane-stress element with thickness (Solid-Quad 4 nodes 182/Plane strs w/thk) was employed. e displacement constraint u y � 0 is added in the lower left corner of the model, and the line constraint u x � 0 is added on the right side of the model. e beam specimen was set as orthotropic materials. e engineering elastic constants in the FEM simulation are shown in Table 1, and the converted constitutive parameters in the elastic matrix are shown in Table 2.
e displacement fields and the strain fields of the threepoint bending test of the orthotropic bimaterial obtained through FEM simulations using the commercial software ANSYS are shown in Figures 2 and 3. Table 3 shows the identification results of the heterogeneous orthotropic bimaterials by using the optimized polynomial VFM, as shown in Section 2.2. To study the identification accuracy of the optimized polynomial VFM, the relative error of every constitutive parameter component and the weighted relative error W re were calculated. Owing to there were 8 unknown constitutive parameters, the weighted relative error employed here can be expressed as follows: It can be seen from Table 3 that the weighted relative error is lower than 1%. e relative error of Q 66 of the two materials is the smallest. Although the relative errors of Q 11 and Q 22 of the two materials are slightly larger than Q 66 , both of them are less than 2%. e relative error of Q 12a and Q 12b is significantly higher than others, because σ 2 must be included in the virtual fields to identify Q 12 , while the stress data density of σ 2 is low in the three-point bending test. e stress data density can be expressed by the following equation: Advances in Materials Science and Engineering According to equation (7), the standard deviations of Q ij are equal to η (i) c. If c is the standard deviation of the strain noise, the coefficients of variations (CV(Q ij )) of each stiffness component are given by e ratio of coefficient of variation to the standard deviation of the strain noise can show the sensitivity to the noise. It indicates the rate at which the coefficient of variation increases with the increase of noise, which can be called the variation coefficient-to-noise ratio, and can be expressed as η ij /Q ij . e lower variation coefficient-to-noise ratio indicates that the corresponding identification constitutive parameters are less susceptible to noise. e variation coefficient-to-noise ratios η ij /Q ij of eight identification results are shown in Table 3. Significant differences on the variation coefficient-to-noise ratio and the coefficient-tonoise ratio of Q 66 is the smallest, the next is Q 11 , and then is the Q 22 , and the biggest is Q 12 . e quantitative results show that if Q 11 , Q 22 , and Q 66 are stable, Q 12 is less robustly identified. is result shows that this configuration is unsuitable for the transverse stiffness.

Comparison with Results of Piecewise Virtual Fields and Polynomial Virtual Fields.
is section will implement the same numerical experiment as Section 3.1 to compare the identified constitutive parameters of the piecewise virtual   10 Advances in Materials Science and Engineering field and polynomial virtual field. At this situation, the bilinear shape function is used to interpolate the four-node element to represent the virtual field instead of the polynomial virtual field. e piecewise function expansion of the virtual displacement field is similar to the interpolation of the real displacement in the FEM. Derived from the relationship between strain and displacement,

12
Advances in Materials Science and Engineering For the quadrilateral element, the virtual displacement can be expressed by the following formula: where u * is the vector containing the virtual displacement of any point (ξ 1 , ξ 2 ) in the element, and u * (a) is the vector containing the virtual displacement of the nodes in the element. e expression of N (a) is defined as where I is the identity matrix. e shape function expression of N (a) can be defined as follows: So, the virtual displacement can be expressed by the following formula: erefore, the unknown coefficient vector Y (equation (17)) in the polynomial virtual fields can be denoted on each node: where n is the total number of nodes. Use the piecewise virtual field instead of the polynomial virtual field to perform the same experiment (Section 3.1). Table 4 shows the values of relative and weighted relative error of A and B. It can be seen from Table 4 that the identification results of piecewise virtual fields are similar to that of polynomial virtual fields. e weighted relative error is 0.97%, which is slightly larger than that of polynomial virtual fields. e relative error of Q 66 of the two materials is the smallest. Although the relative errors of Q 11 and Q 22 of the two materials are slightly larger than Q 66 , both of them are less than 2%. Similar to the results of polynomial virtual fields, the relative error of Q 12a and Q 12b is significantly higher than others, because σ 2 must be included in the virtual fields to identify Q 12 , while the data density of σ 2 is low in the three-point bending test. e variation coefficient-to-noise ratios of eight identification results are also shown in Table 4. Significant differences can be seen on the variation coefficient-to-noise ratios and the quantitative results show that if Q 11 , Q 22 , and Q 66 are stable, Q 12 is less robustly identified. Compared with the polynomial virtual field, the variation coefficient-tonoise ratio of the eight identification parameters calculated by the optimized piecewise virtual field is larger.

Influence of Noise on Identification Results.
e noise is inevitable in the experiments. To investigate the influence of noise on identification results, zero-mean Gaussian white noise with different amplitudes is added in the strain data of the orthotropic bimaterial specimen, and the abovementioned experiments are repeated to obtain the coefficients of variation corresponding to the eight identification results with the increasing noise amplitudes. In order to assess the influence of noise on the identified parameters, it is necessary to discuss the standard deviation of strain noise. It can be seen from equation (27) that the coefficients of variation (CV(Q ij ) a ) are proportional to the uncertainty (c) of strain measurement. Since the coefficient of variation can be directly compared between different constitutive components, this article chose to discuss the coefficient of variation rather than standard deviation. us, coefficients of variation were used to discuss the sensitivity to noise. e coefficients of variation of the identified stiffness components of the two materials were calculated for a range of strain noise standard deviation values. e strain noise standard deviation c varies from 5 × 10 −5 to 1 × 10 −3 by steps of 5 × 10 −5 . For each value of c, the identification is repeated 30 times using the optimized polynomial virtual fields. e coefficients of variation (CV(Q ij )) were calculated by equation (34). Figure 4 indicates the graph plotting, the coefficients of variations for the eight stiffness components of the two materials, as a function of c. e points are fitted by linear regression to calculate the slope, and the slope was compared to the theoretical value of the variation coefficient-to-noise ratio η ij /Q ij .

Advances in Materials Science and Engineering
For both materials A and B, Figure 4 shows that the smallest coefficient of variation value is CV(Q 66 ). e results of Figure 4 are consistent with the variation coefficient-to-noise ratio studies in Section 3.1. Q 66 was the most stable of the identified stiffness components. CV(Q 22 ) was the second and CV(Q 11 ) was the third, and finally CV(Q 12 ), which is also consistent with the results in Section 3.1. e slope of the fitted curve is the fitted variation coefficient-to-noise ratio. e comparisons of the fitted and theoretical value of the variation coefficient-to-noise ratio are listed in Table 5. Figure 4 and Table 5 not only show the linear relationship but also show that the theoretical η ij /Q ij values are consistent with the fitted ones. It suggests that the procedure can perform a priori evaluation of a confidence interval on the identified stiffness components. In addition, due to the hypothesis of equation (16), the straight line is unsuitable for the highest values of noise. Equation (16) shows the noise should be smaller than the signal.

Influence of the Polynomial Degrees.
According to equation (10), when the basic functions to expand the virtual field is selected, the only parameters, the degrees of the x 1 and x 2 monomials, denoted as m and n, need to be selected. m and nare the polynomial degrees that represent the maximum number of monomials. e total number of polynomial degrees is 2(m + 1)(n + 1); for the case of optimized polynomial virtual fields applied to the bimaterial without knowing any material constitutive parameters, the condition of equation (35) needs to be satisfied.
2(m + 1)(n + 1) > n + 10. (36) For the optimized polynomial virtual fields, on the one hand, m and n cannot be selected too small, which will lead to too few polynomial degrees and thus cannot satisfy equation (35). On the other hand, m and n should not be selected too large, which will result in ill-conditioned matrix. Figure 5 shows the evolution of the variation coefficient-tonoise ratio η ij /Q ij as a function of m and n. It can be seen

16
Advances in Materials Science and Engineering from Figure 5 that there are slight variations in the η ij /Q ij with higher gradients for the values related to Q 12a and Q 12a . Due to the three-point bending test configuration, these parameters are relatively less good identifiability. It can also show that the dependance on n is higher than that on m. is is because all the constraints affect the x 2 monomials.

Numerical Verification of Optimized Polynomial Virtual
Fields with DIC Simulations. In practice, to provide a more reliable strain input for the virtual fields method, digital image correlation (DIC) was applied to obtain deformation fields. DIC is an effective optical technique for noncontact full-field deformation measurement for various materials and structures. e DIC computes the displacement of each image point by comparing the gray intensity of images of the test specimen surface in different loading states. e corresponding strain fields are calculated using a numerical differentiation method. In addition to the systematic error of DIC algorithm, aberration, distortion, and out-of-plane displacement will affect the accuracy of displacement measurement.
For the purpose of evaluating the simulated results of constitutive parameters for bimaterials and the influence of  strain input calculated by DIC on the identification of constitutive parameters, the simulated DIC results were used as the input of the VFM to exclude the influence of lens distortion and out-of-plane displacements. e displacement fields of the bimaterial beam specimen obtained through abovementioned FEM simulations were exported and converted to reference and deformed speckle patterns by Gaussian speckle [32], as shown in Figure 6. e displacement fields and strain fields obtained from DIC are shown in Figures 7 and 8, respectively. Table 6 shows the identification results and the corresponding variation coefficient-to-noise ratio and relative error and weighted relative error for both materials. Similar to the results in Table 3, it shows significant differences on the variation coefficient-to-noise ratio and the coefficient-to-noise ratio of Q 66 , which is the smallest, the next is Q 11 , and then is the Q 22 , and the biggest is Q 12 . e lowest value is Q 66 . is result is consistent with studies in Section 3.3. Q 66 was the most stable of the four identified stiffness components in Section 3.3. Q 11 is the second, Q 22 is the third, and Q 12 is the last. is result is also consistent with the studies in Section 3.3. As shown in Table 6, the relative errors of Q 12a and Q 12b are more than 10%, while others are not more than 4%, and the weighted relative error is 1.14%, which is slightly larger than the results in Table 3. is indicates that the strain data calculated by DIC contains noise, resulting in an increase in the error of the identification results. is comparison result reveals the feasibility of the optimized polynomial virtual fields combined with DIC for extracting the constitutive parameters of the unknown bimaterial. It should be pointed out that in the actual DIC calculation, the displacement measurement errors caused by aberration and out-of-plane displacements will lead to additional strain noise, so the relative errors and weighted relative error of the identification results will be larger than the results in Table 6.

Conclusions
In this article, optimized polynomial virtual fields are proposed to extract the constitutive parameters of the heterogeneous orthotropic bimaterials without knowing any material constitutive parameters. Numerical experiments with FEM simulations are employed to investigate the accuracy of the optimized polynomial virtual fields method. e main conclusions are as follows:  (1) e optimized polynomial virtual fields were proposed to identify the constitutive parameters of heterogeneous orthotropic bimaterials under the condition that constitutive parameters of both materials are all unknown from a single test. e results show that the weighted relative error of the constitutive parameter is less than 1%. e stress data density was used to explain the difference in the relative error of the constitutive parameter. For heterogeneous orthotropic bimaterials without knowing any material constitutive parameters, this method is suitable for extracting the constitutive parameters.
(2) To investigate the effect of the noise, a series of numerical experiments with different strain noise amplitudes were used in FEM simulations. e variation coefficient-to-noise ratio and coefficients of variation are applied to quantify the influence of noise on the identification results, and the results suggest that the variation coefficientto-noise ratio can perform a priori evaluation of a confidence interval on the identified stiffness components.
(3) For comparison, piecewise virtual fields were used to extract the constitutive parameters of the heterogeneous orthotropic bimaterials. e results denoted that the optimized polynomial virtual fields have a higher reliability than the optimized piecewise virtual fields. (4) To study the accuracy of the optimized polynomial VFM combined with DIC method, the numerical experiments with DIC simulations are employed. e comparison result shows the feasibility of the optimized polynomial virtual fields combined with DIC for extracting the constitutive parameters of the unknown bimaterials. e result also indicates that the strain data calculated by DIC contains noise, resulting in an increase in the error of the identification results, and the weighted relative error is 1.14%.

Nomenclature σ:
Cauchy stress tensor u * : Virtual displacement vector ε * : Strain tensor T: External force S: Area of the specimen b: Vector of volume force V: Specimen volume a: Distribution of acceleration KA: Virtual displacement field being kinematically admissible Q 11 , Q 22 , Q 12 , Q 66 : In-plane constitutive parameters u * (1) , u * (2) , u * (3) , u * (4) : Four independent KA virtual fields ε * (1) , ε * (2) , ε * (3) , ε * (4) : Four independent KA virtual fields N 1 , N 2 , N 6 : Processes of scalar zeromean stationary Gaussian ε 1 , ε 2 , ε 6 : Strain components c: Measured amplitude of the random variable strain Q app : Approximate parameter components ‖ε i ‖(i � 1, 2, 6): Norm of strain component where σ is the Cauchy stress tensor, u * is the virtual displacement vector, ε * denotes the strain tensor corresponding to u * , T represents the external force associated with the stress tensor σ by the boundary, S is a vector which is the volume force applied over the volume V of the specimen, and a denotes the distribution of acceleration. Such a distribution will cause an additional volumetric force distribution equal to −ρa with D'Alembert's principle, the virtual displacement field being kinematically admissible (KA).
In the case of a quasistatic load, acceleration acan be omitted. e volume force bcan usually also be omitted. Under this circumstance, equation (A.1) can be written as follows: Equation (A.2) denotes the plane-stress problem. For an orthotropic material in a plane-stress state, the principle of virtual work becomes where Q 11 , Q 22 , Q 12 , and Q 66 are the in-plane stiffness components and equation (A.3) is linear. It is suitable for every KA virtual field. is equation uses four independent KA virtual fields ε * (1) , ε * (2) , ε * (3) , and ε * (4) instead of u * (1) , u * (2) , u * (3) , and u * (4) , respectively. So, the relationship between A, Q, and B is as follows: A2.
where N 1 , N 2 , and N 6 denote the processes of scalar zeromean stationary Gaussian generalized by the corresponding strain components ε 1 , ε 2 , and ε 6 , respectively, and c denotes the measured amplitude of the random variable strain. Assume that the noise components are uncorrelated to each other, and presume that the noise between points is also uncorrelated. Equation (A.5) can be rewritten as follows: A3. e Principle of Virtual Work with Ignored Noise. Directly identify Q 11 , Q 22 , Q 12 , and Q 66 using these four special virtual displacement fields ε * (1) , ε * (2) , ε * (3) , and ε * (4) . If the noise is ignored, approximate parameters expressed as Q app , these components are defined by where c can be negligible compared to the ‖ε 1 ‖, ‖ε 2 ‖ and ‖ε 6 ‖, respectively, and the ‖ε i ‖(i � 1, 2, 6) is the norm of strain component L 2 . us, Equations (A.7)-(A.10) mentioned above can be reexpressed by the actual values of the Q ij instead of the approximate counterparts. us, (A.14) Using the rectangle method to discretize the abovementioned integrals,

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 that there are no conflicts of interest regarding the publication of this paper.