Abstract
The primary objectives of this study are twofold. Firstly, the original SPR method of stress recovery has been modified by incorporating the kriging interpolation technique to fit a polynomial to the derivatives recovered at the Gauss points. For this purpose, the -version of finite element analysis is performed to produce the stresses at the fixed Gauss points where the integrals of Legendre polynomials are used as a basis function. In contrast to the conventional least square method for stress recovery, the weight factor is determined by experimental and theoretical variograms for interpolation of stress data, unlike the conventional interpolation methods that use an equal weight factor. Secondly, an adaptive procedure for hierarchical -refinement in conjunction with a posteriori error based on the modified SPR (superconvergent patch recovery) method is proposed. Thirdly, a new error estimator based on the limit value approach is proposed by predicting the exact strain energy to verify the kriging-based SPR method. The validity of the proposed approach has been tested by analyzing two-dimensional plates with a rectangular cutout in the presence of stress singularity.
1. Introduction
The error assessment tools used in finite element analysis are well known and usually classified into two strategies: recovery-based error estimators and residual-type estimators [1, 2]. Stress recovery procedures can be classified as local (i.e., element level), patch-based, and global. To obtain a smooth stress field, averaging either projected or consistent finite element nodal stresses is an example of a patch-based scheme. The ideas of Zienkiewicz and Zhu [3–5] using the superconvergent patch recovery (SPR) are often preferred by researchers since they are robust and simple to use. Some references [6–11] contain extensive reviews of the different proposals for improving the SPR technique. They used the conventional LSM (least square method) to obtain recovered stresses from the -version of finite element solution at the sampling quadrature points. However, the residual-type error estimators have been proposed to evaluate errors for high-order hierarchical elements [12–14]. The residual error for high-order hierarchical elements is the difference between the displacement fields over the original and a refined mesh and is computationally more expensive than the (Zienkiewicz and Zhu) error estimate.
Recently, some researchers [6, 7] used a weighted superconvergent patch recovery technique in which the recovered stresses are calculated by using weighting parameters. In addition to these, the recover-based technique has been also extended to mesh based PUMs and the X-FEM (extended finite element method) [15]. Rodenas et al. [16] also explored the capabilities of a recovery technique based on an MLS (moving least squares) fitting, more flexible than SPR techniques as it directly provides continuous interpolated fields without relying on an FE mesh, to obtain estimates or the error in energy norm as an alternative to SPR. In the context of FEM model, the kriging interpolation technique has been employed as an alternative for estimating the derivative of the unknown variable at any point of interest [17]. This method uses a variogram to express the spatial correlation, and it minimizes the error of predicted values. It estimates the value at a location of interest as a weighted sum of data values at surrounding locations. The weight factors are assigned according to the variogram function that gives a decreasing weight with increasing distance between the given location and one of the surrounding locations. Dai [18] used the kriging interpolation in mesh-free methods in order to compare with the radial point interpolation method (RPIM) based on local supported radial basis function (RBF) and the Galerkin weak form. The literature on kriging interpolation for FEM, however, is very limited [14].
The -adaptive finite element analysis based on the error estimation consists of two stages: a posteriori error estimation and the automatic mesh refinement. The goal is to increase the -level nonuniformly so that the error is within the specified tolerance. The estimated errors in the finite element solution are of primary importance because of the basis for adaptive mesh refinement [19–21]. To minimize the computational cost, an effective and reliable technique of postprocessing is necessary for use in adaptive mesh refinement. It is known that the error estimate has not been directly extended to the -refinement [12, 13, 22, 23] because the high-order shape functions used to interpolate displacements within an element are also used to interpolate recovered stresses.
The objective of this study is to demonstrate the applicability of OK (ordinary kriging) interpolation to the -adaptive refinement of L-shaped domain problem employing the modified SPR method for stress recovery. To verify this method, the limit value approach is proposed to predict the exact strain energy for nonsmooth problems based on the application of the equation of a prior error indicator in the asymptotic range to three FEMs with three successively higher levels of polynomial approximation.
2. Ordinary Kriging Interpolation
The OK method is a geostatistical interpolation technique that requires both the distance and the degree of variation between known sampled data points when estimating values at unsampled locations, in other words, firstly, to calculate the distances between the predicted unknown point and the measured points nearby and, secondly, to derive the weight of each of these surrounding measured points by using the value of the variogram against those distances. The derived weights at unknown points result in optimal and unbiased estimates by minimizing the error variance. defined in (1) is often called a semivariogram or semivariance that can be defined as half the expected squared difference between paired random functions denoted by stresses in FEM, separated by the distance and direction vector called by lag or separation distance [24, 25] such aswhere is the number of pairs of values of which the separation distance is marked with . When the semivariance is plotted against the lag distance or separation distance between points, the plot is called semivariogram. For any given set of spatial data separated from , there will generally be at most one pair that are separated by a given distance . However, the observed data may be scattered near the separation distance . One must necessarily aggregate point pairs with similar distances where is called the allowable limit of separation distance as shown in Figure 1. Thus the separation distance can be allowed to use the similar distance denoted by . Sample data belonging to a certain interval are averaged to find the representative value at a given distance .
The semivariogram model is a function of three parameters, known as the nugget effect, sill, and range. Theoretically, at zero separation distance, the semivariogram value should be zero. However, at an infinitesimally small separation distance, the difference between measurements often does not tend to zero. This is called the nugget effect. Thus, the theoretical semivariogram model with a nugget effect can be fitted where the sill denoted by means the maximum semivariogram value that is the plateau of Figure 2. As the separation distance of two pairs increase, the semivariogram of those two pairs also increases. Eventually, the increase of the separation distance cannot cause the semivariogram increase. The separation distance when the semivariogram reaches plateau is called range . The example of schematic variogram graphs has been plotted in Figure 2 with respect to different variogram models. Several theoretical kriging modules are shown in Table 1.
As mentioned earlier, the OK estimates are linear weighted moving average of available observations or Gauss points in the FEM [23, 24]where and are the weights assigned to the available observations and neighbor data close to the unsampled location , respectively. The weight factors add up to unity to ensure that the estimate is unbiased.When a calculated stress value at any point is and the corresponding true value is , an error variance based on OK technique is as follows:Equation (4) can be written as belowSubstituting (2) in (5) givesThen, (6) can be written as belowEquation (7) can be rewritten asIf variance and covariance in the above equations are marked as belowthen (8) is written as follows:The error variance associated with the OK estimate is called the minimum variance unbiased estimator or best linear unbiased estimator, since the constraint condition defined in (3) should be applied to minimize the variance of estimate errors. Based on the method of Lagrange multipliers, the mathematical form considering (3) and (10) can be expressed aswhere is a Lagrange objective function and is a Lagrange multiplier. In addition, number 2 in the fourth term in (11) is used to derive final equations with simple form. Minimizing the objective function can be carried out by finding the partial derivatives with respect to and such thatEquation (12) can be rearranged to the following form:From (13), weight factors for any unknown stress are calculated by following formFinally, the unknown stresses for points can be obtained from
3. A -Adaptive Refinement Using Modified SPR Technique
Adaptive procedures are to implement iteration analysis based on distinctly different levels of space ()- or function ()-refinements in specified local region to achieve solutions having a certain degree of accuracy in an optimal fashion. Particularly the -adaptive refinement makes it easy to use the initial meshes kept unchanged with selective increase in the polynomial order of shape function that is called “selective or nonuniform -refinement” in FEM. The continuity between meshes with different polynomial order is achieved by assigning zero to the higher-order derivatives associated with edges in common with the lower derivatives. Thus, the higher-order of approximation is degraded to the lower along the interelement boundary. For this purpose, two important algorithms should be established such as an automatic -adaptive mesh refinement scheme and a posteriori error estimator for -refinement strategy.
In this work, the higher-order approximation based on Lobatto shape functions [26] which is often called integrals of Legendre polynomials [14, 23] with hierarchical properties is adopted to obtain displacements as a result of FEM. The SPR technique proposed by Zienkiewicz and Zhu [3, 5] has been adopted after a suitable modification to be compatible with the adaptive p-refinement procedure since the number of sampling Gauss points in each element is increased as the -level becomes higher. The increment of the quadrature point has an effect on the stress norm. According to , represents the local error for a particular element, measured in energy norm. is the recovered stress resultant field or estimated exact stress field over the patch of elements (normally consisting of 4 elements) surrounding the patch node. The estimated exact stress field denoted by is determined by using the ordinary kriging interpolation technique. A posteriori error estimate in a particular energy norm is computed by summing its elemental contribution aswhere is a column vector including the true stresses; including the stresses interpolated by Lobatto shape functions with the displacements obtained by FEM; [] is a constitutive matrix; is a mesh domain. Here a smoothed continuous stress concept is applied for the true stresses that are obtained by the aforementioned OK technique. For the proposed technique, the weight factors depending on the distance between a sampling point and the point corresponding to an unknown value are considered. On the other hand, no weight factor is used in the original SPR technique based on LSM. The energy norm of the stress field itself may also be expressed in terms of stresses as follows:Thus, the relative percentage error can be defined as
In case of the modified SPR method using ordinary kriging, however, the estimated exact stress field denoted by should be calculated at each iteration round of -refinement. In other words, cannot be fixed in the whole process of -adaptive refinement the same as -adaptive refinement since this process is based on a posteriori error estimate. Thus, a new error estimator is proposed on the basis of a prior error estimator to verify the modified SPR method. The limit value approach is proposed to predict the exact strain energy for nonsmooth problems based on the application of the equation of a prior error indicator in the asymptotic range to three FEMs with three successively higher levels of polynomial approximation. In the presence of singularities, the asymptotic convergence behavior of the -version of the FEM permits a close estimate of the exact strain energy by extrapolation that is called the limit value, and hence we can predict the error in the energy norm on the finite element mesh employed. For a two-dimensional problem, under the assumption that the error in the energy norm has entered the asymptotic range where and are the strain energy, the rate of convergence for the p-version of FEM can be derived by the inverse theorem [23, 27] aswhere and are the exact strain energy estimated by the limit value and the approximate strain energy by FEM, is the strength of singularity, and and are the degrees of freedom for the polynomial order and a constant which depends on the mesh, respectively. There are three unknowns , , and in (19). By performing three successive extension processes, , , and , which are in the asymptotic range, we have three equations for computing the unknowns. Cancelling and in (19), the following extrapolation equation can be derived aswhere , , and are the strain energies when the polynomial orders are , , and and represents the limit value in terms of estimated exact strain energy where , , and are the number of degrees of freedom for each analysis. Computational experiences show this estimated limit value to be reliable and accurate for two-dimensional elastostatic problems, especially in the presence of singularity [14, 23] which give no exact solution. Thus, the percentage relative error expressed by energy norm using limit value is defined by (21). It is noted that (21) is not used for local indicator in practice but is used only to validate whether the modified SPR technique is reliable and accurate. If the local error estimate in (16) is large, the polynomial order should be increased to satisfy an acceptable level of accuracywhere is the exact global strain energy and is also the global strain energy calculated in the current -adaptive mesh consisting of nonuniform -distribution at a certain iteration number.
4. Numerical Analysis
The numerical example is shown in Figure 3 that specifies the geometric definition and analysis conditions that are given by cm, cm, N/cm2, , and N/cm2. Due to the symmetry, a quarter of plate is modeled by three elements with fourth-degree polynomials for shape functions of -version finite element analysis. Figure 4 shows a 3D stem plot for distribution of von-Mises stresses at Gauss points. The stress values obtained from finite element analysis are considered as raw data for stress smoothing.
(a) L-shaped domain
(b) Finite element mesh by -FEM
In Figure 5, the experimental semivariogram has been plotted with respect to the separation distance . The allowable limit of separation distance is assumed by 25% of . Thus, the separation distance can be allowed to use the similar distance denoted by . Based on experimental semivariogram, three different theoretical semivariogram modules as shown in Table 1 have been tested. In this study, the Gauss model with % of has been adopted to find the weight factor for the OK process explained in (2).
Before the further analysis of -adaptive refinement, the performance between LSM and OK method is compared with each other in Figures 6 and 7. As described earlier, the FEM raw data represents the computed von-Mises stress at Gauss points obtained by the -version finite element model in Figure 3. Due to the discontinuity of stresses along the element interboundary, it is seen that the stress distribution is not smooth as shown in Figure 6(a). Thus, the stress recovery techniques are applied for stress smoothing. One is LSM based on equal weighted interpolation and the other is OK method by weighted interpolation using variogram model. It is noted that the corner singularity denoted by von-Mises stress is well expressed by OK interpolation comparing with LSM that are shown in Figures 6(b), 6(c), 7(b), and 7(c).
(a) Raw data
(b) OK interpolation
(c) LSM
(a) Raw data
(b) OK interpolation
(c) LSM
To illustrate the applicability of OK interpolation to the -adaptive mesh refinement, two -version models are considered by 3-element and 12-element model. The initial -level of both models begins with one. For nonuniform -distribution, the continuity between elements with different polynomial orders is achieved by assigning zero higher-order derivatives associated with the edge in common with the lower-order derivatives. The iteration step for -adaptivity is proceeded to final adaptive mesh based on a posteriori error estimation. In this study, the modified SPR technique is proposed to estimate the smoothed stress field by projection that is considered as an exact solution to calculate a posteriori error. The final adaptive mesh is automatically determined by the developed computer program for the purpose that are shown in Figures 8 and 9.
(a) OK interpolation
(b) LSM
(a) OK interpolation
(b) LSM
The relative percentage errors have been illustrated in Tables 2 and 3 according to iteration numbers by using the modified SPR method as well as the limit value approach. The first error estimator by the modified SPR method is the objective of this study considering a posteriori error estimator as shown in Table 2. However, the relative percentage errors are not converged gradually, especially when higher-order polynomials are used in the -adaptive mesh regardless of least square method or OK interpolation since the estimated exact stress field can be reproduced at each iteration round. To verify the proposed error estimator, the limit value approach based on a prior error estimate is used to evaluate the relative percentage error in Table 3. Since the exact strain energy is fixed as a certain value, the relative percentage errors are gradually decreased. As shown in Tables 2 and 3, the required number of iterations to determine the final adaptive mesh is 10 by the LSM based -adaptive model and 8 by OK based p-adaptive model. From this result, it is observed that the OK interpolation technique is more suitable for adaptivity procedures than the LSM that has been commonly used in the FE analysis.
In the case of LSM with NDF = 363, the relative percentage error shows 5.48%, but OK yields 4.17% when NDF = 353 as shown in Table 2. However, the relative percentage error based on the limit value is 5.78% for LSM and 5.61% for OK, respectively, as shown in Table 3. It is observed that the exact strain energy can be found by from Figure 10 that has been used to estimate the relative percentage error at each iteration round. In other words, in (21), is the exact global strain energy, and is the global strain energy calculated in the current -adaptive mesh consisting of nonuniform -distribution at a certain iteration round. Thus, (21) cannot be used for local indicator in practice but is only used to check whether the modified SPR method is reliable and accurate. The convergence characteristics of proposed -adaptive models are plotted in Figures 11-12. It is noted that stable convergence pattern by the limit value approach is shown without any numerical oscillations, unlike the modified SPR method.
(a) 3-element model
(b) 12-element model
(a) 3-element model
(b) 12-element model
5. Conclusions
The new -adaptive finite element model with the OK interpolation technique is proposed in this study. This approach shows better performance to determine the final adaptive mesh than the existing SPR method by using LSM to estimated exact stress field by projection. This proposed model is very suitable for stress singularity problems since the higher weight factor is used to interpolate the calculated stresses at the Gauss points near stress singular points by applying theoretical variogram model and OK interpolation. In addition, the proposed new error estimator based on limit value using exact strain energy can be used to a prior error estimate.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of the paper.
Acknowledgments
This work was supported by the 2015 Yeungnam University Research Grant.