Element-Free Galerkin Method Based on Block-Pulse Wavelets Integration for Solving Fourth-Order Obstacle Problem

Abstract

We introduce improved element-free Galerkin method based on block pulse wavelet integration for numerical approximations to the solution of a system of fourth-order boundary-value problems associated with obstacle, unilateral, and contact problems. Moving least squares (MLS) approach is used to construct shape functions with optimized weight functions and basis. Numerical results for test problems are presented in this article to elaborate the pertinent features for the proposed technique. Comparison with existing techniques shows that our proposed method based on integration technique provides better approximation at reduced computational cost.

1. Introduction

Variational inequality theory has turned out to be an effective and potential tool for solution of mathematical models for obstacle, unilateral and contact problems. These problems arise in study of fluid flow through porous media, elasticity, transportation, optimal control, nonlinear optimization, and operations research [1–4]. The area of obstacle problems arising in fluid flow through porous media and elasticity forms an significant basis for the applications of variational inequalities. It has been revealed by Kikuchi and Oden [5] that the problem of equilibrium of elastic bodies in contact with a rigid frictionless foundation can be studied in the structure of variational inequalities. In a variational inequalities formulation, the location of the free boundary becomes an intrinsic part of the solution, and no special techniques are needed to locate it. Various numerical methods were employed to solve such type of problems [1–4]. If the obstacle is known, then the variational inequalities can be characterized by a system of differential equations by using the penalty function method of Lewy and Stampacchia [6]. In this paper, we apply element-free Galerkin method based on block pulse functions to develop a numerical method for obtaining smooth approximations to the solution of a system of fourth-order obstacle boundary-value problems of the type as follows: 𝑢⁽⁴⁾=𝑓(𝑥)𝑔(𝑥)𝑢(𝑥)(𝑥)𝑓(𝑥),𝑎≤𝑥≤𝑐,+𝑓+𝑟,𝑐≤𝑥≤𝑑,,𝑑≤𝑥≤𝑏. (1) With boundary conditions, consider 𝑢(𝑎)=𝑢(𝑏)=𝛼₁,𝑢^′′(𝑎)=𝑢^′′(𝑏)=𝛼₂,𝑢(𝑐)=𝑢(𝑑)=𝛽₁,𝑢^′′(𝑎)=𝑢^′′(𝑏)=𝛽₂, (2) where 𝑓 and 𝑔 are continuous functions on [𝑎,𝑏] and [𝑐,𝑑], respectively. The parameters 𝑟, 𝛼_𝑖, and 𝛽_𝑖 (𝑖=1,2) are real constants and the continuity conditions of 𝑢,𝑢^′′ at 𝑐 and 𝑑.

This problem was solved by many authors using different order of polynomial and nonpolynomial spline. Papamichael and Worsey [4] and Khalifa and Noor [3], have developed first- and second-order methods for solving such type of boundary value problems. Siraj-ul-Islam et al. [7] have established and analyzed smooth approximations for second, third, and fourth order linear/nonlinear boundary-value problems.

Galerkin-finite element method is one of the most popular well-developed methods, has been extensively used in the computational mechanics, due to its robustness, versatility, and convenience. However, the finite element method has inherent problems such as locking and poor derivative solutions [8]. FEM requires considerable amount of labor and time in generation of a predefined mesh/remesh of the problem domain. In structural analysis, the stresses calculated employing FEM algorithms are discontinuous and less accurate. When large deformations are analyzed, distortions in elements cause loss in solution accuracy in the case of large deformations [9]. Other numerical techniques such as the BEM (boundary element method) is limited by the availability of the infinite space fundamental solution for (at least the highest linear) differential operator of the problem.

The physical phenomena modeled by fourth-order obstacle boundary value problems are complicated, particularly, when they involve complex materials (e.g., viscoelastic materials with poisson ratio nearly 0.5) and where load is highly oscillatory. In conventional displacement-based approaches like finite element method (FEM), the interpolation of displacement requires 𝑐¹ continuity to ensure the convergence of FEM procedure for 4th order boundary value problems. This interpolation involves very intricate shape functions. These shape functions occupy large number of degrees of freedom in every element, including nodal rotations, displacement and so forth. For detailed discussion for solving such problems, the reader is referred to works [10] and references therein.

Improved element-free Galerkin (EFG) method provides an alternative choice of numerical tools which presents an attractive option in solution of such problems due to its flexibility in the absence of nodal interconnectivity. A meshfree (or meshless) method is a method used to establish system algebraic equations for the whole problem domain without the use of a predefined mesh for the domain discretization [9]. Several meshless techniques have been developed, for example, the “diffuse element method” [11] the “element-free Galerkin (EFG)” method [12] “meshless local Petrov-Galerkin (MLPG)” method [8], the “smooth particle hydrodynamics (SPH)” method, the “reproducing Kernel particle” method [13] “boundary node” method [14, 15] “boundary point interpolation” methods [16] and so forth [8].

2. Numerical Technique

Meshfree method has been applied to physical problems using Gaussian weight function. Moving least square method is used to generate the shape functions. Galerkin weak form of the boundary value problem is formulated to give the system equations. Penalty/Lagrange method is used to impose the boundary conditions. Numerical integration of the system equation has been performed using block-pulse function.

2.1. Meshless Shape Function

Shape functions for meshless techniques need to satisfy certain conditions such as adherence to partition of unity (∑^𝑛_𝑖=1Φ_𝐼(𝑥)=1), compact domain of influence, and adapting to randomness of nodes, to name a few. Computational efficiency is significantly affected by choice of shape function.

MLS. Moving least square (MLS) scheme was developed by mathematicians for interpolation, data fitting, and surface construction. Meshless techniques maintain the local character of the numerical implementation, by using a local interpolation to represent the trial function with the fictitious values of the unknown variable at some randomly located nodes. The local interpolation techniques generally used are MLS, partition of unity method, reproducing kernel particle method, hp clouds, Shepard function, and so forth. The moving least squares (MLSs) approximation has reasonably high accuracy and can be generalized to work with 𝑛-dimensional problems. It provides continuous approximation for the field function over the entire problem domain [9, 17].

Here, moving least square (MLS) scheme is considered to interpolate the field variable. MLS approximation 𝑢^ℎ(𝑥) of a field variable function 𝑢(𝑥) is defined as 𝑢^ℎ(𝑥)=^𝑚_𝑗=0𝑝_𝑗(𝑥)𝑎_𝑗(𝑥)=𝑝^𝑇(𝑥)𝑎(𝑥)∀𝑥∈Ω, (3) where 𝑚 = number of terms in the basis, 𝑝_𝑗(𝑥) = polynomial basis function,  𝑎_𝑗(𝑥) = nonconstant coefficients, and 𝑝^𝑇(𝑥)=[𝑝₁(𝑥),𝑝₂(𝑥),…,𝑝_𝑚(𝑥)].

𝑝^𝑇(𝑥) is a complete monomial basis vector of order 𝑚; 𝑎(𝑥) is a vector containing coefficients: 𝑎_𝑗(𝑥), 𝑗=1,2,3,…,𝑚. The unknown coefficients of approximation are computed by minimizing the difference between the local approximation at a point 𝑥_𝐼 and the nodal parameter 𝑢_𝐼 for the node 𝐼: 𝑢_𝐼=𝑢(𝑥_𝐼).

Here the sample point 𝑥_𝐼 may be a nodal point under consideration or a quadrature point. The support of the nodal point 𝑥_𝐼 is usually taken to be a circle of radius 𝑟_𝐼, centered at 𝑥_𝐼. The weighted residual functional is evaluated as 𝐽=^𝑛_𝑖=1𝑤𝑥−𝑥_𝐼𝑢^ℎ(𝑥)−𝑢_𝐼², (4) where 𝑤(𝑥−𝑥_𝐼) is weight function associated with the node 𝐼 calculated at point 𝑥. 𝑤(𝑥−𝑥_𝐼)>0 for all points 𝑥 in the support domain of node 𝐼, and 𝑢^ℎ𝑥_𝐼=𝑝^𝑇𝑥_𝐼𝑎(𝑥). (5) Minimization of weighted residual functional results in 𝐀(𝑥)𝐚(𝑥)=𝐁(𝑥)𝐮, (6) where 𝐀(𝑥) is a non-singular matrix. 𝐚(𝑥)=𝐀⁻¹(𝑥)𝐁(𝑥)𝐮,𝐮=𝑢₁,𝑢₂,…,𝑢_𝑛^𝑇. (7) The weight moment matrix 𝐀(𝑥) is 𝐀(𝑥)=𝐏^𝑇𝐖𝐏=^𝑛_𝑖=1𝑤𝑥−𝑥_𝐼𝑝𝑥_𝐼𝑝^𝑇𝑥_𝐼,𝐁(𝑥)=𝐏^𝑇𝐖, (8) where 𝐏 is 𝑝𝐏=^𝑇𝑥₁𝑝^𝑇𝑥₂⋯𝑝^𝑇𝑥_𝑛 (9)𝐖 is the weight function matrix: 𝑤𝐖=𝑥−𝑥₁⋯0⋮⋱⋮0⋯𝑤𝑥−𝑥_𝑛. (10) The following Gaussian weight function is used in the present study:𝑤𝑥−𝑥_𝐼=exp−𝛼_𝑐𝑥−𝑥_𝐼². (11) Shape perimeter (𝛼_𝑐) substantially affects the accuracy of the solution. Weight function for different values of shape parameters and its derivatives are illustrated in Figures 1 and 2 respectively.

Figure 1 

Weight function in the support domain.

Figure 2 

Derivatives of weight function in the support domain.

Obtain value 𝑎(𝑥) from (7), and substituting in (5) 𝑢^ℎ(𝑥)=Φ^𝑇(𝑥)⋅𝐮=^𝑛_𝐼=1𝜙𝑥_𝐼𝑢_𝐼∀𝑥∈Ω, (12)Φ(𝑥) is the nodal shape function Φ(𝑥)=𝜙₁(𝑥),𝜙₂(𝑥),…,𝜙_𝑛(𝑥). (13) For any node 𝐼, Φ_𝐼(𝑥)=𝑝^𝑇(𝑥)𝐴⁽⁻¹⁾(𝑥)𝐵(𝑥)_𝐼. (14) Transform the strong form of obstacle boundary value problems into symmetric variational or weak form as follows: ^𝑏_𝑎𝑎(𝑥)𝑑⁴𝑢𝑑𝑥⁴+𝑔(𝑥)𝑢(𝑥)−𝑓(𝑥)𝑣𝑑𝑥=0, (15) where 𝑢 is the trial function and 𝑣 is the test function. We get the following linear equation system after employing “Integration by parts” technique and by using the shape functions and boundary conditions ^𝑏_𝑎𝑎(𝑥)Φ^′′Φ^′′𝑇+𝑘(𝑥)ΦΦ^𝑇𝑑𝑥𝐮=^𝑏_𝑎𝐾𝑓Φ𝑑𝑥,𝐮=𝐹,𝐾=^𝑏_𝑎𝑎(𝑥)Φ^′′Φ^′′𝑇+𝑘(𝑥)ΦΦ^𝑇𝑑𝑥,𝐹=^𝑏_𝑎𝑓Φ𝑑𝑥. (16) Solving the system of equation, we get that 𝐮 nodal fictitious values are obtained. Required numerical solution can be determined by multiplication of nodal fictitious values and shape functions Φ: 𝐮=𝐮Φ. (17)

2.2. EFG

Element-free Galerkin (EFG) technique was introduced by Belytschko et al. [12]. EFG results in generation of nodes with variable interconnectivity, thus providing huge flexibility in discretization. Solution is expanded on a basis of MLS shape functions instead of piecewise polynomial basis as in case of the Finite element method. System equations are formulated using Galerkin weak form whereas integration is performed using background cells. EFG is conforming due to the use of consistent and compatible shape functions [9].

2.3. Block-Pulse Function-Based Integration

The accuracy of numerical approximation is dependent on the integration technique. Quadrature rule is commonly used for numerical integration. In this method, interpolating polynomials are employed to find weights for nodes in the domain. Conventional Gaussian quadrature scheme involves large number of points for accurate numerical solution. In other numerical integration techniques such as Newton-Cotes quadrature rule, fine discretization leads to erroneous results due to high degree of polynomial interpolation [18]. Accuracy of Gaussian quadrature rule depends on the selection of nodes and corresponding weights. Method of undetermined coefficients may be used for quadrature rule which requires tabulation of nodes and weights before numerical integration.

In order to avoid these difficulties, Block-Pulse functions (BPF) have been used in this study to find numerical solution of system integral equations. Wavelet transforms were primarily developed for signal analysis, but they have also been used in applications like image compression, data compression, de-noising data and numerical approximation. Various types of wavelets include Battle-Lemarie, B-spline, Chebyshev, Daubechies, Legendre and Block-Pulse functions [18].

Numerical Integration based on wavelets is becoming popular in recent years [18]. In this study, Block-Pulse functions have been used due to their suitable properties and exceptional performance in numerical integration [18]. Their performance particularly, in approximation of highly oscillatory and improper integral is more accurate as compared to existing method in literature [18–23].

Block-Pulse wavelet is one of the simplest wavelet. BPF with 𝑚 number of sets is defined as [18]: 𝜙_𝑖(𝑡)=()𝑇1for𝑖−1𝑚≤𝑡<𝑖𝑇𝑚0otherwise. (18) Here 𝑡∈[0,𝑇] and 𝑖=1,2,…,𝑚 and ℎ=𝑇/𝑚.

Let us consider the integrand ∫^𝑏_𝑎𝑔(𝑥)𝑑𝑥, for 𝑥=(𝑏−𝑎)𝑡+𝑎, we get: ^𝑏_𝑎𝑔(𝑥)𝑑𝑥=(𝑏−𝑎)¹₀𝑔((𝑏−𝑎)𝑡+𝑎)𝑑𝑡. (19)

Theorem 1. The integral can be approximated as [18] ¹₀𝑔(𝑡)1𝑑𝑡≈𝑚^𝑚_𝑖=1𝑔_𝑖. (20)

Proof. Consider ¹₀𝑔(𝑡)𝑑𝑡≈^𝑚_𝑖=1𝑔_𝑖¹₀𝜙_𝑖(𝑡)1𝑑𝑡=𝑚^𝑚_𝑖=1𝑔_𝑖. (21) To calculate the coefficients 𝑔_𝑖 we consider the nodal points, 𝑡_𝑘=2𝑘−12𝑚,𝑘=1,2,…,𝑚. (22) The equation can be discretized to give 𝑔𝑡_𝑘=^𝑚_𝑖=1𝑔_𝑖𝜙_𝑖𝑡_𝑘=𝑔_𝑘,𝑘=1,2,…,𝑚. (23) Thus the integral can be approximated based on 𝑚-set block-pulse function as ¹₀𝑔(𝑡)1𝑑𝑡≈𝑚^𝑚_𝑖=1𝑔2𝑖−12𝑚 (24) or ^𝑏_𝑎𝑔(𝑥)𝑑𝑥≅𝑏−𝑎𝑚^𝑚_𝑖=1𝑔𝑎+(𝑏−𝑎)2𝑖−1𝑚. (25) Optimal weights are calculated using built-in procedure in terms of wavelets. There is no need for tabulation of weights, and no intermediate technique is required for integration. BPF-based integration eliminates the need for quadrature nodes. Integrands are treated explicitly without solving nonlinear system which results from unknown nodes and weights.

2.4. Consistency

A numerical method is called convergent if it approaches the exact solution as the discretization is refined. Consistency and stability guarantee convergence. The choice of monomial basis vector 𝑝(𝑥) influences the consistency of the EFG method. MLS approximation reproduces all components that appear in 𝑝(𝑥). EFG method is consistent with order 𝑘 if all monomials up to order 𝑘 are included in 𝑝(𝑥) [24]. If 𝑘=1, this implies that 𝑝(𝑥) is linear complete and if 𝑘=2, 𝑝(𝑥) is quadratic complete, and so forth. For consistency test we consider MLS approximation which reproduces all components of polynomials. Consider ^𝑛_𝐼=1Φ^𝑘_𝐼(𝑥)𝑢𝑥_𝐼=𝑢^𝑘(𝑥)(𝑥)∀𝑥∈Ω,If𝑢=1,then^𝑛_𝐼=1Φ_𝐼(𝑥)(𝑥)=1,If𝑢=𝑥,then^𝑛_𝐼=1Φ_𝐼(𝑥)𝑥_𝐼=𝑥. (26)

3. Application and Test Problems

We consider the following example fourth-order obstacle boundary value problem (BVP) for finding 𝑢 such that 𝑢⁽⁴⁾(𝑥)(𝑥)[],𝑢(𝑥)(𝑥)[],𝑢≥=𝑓,onΩ=−1,1≥=𝜓,onΩ=−1,1₍₄₎(𝑥)(𝑥(𝑥[],𝑢()(1)−𝑓)𝑢−𝜓)=0,onΩ=−1,1−1=𝑢=0,𝑢^′′(−1)=𝑢^′′(1)=0, (27) where 𝑓, a given force acting on the beam and 𝜓(𝑥), is the elastic obstacle. Now the problem will be discussed in the framework of variational inequality approach. Define set 𝐾, as 𝐾=𝑣∶𝑣∈𝐻²₀(Ω)∶𝑣≥𝜓,onΩ. (28)𝐾 is a closed convex set in 𝐻²₀(Ω), (Sobolev space), which is in fact a Hilbert space. It can be easily shown that energy functional associated with the obstacle problem (27) is 𝐼[𝑣]=¹₋₁𝑑⁴𝑢𝑑𝑥⁴−2𝑓(𝑥)𝑣(𝑥)𝑑𝑥,∀𝑣∈𝐻²₀(Ω),𝐼[𝑣]=¹₋₁𝑑²𝑢𝑑𝑥²²𝑑𝑥−2¹₋₁𝑓(𝑥)𝑣(𝑥)𝐼[𝑣](),𝑎()𝑑𝑥,=𝑎𝑢,𝑣−2𝑓,𝑣𝑢,𝑣=¹₋₁𝑑²𝑢𝑑𝑥²𝑑²𝑣𝑑𝑥²=𝑑𝑥,𝑓,𝑣¹₋₁𝑓(𝑥)𝑣(𝑥)𝑑𝑥, (29)(30)(31) it can easily proved that the form 𝑎(𝑢,𝑣), defined by (30), is bilinear, symmetric, and positive (coercive) and the functional 𝑓, is a linear continuous functional [1–3]. Minimum 𝑢 of the functional 𝐼[𝑣] defined by (29) on the closed convex set 𝐾 in ∈𝐻²₀(Ω) can be characterized by the variational inequality 𝑎(𝑢,𝑣−𝑢)≥𝑓,𝑣−𝑢,∀𝑣∈𝐾. (32) Thus we conclude that the obstacle problem (27) is equivalent to solving variational inequality problem (32). This equivalence has been used to study the existence of a unique of (27); see [1–3], and using the idea of Lewy and Stampacchia, [6] can be written as 𝑢⁽⁴⁾+𝑣𝑢−𝜓𝑢−𝜓=𝑓, (33) where 𝜓 is the obstacle function and 𝑣(𝑡) is a penalty function defined by 𝑣(𝑡)=4𝑡≥00𝑡≺0. (34) We assume that the obstacle function, 𝜓(𝑥), is defined by 𝜓(𝑥)=−141for−1≤𝑥≤−2,121≤𝑥≤14−121≤𝑥≺2. (35) From (33)–(35), we obtain the following system of equations 𝑢⁽⁴⁾=1𝑓for−1≤𝑥≤−2,121≤𝑥≤11−4𝑢+𝑓−21≤𝑥≺2. (36) With the boundary conditions, consider 𝑢()−1−1=𝑢21=𝑢2(1)𝑢=𝑢=0,_′′(−1)=𝑢^′′−12=𝑢^′′12=𝑢^′′(1)=0. (37)

Problem 1. In this example, we consider system of differential equation (38), when 𝑓=0,𝑢⁽⁴⁾=10for−1≤𝑥≤−2,121≤𝑥≤11−4𝑢−21≤𝑥≺2,𝑢()−1−1=𝑢21=𝑢2(1)𝑢=𝑢=0,_′′(−1)=𝑢^′′−12=𝑢^′′12=𝑢^′′(1)=𝜀,𝜀⟶0. (38) The analytical solution for this boundary value problem is 𝑢(𝑥)=−23𝑥²−32𝑥²−13−112𝑥4𝜀1−1≤𝑥≤−210.5−𝜙₃𝜙₁sin(𝑥)sinh(𝑥)+𝜙₂(𝑥)(𝑥−1coscosh)21≤𝑥≤2−23𝑥²+32𝑥²−13+112𝑥4𝜀12≤𝑥≤1, (39) where 𝜙₁=sin(1/2)sinh(1/2),  𝜙₂=cos(1/2)cosh(1/2),  𝜙₃=cos(1)+cosh(1). Numerical and analytical solution of Problem 1 is shown in Figure 3.

Figure 3 

Numerical and analytical solution plot (16 nodes), with constant basis-Problem 1.

Problem 2. For 𝑓=1, 𝑢⁽⁴⁾=11for−1≤𝑥≤−2,121≤𝑥≤12−4𝑢−21≤𝑥≺2. (40) With the same boundary conditions as in previous example, consider 𝑢(𝑥)=124𝑥⁴+18𝑥³+18𝑥²+3+164𝑥1192−1≤𝑥≤−210.5−𝜙₃𝜙₁sin(𝑥)sinh(𝑥)+𝜙₂(𝑥)(𝑥−1coscosh)21≤𝑥≤2124𝑥⁴−18𝑥³+18𝑥²−3+164𝑥11922≤𝑥≤1. (41) Numerical and analytical solution of Problem 2 is shown in Figure 5.

Test problems were solved for different number of nodes and constant basis with the help of improved EFG based on numerical integration using Block pulse functions/wavelets (BPF). Figures 4 and 6 show the errors in solution of test Problems 1 and 2, respectively. Maximum absolute errors in solution are tabulated and compared with the results of exact solution, which show better accuracy.

Figure 4 

Max. absolute error plot (16 nodes), constant basis-Problem 1.

Figure 5 

Numerical and analytical solution plot (12 nodes), with constant basis-Problem 2.

Figure 6 

Max. absolute error plot (12 nodes), constant basis-Problem 2.

3.1. Results and Discussion

Maximum errors obtained using our technique based on block pulse function (BPF) for our test problems are shown in Tables 1 and 2. Comparison of maximum absolute errors with exact solution demonstrates the accuracy of our method. Results using quadrature rule-based integration are also calculated to show the performance of BPF-based integration technique. Maximum absolute errors for test problems are shown in Tables 1 and 2. It is clearly seen that integration of system equations using BPF has produced better results as compared with quadrature rule method for numerical integration.

Table 1 

Observed maximum errors for Problem 1, using BPF (25).

Table 2 

Observed maximum errors for Problem 2, using BPF (25).

Effect of Integration Points. It has been observed, that the accuracy of meshless computational technique depends upon the size of influence of domain, penalty factor, and Integration points. In this paper, we have studied that the accuracy of numerical method generally improves as we increase the number of integration points; however, if the number of gauss points exceeds some limiting value, it may have some adverse effects on the accuracy. It has been determined that arrangement of background cells and ratio of total number of Gauss points to the total number of nodes, plays very important roles in accuracy. To avoid all these complications to improve the accuracy of the meshless method, we purpose new numerical integration technique as mentioned above.

4. Conclusion

In this study, improved numerical technique “element free Galerkin (EFG) technique using Block-Pulse function-based integration” has been presented. The results obtained by the suggested method for given problems as mentioned ((38) and (40)) exhibit its ability to provide improved solutions. The improved EFG technique showed fast convergence and provided better results at reduced number of nodes, as mentioned in Tables 1 and 2. It has capability of handling of improper, highly oscillatory function. Details of the results are shown the figures and tables.