Tuesday, 20 April 2021

Lupine Publishers| A Reduced Lagrange Multiplier Method for Dirichlet Boundary Conditions in Isogeometric Analysis

 Lupine Publishers| Trends in Civil Engineering and its Architecture (TCEIA)


Abstract

Although the well-known standard Lagrange multiplier method (LMM) can even produce higher accuracy and easier implementation than other conventional schemes (e.g., the modified variational principle, the Nitsche's method), however it inherently owns many difficulties in solving the system of discretized equations, mainly caused by new unknown Lagrange multipliers. The LMM naturally increases the problem size and leads to a poorly conditioned matrix equation. The singularity is also often encountered because of inappropriate selection of interpolation space for the Lagrange multiplier. In this paper, we propose an improved method, called reduced Lagrange multiplier method, which can overcome such drawbacks raised by the LMM in treating the Dirichlet-type boundary conditions in terms of Isogeometric Analysis. By simply splitting the system equations into boundaries and interior groups, the size of system equations derived from the LMM is reduced; no additional unknowns have been added into the resulting system of equations; the Lagrange multiplier is hence disappeared; and more importantly the singular problem mentioned is avoided. The accuracy and convergence rates of the proposed method are studied through three numerical examples, exhibiting all the desirable features of the method. Optimal convergence rate and high accuracy for the present method is found.

Keywords: NURBS; Isogeometric analysis; Dirichlet boundary conditions; Lagrange multiplier method; Finite element method

Abbrevations: LMM: Lagrange Multiplier Method; IGA: Iso-Geometric Analysis; CAD: Computer-Aided DBCs: Design Boundary Conditions; DM: Direct Method; LSCM: Least-Squares Collocation Method; RLMM: Reduced Lagrange Multiplier Method

Introduction

The Isogeometric analysis (IGA) [1] employs the same computer- aided design (CAD) spline basis functions (e.g., NURBS, T-splines) to describe exact geometries and approximate physical responses. The IGA offers possibilities of integrating finite element analysis into CAD design tools and avoids mesh generation process and subsequent communication with a CAD model during refinement. The IGA has shown to be a highly accurate and robust approach for numerical simulation with exact geometry representation even with coarse mesh [2]. Moreover, the high-order smoothness of the NURBS basis functions allows for a direct construction of rotation- free plate/shell formulations [3,4] and is attractive for solving PDEs that incorporate fourth-order (or higher) derivatives of the field variable, such as the Hill-Cahnard equation [Gomez et al., 2008] [5] and gradient damage models [6].

The IGA has shown to be an effective approach in solving engineering problems as a large number of works on the theory and application of IGA have been conducted and reported in literature, including mathematical properties [7,8], integration method [9,10], spline techniques [11], fluid mechanics [12], fluid- structure interaction [13], plate/shell structures [14-17], electromechanic [18], magneto-electro elastic [19], shape optimization [20], fracture mechanics [21], and unsaturated flow problem [22]. Despite the efforts dedicated to the IGA in recent years and the significant progress has been made, some relevant technical aspects still require further studies; and among others the imposition of the Dirichlet-type essential boundary conditions (DBCs) in the IGA is an important technical issue. Given the non-interpolating nature of the NURBS basis, similar to moving least-squares [23], the properties of Kronecker delta are not satisfied, directly imposing DBCs is difficult, and special techniques are hence required.

Hughes et al. [1] addressed a direct method (DM) for DBCs enforcement to the control points. However, Wang & Xuan [24] found that the DM may induce significant errors with deteriorated convergence rates, and a special DBC treatment technique in terms of mesh less computation for contact problems [25] is presented, which is accomplished by separating the interior and boundary points on the basis of the mixed transformation method. They numerically illustrated that the improved approach produces favorable solution accuracy and optimal convergence rates compared with DM. However, the selection of boundary interpolation points should be emphasized in order to avoid a singular transformation matrix. Koo et al. [26] recently extended this method to iso geometric shape design.

Costantini et al. [27] applied the quasi-interpolation method to set inhomogeneous DBCs in the IGA on the basis of generalized B-splines to eliminate the solution of large global linear systems. Lu et al. [28] presented a technique to blend NURBS with Lagrangian interpolation in the essential boundary and allow the direct imposition of prescribed values. Other strategies based on a modified weak-form equation, such as Nitsche's method [29] and Lagrange multiplier method (LMM) have been adopted in NURBS-based IGA. Recently, De Luycker et al. [30] proposed a least- squares collocation method (LSCM) for imposing arbitrary DBCs in extended IGA. Govindjee et al. [31] later improved the LSCM by offering an efficient local least-squares algorithm to global least squares. However, LSCM is affected by the number of interpolating points. Thus, an important issue on how to choose the location and number of interpolating points precisely arises. The LMM is widely used in weak-form methods for DBCs enforcement [32].

The LMM offers solutions higher accuracy than the modified variational principle [33]. The implementation of LMM is simpler and more straightforward than the Nitsche's method [32]. However, the LMM introduces a new unknown Lagrange multiplier, which makes it difficult in selection of interpolation space for the Lagrange multiplier so that singular issues in solving the system of equations are avoided. In this paper, we present a new improved method, reduced Lagrange multiplier method (RLMM) for effective treatment of the DBCs in terms of IGA. Given that the NURBS control points can be clearly partitioned into boundary and interior groups and that the NURBS basis functions of interior control points vanish at the boundary [24], we divide the system equation of the LMM into boundary and interior groups.

Thus, the size of the equation system is now reduced and the shortcomings of the LMM can hence be overcome. Furthermore, the choice of the interpolation space for the Lagrange multiplier is studied numerically and is discussed in numerical results. The rest of this paper is organized as follows. Section 2 presents a brief on the NURBS-based IGA. Section 3 details the RLMM for essential boundary condition treatments in IGA. Numerical examples are considered in Section 4 to verify the accuracy and effectiveness of the present method. The conclusions are drawn in Section 5.

NURBS-Based Isogeometric Analysis

NURBS basis functions

In a 1D parametric space ξε[0,1], a knot vector κ(ξ) to construct the B-spline basis functions is a set of non-decreasing numbers written as Lupinepublishers-openaccess-journals-Civil-engineering-Architechture where and are the number of basic functions and the order of basic functions, respectively. Given a knot vector κ (ξ), the B-spline basis function of degree p, written as Ni,p (ξ), is defined recursively as follows [34] Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

and

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

The NURBS basis function Ri,p (ξ) in the framework of partition of unity is constructed by a weighted average of the B-spline basis functions [35]:

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

where is the ith weight.

Similarly, the bivariate NURBS basis functions are expressed in the following form:

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

Where represents the 2D weight; Ni,p (ξ) and Nj,q (η) are the B-sp line basis functions defined on knot vectors κ(ξ) and κ(η) , respectively, following the recursive formula shown in Eqs. (1) and (2).

The NURBS basis functions are non-interpolatory in the interior parametric space interval. Therefore, essential boundary conditions cannot be applied simply by enforcing individual values to the control points.

By using the NURBS basis functions, a NURBS surface of order in the direction and order in the direction can be constructed as follows:

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

where Pij are the coordinates of the control points in two dimensions.

Governing equations and discretization

Without loss of generality, we consider a 2D isotropic elastic problem whose problem domain is with a boundary ( and ), the corresponding governing equations are expressed as follows:

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

where ∇ denotes the gradient operators; ó is stress tensor; b and t are the body and surface force vectors, respectively; n is the normal vector on Neumann boundary Γt ;g is the known displacement on Γu .

The total potential energy of the system can be constructed as follows:

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

In the NURBS-based IGA, an isoparametric concept is adopted and the NURBS basis functions are used for geometric modeling and analysis:

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

After the control points are rearranged with a unified subscript , the field approximation can then be written as follows:

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

where Î is the parametric coordinate; RI (Î) denotes the NURBS basis function at control point I ; NP is the total number of control points.

On the basis of the principle of minimum potential energy and by substituting the approximated displacements of Eq. (10) into Eq. (7), a system of discretized algebraic equations can be obtained as follows:

KU = f , (11)

where is the vector of the unknown variable at control points; and are the global stiffness matrix and external force vector at control points, respectively.

The element contributions to and are as follows:

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

where D is the Hooke's elasticity matrix.

Reduced Lagrange Multiplier Method For Dirichlet Boundary Conditions

The NURBS basis functions that are often used in the IGA do not satisfy the Kronecker delta function, and the direct imposition of uI (x) = gI (x) on the control points may induce significant errors with deteriorated convergence rates [12]. The novel RLMM is presented in the following.

Reduced lagrange multiplier method

The LMM is easy to understand and implement because it has been widely used for the enforcement of DBCs. However, the LMM increases the problem size and leads to a poorly conditioned matrix equation. To overcome these shortcomings, we modify the standard LMM such that the size of the equation system is reduced. We named the modified method the RLMM. By applying the LMM, the constrained weak form is obtained as follows:

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

The matrix form of the LMM is derived as follows [Ghorashi et al., 2012]:

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

The Lagrange multiplier is approximated by the NURBS basis functions as follows:

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

The degrees of freedom (DOF) u for the control points discretized in the structure can be categorized into two groups: the constrained DOF uB that is the collection of control points ( ) for boundary description, and the unconstrained DOF uA that is the collection of control points ( ) for interior domain description. The NURBS control points can be clearly partitioned into boundary and interior groups because the corresponding NURBS basis functions associated with the interior control points vanish at the boundary. Thus, field approximation can be partitioned as follows:

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

Given that the basis functions of interior control points ( NA) vanish identically at the boundary Γu ,CA =0, Eqs. (22a) and (22b) can then be respectively written as follows:

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

The results systems Eqs. (23a) and (23b) show that the additional unknowns, the Lagrange multiplier are disappeared, the singular problem can be avoid and the size of the equation system is actually reduced to an (n-NB) × (n-NB) matrix problem. Besides, as the NURBS basis functions are used as the interpolating functions of finite element analysis and the Lagrange multiplier (Eq.(17)), so the coefficients Lupinepublishers-openaccess-journals-Civil-engineering-Architechture is similar to a 1D mass matrix along the boundary, and the full quadrature (the number of Gauss points in each direction is p+1) is used to sure that it has an inverse.

Further discussions to the method

Several possibilities are available for the choice of the interpolation space for the Lagrange multiplier, such as (1) the Lagrange shape functions or (2) the same shape functions used in the interpolation of U (i.e., the NURBS basis functions in the IGA). To verify the Babuska-Brezzi stability condition and impose accurate essential boundary conditions, choosing the interpolation of the Lagrange multiplier is important. In the IGA, the matrix form of the LMM equation is expressed as follows:

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

The combination of Eqs. (24a) and (24b) perform a block-wise Gaussian elimination to obtain the pressure Schur complement as follows:

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

The question now is whether CK-1CT is invertible. K-1 is usually invertible for Laplace- and Poisson-type problems; therefore, we only need to consider whether CCT can be invertible. The Babuska-Brezzi stability condition also ensures CCT to be invertible. In the IGA, CA= 0; thus, CCT is reformed to becBcTB . If the NURBS basis functions are used as the interpolating functions of the Lagrange multiplier (Eq.(17)), then Lupinepublishers-openaccess-journals-Civil-engineering-Architechtureis similar to a 1D mass matrix, and CBCTB can be invertible. Therefore, the NURBS basis function is a stable space for the Lagrange multiplier to impose essential boundary conditions in the IGA. The NURBS basis function can obtain an acceptable solution on the numerical examples (Section 4), and the resulting system of equations will not be singular even if the number of Lagrange multipliers is large.

Numerical Examples and Discussion

In this section, three 2D examples are considered to show the accuracy and performance of the RLMM approach. The number of Gauss points in each direction is p+1 in all studies, with p being the order of the NURBS basis function. For convergence analysis, the h-refinement strategy is applied and knots are added to the middle of the knot spans. Moreover, three different measurements of error are calculated to assess the accuracy of the present method. The L2 error norm is expressed as follows:

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

where and are the exact and approximate displacements, respectively.

The H1 error norm is expressed as follows:

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

The strain energy error norm is expressed as follows:

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

where and are the exact and approximate strains, respectively; and are the exact and approximate stresses, respectively.

The exact strain energy norm is defined as follows:

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

Laplace problem on a square domain

The first example considers the Laplace problem. The governing equation and the boundary conditions are expressed as follows:

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

Figure 1: Initial geometry: physical mesh (a) and control mesh (b).

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

As shown in Figure 1, the initial geometry is constructed by a NURBS surface of order p = 3 for both directions with unit weights. The initial parameter space is given by knot vectors κξ = {0,0,0,0,0.5,1,1,1,1} and κη = {0,0 ,0,0,0.5 ,1,1,1,1} . For the convergence analysis, h-refinement is used to refine the mesh by knot insertion. The refined meshes have 16, 64, 256, and 1024 elements (Figure 1).

The convergence results of the RLMM, LMM and the DM are shown in Figs. 2 and 3, respectively. The optimal convergence rates are obtained for the present method and LMM while the results based on the RLMM are identical to those from the LMM. A quadratic convergence rate in the L2 error norm and a cubic convergence rate in the H1 error norm are obtained that are comparable with the quadratic convergence rate in both cases for the direct imposition of boundary conditions. We can also predict that the RLMM can compute faster than traditional LMM because only an (n-NB)x(n- NB) matrix problem (Eq.23) needs to be solved in the former rather than an ill-conditioned (n+NB)x(n+NB) matrix problem [Eq. (16)] in the latter.

Table 1: Comparison of the DOFs and condition number used for the computation among different methods.

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

The comparison of different methods on the number of degrees of freedom and matrix condition number are presented in Table 1. Divided the total DOF n into boundary DOF NB and interior DOF NA, it can be seen from Eq.(23) that the RLMM owns the same number of DOF as the DM, but the boundary DOF NB is imposed directly in DM while it obtained by solving the system of equation with smaller size Eq. (23a) in the RLMM. The total number of unknowns in the LMM will be n+NB as Lagrange multiplier unknowns are equal to boundary DOF NB. On the other hand, the matrix condition number of RLMM and DM are the same as shown in Eq. (23b). As shown in Table 1, the matrix condition number of LMM is three orders of magnitude larger than the RLMM and DM which is even worse for complicated problems. In conclusion, the RLMM has the same matrix condition number as the DM but needs to solve an Eq. (23a) to get the boundary unknowns NB. Compared with Eq. (23b), Eq. (23a) is smaller and would not waste too much time. However, compared with the DM and LMM, additional unknowns and bad matrix condition number are found for the LMM.

Figure 2: Comparison of L2 error norms for Laplace problem on a square domain.

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

The RLMM yields the same accuracy as the LMM while avoiding bad matrix condition number as the DM and also the singularity often met in the LMM for large number of degrees of freedom is no longer encountered (Figures 2 & 3). The absolute error Lupinepublishers-openaccess-journals-Civil-engineering-Architechture distributions utilizing 256 elements are depicted in Fig. 4 to show the accuracy of the proposed approach. The following can be observed from the results: (1) the error domain derived from the RLMM is small, whereas less accuracy at the boundary can be found for the DM; (2) according to the error distributions, although the maximum value of error given by the DM is smaller than that by RLMM, the error region of DM is larger than the those of RLMM (Table 1) (Figure 4).

Figure 3: Comparison of H1 error norms for Laplace problem on a square domain.

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

Figure 4: Distributions of absolute errors obtained by (a) DM, (b) RLMM.

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

A Cantilever beam

A cantilever beam with length L = 48 units and height H = 12 units subjected to a parabolic traction at the free end as shown in Fig. 5 is considered. The material properties of the beam are taken as follows: Young's modulus E = 3.0x107 units, Poisson’s ratio v = 0.3. Plane stress conditions are assumed. The imposed parabolic traction is expressed by the following:

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

where I = H3 /12 denotes the moment of inertia, the shear force applied is P = 1000 units, and the left edge of the beam is clamped. The analytical displacements of this problem are given by the following [35,36]:

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

The order of the NURBS basis functions is p = 2 and the initial parametric space is defined by knot vectors . κξη,={0000.5’U’1} The weights associated with all the basic functions are set to one; thus, the NURBS basis functions are identical to the corresponding B-spline functions. The meshes used for the convergence study are 4 x 4, 8 x 8, 16 x 16, and 32 x 32 elements. The results of the L2 and strain energy error norms shown in Figs. 6 and 7 clearly confirm that the proposed method has smaller solution errors than the DM and preserve the optimal convergence rates (L2 error norm: 3.016, H1 error norm: 2.014) (Figure 5-7).

Figure 5: Description of cantilever beam problem.

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

Figure 6: Comparison of L2 error norms for cantilever beam problem.

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

Figure 7: Comparison of strain energy error norms for cantilever beam problem.

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

Infinite plate with a circular hole

The last example deals with a plane stress infinite plate with a central circular hole under the x-direction traction Tk =1 unit (Figure 8a). Owing to geometrical symmetry, a quarter of the domain is simulated (Figure 8b). To avoid the singularity in the geometric representation, we further consider a finite annular plate in the numerical model (Figure 8c). The exact displacements using Eq. (36) are imposed on the inner and outer edges of the annular plate. The analytical displacement solutions are expressed as follows [35]:

Figure 8: Schematic description of (a) an infinite plate with a circular hole, (b) a quarter of the geometry, and (c) a finite annual plate used for numerical implementation model.

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

The geometry and material parameters are the radius of the circular hole a = 1 unit, the outer radius b = 4 units, Young's modulus E = 105 units, and Poisson's ratio v = 0.3. The initial geometry for this problem is shown in Figure 9. The refined meshes used in the convergence analysis are 4 x 4, 8 x 8, 16 x 16, and 32 x 32 elements. The results of L2 and the strain energy error norms are analyzed and respectively plotted in Figures 10 & 11. The proposed method has less solution errors and achieves optimal convergence rates, thus outperforming DM (Figure 9-11).

Figure 9: Initial geometry of (a) physical mesh and (b) control mesh.

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

Figure 10: Comparison of L2 error norms for infinite plate with a circular hole.

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

Figure 11: Comparison of strain energy error norms for infinite plate with a circular hole.

Lupinepublishers-openaccess-journals-Civil-engineering-Architechture

Conclusion

The Dirichlet boundary conditions (DBCs) in IGA are studied by developing an accurate treatment method, the RLMM, for overcoming the difficulties in the enforcement of DBCs. Compared with the standard LMM, the size of equation system is reduced as the resulting system of equations have made no additional unknowns, the Lagrange multipliers are thus disappeared and the singular problem in LMM is avoided. The benefit of the proposed approach can gain even greater if the number of degrees of freedom of the discretization system is relatively large. The choice of the interpolation space for the Lagrange multiplier is addressed in detail. As a result, the NURBS basis function is found to be a stable space for the Lagrange multiplier for which the DBCs in the IGA is imposed.

Nevertheless, the realization of the improved method in the treatment of DBCs in the IGA is straightforward, and can be integrated into existing IGA codes with less effort. Thus, the improved method can be ideal candidate for practical problems. Based on the numerical investigation of the convergences, the new method yields high accuracy on the solutions and optimal convergence rates, confirming the effectiveness of the enhanced treatment technique of the DBCs in the IGA. The level of accuracy and convergence rates achieved by the RLMM is better than the DM.

Acknowledgment

This work was supported by Natural Science Foundation of Hunan Province of China (2017JJ3306), Foundation of the Education Department of Hunan Province of China (17C1532), and Research Foundation of Xiangtan University (16QDZ15). The financial supports are gratefully acknowledged.

Read More Lupine Publishers Trends in Civil Engineering and its Architecture (TCEIA) Please Click on Below Link: https://lupinepublishers-civilengineering.blogspot.com/





No comments:

Post a Comment