# Adaptive meshless local Petrov-Galerkin method with variable domain of influence in 2D elastostatic problems.

IntroductionIn recent years, meshless methods have been developed as alternative numerical approaches in efforts to eliminate known drawbacks of the Finite Element Method (FEM). The main objective in developing meshless methods was to eliminate, or at least reduce, the difficulty of meshing and remeshing of complex structural elements. The nature of various approximation functions employed by meshless methods allows the definition of problem domains by simply adding or deleting nodes where desired. Nodal connectivity to form an element as in FEM method is not needed, only nodal coordinates and their domain of influence (DOI) are necessary to descretize the problem domain. Meshless methods may also reduce other problems associated with the FEM, such as solution degradation due to locking and severe element distortion [1].

One of these meshless method is the Meshless Local Petrov-Galerkin (MLPG). This method is believed to have a good future due to its generality in choosing the form of test and trial functions and also that it is similar to the well established Element Free Galerkin (EFG) method. Atluri et al. [1] proposed a new integration method in a local domain, based on a Local Symmetric Weak Form (LSWF).

Therefore, the MLPG method is a truly meshless method, and all other meshless methods can be derived from it, as special cases, if trial and test functions and the integration method are chosen appropriately.

Pudjisuryadi and Barry [2] presented MLPG method using polygonal sub-domains constructed from several triangular patches rather than the typically used circular sub-domains. Moving least-squares (MLS) approximation is used to construct the trial displacements and linear, Lagrange interpolation functions are used to construct the test functions. An adaptive technique to improve the accuracy of approximate solutions is developed using the effective stress gradient as the error indicator. The nodal sub-domains were held constant throughout the subsequent adaptive analysis.

In the adaptive technique, problem domain will be refined by placement of new nodes in area which local error exceeds the prescribed level. If the nodal configuration becomes very dense in an area, and domain of influence size is set to be constant, too many nodes will be taken into account in approximating the trial function. This will result in higher computational cost and in-effective adaptive technique (local characteristic of the approximation is lost). To alleviate this problem, variable domain of influence is proposed in this study.

The standard MLPG formulation and MLS approximation can be seen in some references [1,2,3,4], while the effective stress gradient error indicator, variable domain of influence, numerical results and conclusion are discussed in these following sections.

Effective Stress Gradient Indicator

The measure used in this work to determine if a patch exceeds a certain error level is the Effective Stress Gradient Indicator [5], EK, which is calculated as:

[E.sub.K] = [h.sub.K] [g.sub.K] (1)

where [h.sub.K] is the characteristic size of the patch (the edge of the triangle patch with the smallest length), and

[g.sub.K] = max |[[sigma].sub.e] ([P.sub.K]) - [[sigma].sub.e] ([P.sub.Ki])/d([P.sub.K][P.sub.Ki]) (2)

with

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (3)

is the effective von Mises stress at the center of the patch, and d([P.sub.K][P.sub.Ki]) is the distance from the central point of the patch ([P.sub.K]) to the central point of the ith adjacent patch ([P.sub.Ki]), as shown in Figure 1. While [[sigma].sub.x], [[sigma].sub.y], [[sigma].sub.z], [[tau].sub.xy], [[tau].sub.xz] and [[tau].sub.yz] are the normal and shear stresses in a three dimensional body.

[FIGURE 1 OMITTED]

Variable Domain of Influence

The MLS approximation of trial function in an arbitrary location in a problem domain is well defined, if there are at least a minimum number of nodes which domains of influence cover that location of interest (has non-zero weight function on that location). Such minimum number of nodes depends on the polynomial basis function used in the MLS approximation. In 2D spatial coordinates [x, y], linear and quadratic approximation of trial functions are as follows:

u(x) = [a.sub.0] + [a.sub.1]x + [a.sub.2]y (4)

u(x) = [a.sub.0] + [a.sub.1]x + [a.sub.2]y + [a.sub.3][x.sup.2] + [a.sub.4]xy + [a.sub.5][y.sup.2] (4)

where u(x) is the approximation of displacement at an arbitrary location x, while ai are the constants that define the approximation. From Equations 4 and 5, it can be explained, that three and six nodes are the minimum requirements to define the approximation constants in linear and quadratic basis functions respectively. In order to make the MLS approximation effective throughout the subsequent adaptive, the number of nodes to determine the constants should not be too few or too many. In this study, the minimum number of nodes is set to be four times the minimum requirements (12 and 24 nodes for linear and quadratic basis functions respectively), thus the size of domain of influence will automatically be variable depending on the nodal configuration. Figure 2 shows the typical variable domain of influence used in MLS approximation with linear basis function. In the figure, black dots show the nodes, while crosses indicate the point where MLS approximation is evaluated. It can be seen that 12 and 15 nodes are inside the domain of influence of first and second point of interest (full line circle). The domain of influence used for second point of interest is forced to be larger, since smaller circle (dashed line circle) does not include enough number of nodes (only include 11 valid nodes since four nodes a, b, c, and d on the perimeter of the circle have zero weight function). Moreover, for non-convex problem domain, visibility criterion is employed [6], to take into account the material discontinuity. One of the simplest visibility criterion, that is the influence of a node to a location that is not visible is neglected (weight function is set to zero), is used in this study. Visibility is a condition where a straight line connecting the node and the point of interest does not cross the problem boundary. Domain of influence of such case can be seen in Figure 3. Again, cross sign indicates the point of interest, while black dots show the nodes. It can be seen that dashed line circle has already covered 14 nodes, but two of them (nodes d and e) have zero weight function and another two (nodes b and c) fail the visibility criterion, resulting in only 10 valid nodes. Thus, the domain of influence should be increased (full line circle) to include enough number of valid nodes. Note that in the new domain of influence, another node (a) should be neglected.

[FIGURE 2 OMITTED]

[FIGURE 3 OMITTED]

Patch Tests and Numerical Examples

Patch test is a procedure to assure that an element formulation is stable. It is originally used in the finite element method. Patch test is a simple test that can be performed numerically to check the validity of an element formulation and its implementation in a program. When an element formulation is used to model a structure, mesh refinement will produce a sequence of approximate solutions that converges to the exact solution [7].

[FIGURE 4 OMITTED]

This procedure is adopted to check the validity of the proposed meshless method. Linear, quadratic, and cubic patch tests, with configuration as seen in Figure 4, are done to verify the proposed method. Convergence is expected as finer nodal configuration is used to model a structure. Moreover, exact solution should be achieved if the basis function used in the MLS approximation has the same or higher degree than the patch test problem (linear patch tests which are solved using linear or quadratic basis functions, and quadratic patch test which is solved using quadratic basis function).

A unit square is chosen to represent the problem domain. The boundary conditions applied in each of the three cases as well as the prescribed tractions can be seen in Figure 4. The exact solutions of the displacement fields for linear, quadratic and cubic patch tests are given in Equations 6, 7, and 8 as follows:

u(x, y) = x/E & v(x, y) = -vy/E (6)

u(x, y) = xy/E & v(x, y) = -[x.sup.2]/2E - [vy.sup.2]/E (7)

u(x, y) = [x.sup.2]y/2E + [vy.sup.3]/6E - [y.sup.3] (1 + x)/3E &

v(x, y) = -[vxy.sup.2]/2E - [x.sup.3]/6E (6)

where:

u(x, y) = displacement component in x direction.

v(x, y) = displacement component in y direction.

E = modulus of elasticity.

v = Poisson ratio.

An initial 4x4 nodal distribution is used in each patch test and the EK error indicator limit is chosen to be 0.1 (F/[L.sup.2]). Linear, quadratic and cubic patch tests are named C1, C2 and C3, respectively. MLS approximation using linear and quadratic basis function are named B1 and B2 respectively. [L.sup.2] Norm of displacement error, [parallel] [e.sup.u] [parallel] calculated over the entire problem domain . is used to examine the convergence in patch tests ([parallel] [e.sup.u] [parallel] can be used only for problems which exact solutions are known), and it is computed as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (9)

in Equation 9, subscripts ex and app indicate exact and approximation and . indicates the problem domain (union of the triangular patches). The results of the tests are summarized in Table 1. It can be seen that the method proposed satisfy the patch tests, that is exact solution is found directly (practically zero value of [parallel] [e.sup.u] [parallel]) if the degree of basis function is the same or higher then the exact solution (C1B1, C1B2, and C2B2), while other cases (C2B1, C3B1, and C3B2) show convergence. It should be noted that high [E.sub.K] indicator does not always show poor solution, but it only shows the presence of high stress gradient. In case C2B2, exact solution is found in 1st step, but relatively high stress gradient (larger than [E.sub.K] limitation) is still detected since coarse nodal configuration is used. In this case, the analysis will continue to the next step until the prescribed error limit (=0.1 F/[L.sup.2]) is met.

[FIGURE 5 OMITTED]

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

[FIGURE 8 OMITTED]

[FIGURE 9 OMITTED]

[FIGURE 10 OMITTED]

The final nodal configurations of the patch tests after sequences of adaptive procedure can be seen in Figure 5 to Figure10. In the figures, the nodes are plotted in the deformed shape of the patch. The circle (o) and plus (+) signs show the approximate and exact nodal displacement respectively. The two signs almost coincide to each other since the [L.sup.2] Norm of displacement error, [parallel] [e.sup.u] [parallel] in each test is small enough.

Two additional numerical tests are performed to further evaluate the proposed method. One is a cantilever beam subjected to point load at the tip, and the other is plate with center hole under subjected to uni-axial traction. The illustration of the two tests can bee seen in Figures 11 and 12. In these two additional tests, only linear basis function is used, since reasonably good reason still can be achieved with less computational effort (see Table 1).

A cantilever beam, as shown in Figure 11 is analyzed using an initial 6x4 uniform nodal distribution. The modulus of elasticity, E and the Poisson ratio, v are assumed as 1000 (F/[L.sup.2]) and 0.25, respectively, the maximum allowable effective stress gradient ([E.sub.K]) is chosen as 5.0 (F/[L.sup.2]), half of given uniform traction, P.

[FIGURE 11 OMITTED]

The analytical solution for the cantilever beam problem [8] is given in Equation 10,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (10)

A plate with a circular center hole, as illustrated in Figure 12, is analyzed using an initial 7x7 nodal distribution. The modulus of elasticity, E and Poisson ratio, v are assumed as 10 (F/[L.sup.2]) and 0.25, respectively, the maximum allowable effective stress gradient, [E.sub.K] is set as 0.5 (F/[L.sup.2]), half of given uniform traction, P.

[FIGURE 12 OMITTED]

[FIGURE 13 OMITTED]

Only the right, upper quarter of the plate is analyzed due to problem symmetry. The corresponding symmetry boundary conditions consist of rollers as shown in Figure 12 since the plate fibers on the axis of symmetries (x=0 and y=0) cannot move perpendicular to an axis of symmetry. The plate is subjected to a uniform load in the x direction as shown above. The quarter plate has a dimension of 4x4 units, and the quarter hole has a radius (a) of 1 unit (L). The analytical solution is given by Atluri et al. [9] in Equation 11:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (11)

where r is the distance from the center of the hole, and v is the angle measured from the positive x axis, counter-clockwise, as illustrated in Figure 12.

[FIGURE 14 OMITTED]

The summary of the test results are shown in Table 2. The final nodal configurations of the two tests are shown in Figures 13 and 14. In the figures, circle (o) and plus (+) signs show the approximate and exact nodal displacement respectively. The signs almost coincide to each other in the cantilever beam test. Noticeable displacement error is shown in the plate with a center circular hole test since the L2Norm of displacement error, [parallel] [e.sup.u] [parallel] is still relatively high. Of course, better result can always be obtained by reducing the effective stress gradient indicator, [E.sub.K] limit.

Conclusion

In this study, an adaptive meshless local Petrov-Galerkin (MLPG) method that employs polygonal sub-domains constructed from several triangular patches rather than the typically used circular subdomains is presented. Variable domain of influence and effective stress gradient error indicator are used in the proposed adaptive technique. Linear, quadratic and cubic displacement patch tests are done to verify the method. Two additional tests, a cantilever beam subjected to point load at the tip, and a plate with a hole subjected to uniform tension are also performed.

From the results of the five test cases, the achievements and conclusions of this work are summarized as follows:

1. Refinement patterns using the local effective stress gradient indicator ([E.sub.K]) for adaptive MLPG analysis takes place in the regions with high stress gradients, which usually also have the highest values of stress.

2. The local effective stress gradient indicator is applicable to all cases since it does not require any exact solutions. This is a very important advantage of the [E.sub.K] error indicator.

3. One drawback of the local effective stress gradient indicator is that adaptive refinement will occur indefinitely in regions which have an infinite stress solution (very high effective stress gradient). In order for the local [E.sub.K] indicator to decrease, the rate of decreasing effective stress difference must be faster than the rate of decreasing patch size (see Equation 1).

Note: Discussion is expected before November, 1st 2008, and will be published in the "Civil Engineering Dimension" volume 11, number 1, March 2009.

Received 31 May 2008; revised 1 August 2008; accepted 2 August 2008.

References

[1.] Atluri, S.N., Kim, H.G., Cho, J.Y, A Critical Assessment of the Truly Meshless Local Petrov-Galerkin (MLPG), and Local Boundary Integral Equation (LBIE) Methods, Computational Mechanics vol. 24, 1999, pp, 348-372.

[2.] Pudjisuryadi, P., Barry, W.J., An Adaptive Technique for 2D Elastostatic Analysis by the Meshless Local Petrov Galerkin Method, International Conference on Advancement in Design, Construction, Construction Management and Maintenance of Building Structures, Bali, 2002, pp, I-125-I-136.

[3.] Pudjisuryadi, P., Moving Least Squares Approximation to be used with Meshless Numerical Analysis Methods, Civil Engineering Dimension, vol.4, no.1, 2002, pp, 47-50

[4.] Pudjisuryadi, P., Introduction to Meshless Local Petrov-Galerkin Method, Civil Engineering Dimension, vol.4, no.2, 2002, pp, 112-116.

[5.] Moshfegh R., Li X., Nilsson L., Gradient-Based Refinement Indicators in Adaptive FEM Analysis with Special Reference to Sheet Metal Forming, Engineering Computation, vol.17, no.3, 2000, pp, 910-932.

[6.] Belytschko, T., Krongauz, Y., Organ, D., Fleming, M., Krysl, P., Meshless Methods: an Overview and Recent Developments, Computer Method in Applied Mechanics and Engineering 139, Nos. 14, 1996, pp, 3-47.

[7.] Cook, R.D., Malkus, D.S., Plesha, M.E., Concepts and Applications of Finite Element Analysis, third edition, John Wiley & Sons, New York, 1988.

[8.] Timoshenko, S.P., Goddier, J.N., Theory of Elasticity, third edition, McGraw-Hill Inc., New York, U.S.A., 1987.

[9.] Atluri, S.N., Zhu, T., New Concepts in Meshless Methods, International Journal for Numerical Methods in Engineering 47, 2000, pp, 537-556.

Pamuda Pudjisuryadi (1)

(1,) Department of Civil Engineering, Petra Christian University, Surabaya, Indonesia E-mail: pamuda@petra.ac.id

Table 1. Patch tests employing the EK indicator (0.1 F/L2)) and variable DOI Case Step Number of ||e|| (L) Max [E.sup.K] nodes/patches (F/[L.sup.2]) C1B1 1 16/18 6.24 x [10.sup.-14] 6.52 x [10.sup.-13] C1B2 1 16/18 1.38 x [10.sup.-15] 9.89 x [10.sup.-15] C2B1 1 16/18 2.55 x [10.sup.-03] 3.34 x [10.sup.-01] 2 49/72 3.51 x [10.sup.-04] 1.69 x [10.sup.-04] 3 169/288 6.30 x [10.sup.-05] 8.56 x [10.sup.-02] C2B2 1 16/18 2.34 x [10.sup.-15] 3.33 x [10.sup.-01] 2 49/72 6.66 x [10.sup.-13] 1.67 x [10.sup.-01] 3 169/288 7.57 x [10.sup.-10] 8.33 x [10.sup.-02] C3B1 1 16/18 2.27 x [10.sup.-03] 5.53 x [10.sup.-01] 2 49/72 4.69 x [10.sup.-04] 2.92 x [10.sup.-01] 3 152/261 1.52 x [10.sup.-04] 1.45 x [10.sup.-01] 4 323/587 1.86 x [10.sup.-04] 1.22 x [10.sup.-01] 5 345/630 2.47 x [10.sup.-04] 1.06 x [10.sup.-01] 6 350/640 2.43 x [10.sup.-04] 9.65 x [10.sup.-02] C3B2 1 16/18 1.17 x [10.sup.-03] 5.23 x [10.sup.-01] 2 49/72 3.57 x [10.sup.-04] 2.77 x [10.sup.-01] 3 152/261 3.06 x [10.sup.-04] 1.51 x [10.sup.-01] 4 329/598 2.75 x [10.sup.-04] 8.73 x [10.sup.-02] Table 2. Cantilever beam and plate with a center hole tests (Variable DOI) Case Step Number of e (L) nodes/patches Cantilever 1 24/30 2.18 x [10.sup.-03] Beam Test 2 75/117 7.71 x [10.sup.-04] (max [E.sub.K]=5.0) 3 209/358 4.43 x [10.sup.-04] 4 485/881 5.07 x [10.sup.-04] Plate With Center 1 48/70 1.67 x [10.sup.-01] Hole Test 2 147/252 8.99 x [10.sup.-02] (max [E.sub.K]=0.5) 3 199/350 8.49 x [10.sup.-02] 4 251/445 8.44 x [10.sup.-02] Max EK (F/[L.sup.2]) Case Cantilever 15.8 Beam Test 10.1 (max [E.sub.K]=5.0) 5.44 3.55 Plate With Center 11.0 Hole Test 0.80 (max [E.sub.K]=0.5) 0.55 0.37

Printer friendly Cite/link Email Feedback | |

Author: | Pudjisuryadi, Pamuda |
---|---|

Publication: | Civil Engineering Dimension |

Date: | Sep 1, 2008 |

Words: | 3298 |

Previous Article: | Fractional critical shear stress at incipient motion in a bimodal sediment. |

Next Article: | Modeling the influence of project manager trustworthy leadership behavior upon construction team trust. |