# p-Adaptive Refinement Based on Stress Recovery Technique Considering Ordinary Kriging Interpolation in L-Shaped Domain

Authors: Kwang S. Woo; Jae S. Ahn
Journal: Mathematical Problems in Engineering
 (2017)

---

## 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, thep-version of finite element analysis is performed to produce the stresses at the fixed 10×10 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 p-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.


---

## Body

## 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 h-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 Z/Z (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].Thep-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 p-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 Z/Z error estimate has not been directly extended to the p-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 thep-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.γ(h) defined in (1) is often called a semivariogram or semivariance that can be defined as half the expected squared difference between paired random functions Sxi denoted by stresses in FEM, separated by the distance and direction vector called by lag or separation distance h [24, 25] such as(1)γh=12n∑i=1nSxi-Sxi+h2,where n is the number of pairs of values of which the separation distance is marked with h. 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 S(xi+h) separated from S(xi), there will generally be at most one pair that are separated by a given distance h. However, the observed data may be scattered near the separation distance h.One must necessarily aggregate point pairs Sxi,Sxi+h with similar distances h±Δh where Δh is called the allowable limit of separation distance as shown in Figure 1. Thus the separation distance h can be allowed to use the similar distance denoted by h±Δh. Sample data belonging to a certain interval h±Δh are averaged to find the representative value at a given distance h.Figure 1
Definition of the allowable limit of separation 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 byC0 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 a. 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.Table 1
Different models of variograms.
Model
Theoretical semivariogramγ(h)
Linear
C
o
h
a for h≤a Co for h>a
Spherical
C
o
1.5
h
a
-
0.5
h
a
3for h≤a Co for h>a
Exponential
C
o
1
-
exp
⁡
-
3
h
a
2
Gaussian
C
o
1
-
exp
⁡
-
3
h
a
2
Polynomial
k
o
+
k
1
h
+
k
2
h
2
+
⋯
+
k
n
h
nFigure 2
Comparative plot of theoretical semivariogram models.As mentioned earlier, the OK estimates are linear weighted moving average of available observations or Gauss points in the FEM [23, 24](2)S∗xo=∑i=1nλiSxi,where λi and S(xi) are the weights assigned to the available observations and neighbor data close to the unsampled location x0, respectively. The weight factors add up to unity to ensure that the estimate is unbiased.(3)∑i=1nλi=1.When a calculated stress value at any point is s0 and the corresponding true value is s0∗, an error variance based on OK technique is as follows:(4)σOK2=Es0-s0∗2.Equation (4) can be written as below(5)σOK2=var⁡s0+var⁡s0∗-2cov⁡s0,s0∗.Substituting (2) in (5) gives(6)σOK2=var⁡s0+var⁡∑i=1nλisi-2cov⁡s0,∑i=1nλisi.Then, (6) can be written as below(7)σOK2=var⁡s0+cov⁡∑i=1nλisi,∑j=1nλjsj-2cov⁡s0,∑i=1nλisi.Equation (7) can be rewritten as(8)σOK2=var⁡s0+∑i=1nλi∑j=1nλjcov⁡si,sj-2∑i=1nλicov⁡s0,si.If variance and covariance in the above equations are marked as below(9)σ2=var⁡s0;σij2=cov⁡si,sjthen (8) is written as follows:(10)σOK2=σ2+∑i=1n∑j=1nλiλjσij2-2∑i=1nλiσ0i2.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 as(11)Lλ1,λ2,…,λn,μ=σ2+∑i=1n∑j=1nλiλjσij2-2∑i=1nλiσ0i2-2μ∑i=1nλi-1,where L(λ1,λ2,…,λn,μ) 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 λl(l=1,2,…,n) and μ such that(12)∂L∂λl=0,∂L∂μ=0.Equation (12) can be rearranged to the following form:(13)∑i=1nλiσil2-σ0l2-μ=0,∑i=1nλi-1=0.From (13), weight factors for any unknown stress are calculated by following form(14)λ1λ2⋮λnμ=σ112σ122⋯σ1n2-1σ212σ222⋯σ2n2-1⋮⋮⋮⋮⋮σn12σn22⋯σnn2-111⋯10-1σ012σ022⋮σ0n21.Finally, the unknown stresses for m points can be obtained from(15)sasb⋮sm∗=σ112σ122⋯σ1n2-1σ212σ222⋯σ2n2-1⋮⋮⋮⋮⋮σn12σn22⋯σnn2-111⋯10-1n+1×n+1σa12σb12⋯σm12σa12σb12⋯σm12⋮⋮⋮⋮σan2σb12⋯σm121111n+1×mTs1s2⋮sn0n+1×1.


## 3. Ap-Adaptive Refinement Using Modified SPR Technique

Adaptive procedures are to implement iteration analysis based on distinctly different levels of space (h)- or function (p)-refinements in specified local region to achieve solutions having a certain degree of accuracy in an optimal fashion. Particularly the p-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 p-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 p-adaptive mesh refinement scheme and a posteriori error estimator for p-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 adaptivep-refinement procedure since the number of sampling Gauss points in each element is increased as the p-level becomes higher. The increment of the quadrature point has an effect on the stress norm. According to Z/Z, er 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 as(16)er=∫ΩS∗-SpTD-1S∗-SpdΩ,where S∗ is a column vector including the true stresses; Sp including the stresses interpolated by Lobatto shape functions with the displacements obtained by FEM; [D] is a constitutive matrix; Ω is a mesh domain. Here a smoothed continuous stress concept is applied for the true stresses S∗ 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:(17)r^=∫ΩZpTD-1ZpdΩ1/2.Thus, the relative percentage error can be defined as(18)ηΩ=er2er2+r^21/2.In case of the modified SPR method using ordinary kriging, however, the estimated exact stress field denoted byS∗ should be calculated at each iteration round of p-refinement. In other words, S∗ cannot be fixed in the whole process of p-adaptive refinement the same as h-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 p-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 Uex and Ufe are the strain energy, the rate of convergence for thep-version of FEM can be derived by the inverse theorem [23, 27] as(19)Uex-Ufe≤kNp2α,where Uex and Ufe are the exact strain energy estimated by the limit value and the approximate strain energy by FEM, α is the strength of singularity, and Np and k are the degrees of freedom for the polynomial order p and a constant which depends on the mesh, respectively. There are three unknowns Uex, k, and α in (19). By performing three successive extension processes, p-2, p-1, and p, which are in the asymptotic range, we have three equations for computing the unknowns. Cancelling α and k in (19), the following extrapolation equation can be derived as(20)log⁡UexL-Up/UexL-Up-1log⁡UexL-Up-1/UexL-Up-2=log⁡Np-1/Nplog⁡Np-2/Np-1,where Up, Up-1, and Up-2 are the strain energies when the polynomial orders are p, p-1, and p-2 and UexL represents the limit value in terms of estimated exact strain energy where Np, Np-1, and Np-2 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 er in (16) is large, the polynomial order should be increased to satisfy an acceptable level of accuracy(21)ρΩ=UexL-UpUexL1/2,where UexL is the exact global strain energy and Up is also the global strain energy calculated in the current p-adaptive mesh consisting of nonuniform p-distribution at a certain iteration number.


## 4. Numerical Analysis

The numerical example is shown in Figure3 that specifies the geometric definition and analysis conditions that are given by a=50 cm, t=1.0 cm, E=2×107 N/cm2, ν=0.3, and σ=10 N/cm2. Due to the symmetry, a quarter of plate is modeled by three elements with fourth-degree polynomials for shape functions of p-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.Figure 3
Plate with a square cutout.
(a)
L-shaped domain
(b)
Finite element mesh byp-FEMFigure 4
Raw data fromp-FEM analysis when p-level = 4.In Figure5, the experimental semivariogram has been plotted with respect to the separation distance h. The allowable limit of separation distance Δh is assumed by 25% of h. Thus, the separation distance h can be allowed to use the similar distance denoted by h±Δh. Based on experimental semivariogram, three different theoretical semivariogram modules as shown in Table 1 have been tested. In this study, the Gauss model with Δh=25% of h has been adopted to find the weight factor for the OK process explained in (2).Figure 5
Comparison of semivariogram models usingΔh/h=0.25 when p-level = 4.Before the further analysis ofp-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 p-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).Figure 6
Stress contours by different interpolation techniques whenp-level = 4.
(a)
Raw data
(b)
OK interpolation
(c)
LSMFigure 7
Stress plots by different interpolation techniques whenp-level = 4.
(a)
Raw data
(b)
OK interpolation
(c)
LSMTo illustrate the applicability of OK interpolation to thep-adaptive mesh refinement, two p-version models are considered by 3-element and 12-element model. The initial p-level of both models begins with one. For nonuniform p-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 p-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.Figure 8
The finalp-adaptive meshes by 3-element model.
(a)
OK interpolation
(b)
LSMFigure 9
The finalp-adaptive meshes by 12-element model.
(a)
OK interpolation
(b)
LSMThe relative percentage errors have been illustrated in Tables2 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 p-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 p-adaptive model and 8 by OK basedp-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.Table 2
The relative percentage errors of 12-element model by the modified SPR method.
Number of iterations
Least square method
Ordinary kriging interpolation
NDF
η
Ω (%)
NDF
η
Ω (%)
12345678910
3689135190224272282314343363
20.4913.828.357.897.216.226.196.297.505.48
3696149214234264301353——
33.1215.117.114.004.303.013.544.17——Table 3
The relative percentage errors of 12-element model by the limit value.
Number of iterations
Least square method
Ordinary kriging interpolation
NDF
ρ
Ω (%)
NDF
ρ
Ω (%)
12345678910
3689135190224272282314343363
33.2615.6611.999.518.056.816.306.005.785.60
3696149214234264301353——
33.2613.3310.507.806.906.205.605.20——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 Table2. 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 Uex=7.7788×10-3 from Figure 10 that has been used to estimate the relative percentage error at each iteration round. In other words, in (21), UexL is the exact global strain energy, and Up is the global strain energy calculated in the current p-adaptive mesh consisting of nonuniform p-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 p-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.Figure 10
Estimation of exact strain energy by limit value approach.Figure 11
Convergence characteristics ofp-adaptive model using modified SPR method.
(a)
3-element model
(b)
12-element modelFigure 12
Convergence characteristics ofp-adaptive model using limit value approach.
(a)
3-element model
(b)
12-element model


## 5. Conclusions

The newp-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 Z/Z 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.



---

*Source: 1790256-2017-08-01.xml*