Abstract
The objective of this article is to discuss the existence and the uniqueness of a weighted extended B-spline- (WEB-spline-) based discrete solution for the 2D Navier-Lamé equation of linear elasticity with a different type of mixed boundary condition called boundary condition. Along with the usual weak mixed formulation, we give existence and uniqueness results for weak solution. Then, we illustrate the performance of Ritz–Galerkin schemes for a model problem and applications in linear elasticity. Finally, we discuss several implementation aspects. The numerical tests confirm that, due to the new integration routines, the weighted B-spline solvers have become considerably more efficient.
1. Introduction
The finite element method has become the method of choice for solving many types of partial differential equations in engineering and physical sciences. Important applications include structural mechanics, fluid flow, thermodynamics, and electromagnetic fields [1], using the approximation of Lagrange [2]. This type of approximation has experienced a great restriction in the level of the geometric domain, especially in the case of complicated boundaries such as that in the form of curvilinear graphs. A new way of approximation came a few years ago called B-spline that will solve this problem with excellence.
B-splines have become standard tools in approximation, computer graphics, design, and manufacturing. Recently, they have also been used to construct basis functions for finite element methods. The resulting techniques combine the computational efficiency of regular grids with the geometric flexibility of classical finite elements. Some key advantages are the free choice of order and smoothness, a simple data structure with one parameter per grid point, and the exact representation of boundary conditions. A type of B-spline is called weighted extended B-spline (WEB-spline) [3]. As will be described in the next section, a weight function, , is being multiplied to the B-spline functions to ensure that all WEB-splines satisfy the boundary conditions exactly. The books [4, 5] provide a comprehensive treatment of the relevant theory. Which technique is best suited for a particular problem depends to some extent on the topological form of the simulation domain and its representation. Weighted methods are a good choice for problems with natural (Neumann) boundary conditions or if the part of the boundary, where essential (Dirichlet) boundary conditions are prescribed, has a convenient implicit description. Isogeometric methods can handle domains well which are parametrized over rectangles or cuboids or which can be expressed as union of few such parametrizations, e.g., NURBS representations for CAD/CAM applications. There are also problems where a combination of both techniques might be useful [6].
In this project, we use the implementation of weighted finite element methods to solve the Navier-Lamé system with a different type of mixed boundary condition [7] that generalizes the well-known basis conditions, especially the Dirichlet and the Neumann conditions. This boundary condition helps us to implement a single MATLAB code that summarizes any kind of boundary conditions that we can meet. At the programming level of the method, we can reap several programs in one. The outline of this article is as follows. In the next section, we introduce the modeling of the Navier-Lamé equation and we will give the variational problem that corresponds to this equation. We review the preliminary aspects on WEB-spline-based methods for Navier-Lamé equations. In Section 3, we prove inf-sup conditions for the variational problem in two-dimensional case using WEB-spline-based mesh-free method, which are necessary for the well posedness of the problem.
2. Governing Equation
Linear elasticity is the mathematical study of how each point of the solid object is displaced. The latter produces a deformation while the object becomes subjected to internal stresses due to the prescribed loading conditions, knowing that the linear elasticity models materials in a continuous state. Linear elasticity is a simple case of the nonlinear elasticity theory and is a branch of continuous domain mechanics. The basic theory of linear elasticity is based on the following: the strains (or stress) are infinitesimal or small and the relationships between the components of the stress and the strain are linear, and the linear elasticity is valid only for the states of stress which do not produce a yield. These assumptions cited above are reasonable for the analysis of many engineering materials and in technical design. It is often that this analysis is done using the finite element method.
Let us consider to be a bounded Lipschitz domain with boundary condition which will be presented in a new form that generalizes the Neumann and Dirichlet boundary conditions. Given , , two matricial functions that are invertible or zero. as well as the positive parameters and . When the shear modulus is small, then this kind of problems is called singularly perturbed problems, where the uniform mesh or L2 norm-based error analysis does not converge to the original problem, and hence, one must use adaptive mesh with uniform error analysis. This numerical analysis is given in the papers [8–14].
When solid objects are subjected to external or internal loads, they deform and lead to stress. If the deformation of the solid is relatively small, linear relationships between the components of stress and strain are maintained. Consequently, linear elasticity theory is valid. In practice, linear elasticity theory is applicable to a wide range of natural and engineering materials and thus extensively used in structural analysis and engineering design. The equation of Navier-Lamé below is governed as follows.
A solid object is deformed under the action of forces applied. Any fixed point of the domain and when surface forces are applied to the domain, this point moves to another point . Then, the vector is called displacement. When the movement is small and the solid is elastic, then Hooke’s law gives a relationship between the stress tensor and the strain tensor. is the stress tensor, is the strain tensor, is the identity matrix, and is the shear modulus (or rigidity), where is Lam’s first parameter. The Navier-Lamé equation is given by the law of conservation moment where is the acceleration and is the density of the material. On the other hand,
Then, we have With a simple calculation, we find that
Then, we get
If the solid is in dynamic equilibrium, then we have ; are the external forces applied to the solid. Finally, we find out the equation
We refer the reader to [15, 16] for more information of the elasticity problems.
We create a new unknown that is equal to the divergence of the displacement. The equation of Navier-Lamé becomes
Our mathematical model is the Navier-Lamé system with the boundary condition noted such that is called the Dirichlet matrix and is the Neumann matrix.
There are two strictly positive constants and , such that
is a matrix norm that will be defined below.
If , then is the Neumann boundary; and if , then is the Dirichlet boundary.
We need functional spaces and norms
The variational formulation of the Navier-Lamé problem (6) is as follows.
Find such that
The weak formulation (14) may be restated as follows.
Find
With the bilinear forms
3. Solvability of the Generalized Saddle Point System
In this section, we introduce some existing saddle point theory. Let and be two finite- or infinite-dimensional Hilbert spaces equipped with the inner products and and the induced norms and , respectively. Let , and , be bilinear forms on , and , , respectively, which are bounded; i.e., there are positive constants , , , and such that
Furthermore, we define the following Hilbert spaces:
In addition to assumption (17), we assume that
, for , satisfying the inf-sup condition such that there exists a constant
The bilinear form satisfies the weak coerciveness such as there exists a constant such that
Find where , for and are bounded bilinear forms, where and are bounded linear functionals on and , respectively. Reference [17] contains an existence and uniqueness, for the saddle point problem.
Find such that
Lemma 1. Let be a linear operator, defined by such that for all . Therefore, is continuous and we have .
Proof. is closed in . In fact, let be a sequence of such as on ; then, is cauchy on .
By definition of the coercivity of the bilinear form , we have
We deduce that is cauchy in ; then, there exist such as on .
By the continuity of the operator , we have . The uniqueness of the limit gives , so we obtain that is a closed set of ; it gives that and .
Suppose that ; this implies that there exists nonzero. On the other hand, we have for all (because ); this implies that but ; this contradicts the fact that is nonzero. Finally, .
Theorem 2. With assumptions (17)–(21), the generalized variational problem (22) has a unique solution for any and as long as Further, the following stability estimates hold: where is the solution to (23) and thus has the bounds
Proof. is a closed set included in , so is a Hilbert space.
Problem (22) became
, , and satisfy the conditions of theorem 3.1 of the article in [17]; then, problem (29) has a unique solution (if is reduced to , just take ). We are now looking for such that
For this , the displacement solution of (29) is fixed; we consider the following linear form defined by . This form inherits its continuity from the continuity of and . We define the space such that is a Hilbert space with the norm and inner product inherited from . Problem (30) is equivalent to the problem
According to the Lax Milgram theorem, problem (31) admits a unique solution . Assume that (31) is verified, let us show that (30) is verified.
is closed in , so we have ; then, we obtain for all . Therefore, we have for and ; it implies that .
Conversely, let be the solution of (30).
The linear form is continuous on that is a Hilbert space with the scalar product on .
From the Riesz theorem, there exists a unique such that
This defines an operator which is linear and continuous on such that . From Lemma (1), we have .
If , Riesy implies that there exists such that for all . Since , for this , there exists such that for all
Finally, it has been shown that problems (30) and (31) are equivalent.
Lax Milgram guarantees that problem (31) admits a unique solution . Just take .
Theorem 3. With assumptions (17)–(21), the generalized variational problem (15) has a unique solution for any as long as Further, the following stability estimates hold: where is the solution to (23) and thus has the bounds
Proof. It is sufficient to show that the bilinear forms , , , and have checked the conditions (17)–(21), applying Theorem 2 in the case of .
Remark 4. We rewrite the system (15) in a single equation. For that, let us define the operators , , , and by , , , and so that (15) becomes
The well posedness [17–19] of the above equation is written in the following form. For the variational problem (15), the mapping defines an isomorphism if and only if , , , and satisfy assumptions (17)–(21).
4. WEB-Spline Process
Several types of mesh-less approximations were proposed for various applications. A central problem of the mesh-less Galerkin method is to incorporate boundary conditions of Dirichlet type. The weighted extended B-spline approximation takes care of not only the boundary constraints but also the issue of well conditioning of the Galerkin systems [3, 20]. This essential new feature of constructing well-conditioned basis is of cardinal importance in this work. We are interested in applying the approximation properties of the WEB-spaces for discretizing the Navier-Lamé equations as it have been done for Stokes and Navier-Stokes problems [21, 22].
The standard uniform B-spline of degree is defined by the recursion [5] starting from , the characteristic function of the unit interval between zero and one. Figures 1–3 show the uniform B-splines of degrees one, two, and three. These are also known as linear, quadratic, and cubic B-splines, respectively. In this paper, we use the following notational conventions [5]. For functions and , we write if with a constant which does not depend on the grid width , indices, or arguments of functions. The symbols and . are defined analogously. We describe the procedure of constructing the WEB-splines and discuss the approximating properties of the web-space as a finite element space. For , and define the scaled translates:
is the univariate B-spline of degree of support ; the B-splines are polynomials by pieces on the -grid with vertices and scaled so that the -norm, , is independent of .
Knowing that the tensor product B-spline is the extension of B-spline to higher dimensions.
In general, it can be defined as follows: such that the B-spline is an -variant of degree in the th variable, index and the width of the grid .
With the convention that unless otherwise indicated, for the sake of the problem, the element can be taken as an integer, rather than an integer vector.
For , we partition the grid cells with into interior, boundary, and exterior cells, depending on whether , the interior of intersects , or .
For (the relevant index set for ), if has at least one grid cell completely inside , then is named as inner B-spline; otherwise, it is outer.
The corresponding subsets of are and (c.f Figure 4): .
More details about the construction of WEB-spline basis are given in [5].
Now, it is tempting to use as a finite element approximation space. At first sight, this does not seem feasible since -splines do not conform to the boundary conditions. But this difficulty can be resolved easily by multiplying by a smoothed distance function . Then, the weighted -spline space (WB-space) spanned by weighted -splines is a possible finite element subspace for Dirichlet boundary value problems yielding optimal order approximations. But, the condition number of the Galerkin matrix can become extremely large. This is due to the outer -splines which have only very small support inside . One might think that these basis functions do not contribute much to the approximation power and can simply be omitted. Unfortunately, this is not the case. A suitable solution to the problem of controlling the unstable outer -splines is provided by adjoining them appropriately with the inner -splines. It has been done in such a way that the approximation power of the finite element subspace is retained.
Definition 5 (see [20]). For , WEB-spline is defined as where denotes the center of a grid cell, which is completely inside the domain . The coefficients satisfy , for and are chosen so that all weighted polynomials of order are contained in the WEB-space .
Theorem 6 (see [20]). For an outer index , let be an -dimensional array of inner indices closest to assuming that is small enough so that such an array exists. Then, the coefficients are admissible for constructing WEB-splines according to Definition 5.
5. Boundary Conditions with Resolution Adapted
The imposition of inhomogeneous Dirichlet boundary conditions is essential in numerical analysis of a structure. It is especially difficult and no longer straightforward whenever nonconformal mesh is used to discretize a structure. One of the contributions of this paper is to develop a weighted extended spline basis with high computing accuracy and appropriate with the mixed formulation of boundary conditions that generalize all the cases that can be encountered on the edge of .
With and are two reversible square (in case 2D) matrices, and where the function is a priori known function of spatial coordinates. The inhomogeneous Dirichlet condition is expressed in the case of ; then, we obtain . In general, a boundary condition function may not be known explicitly or is not prescribed globally. More typically, the boundary conditions are prescribed in a piecewise manner, i.e., the boundary conditions are specified by individual functions on each portion of boundary, , , and , . If each Dirichlet boundary is geometrically represented by implicit function so that on [23],
In fact, as the weighting function and boundary value function are expressed in the form of implicit functions, no unique expressions exist but typical forms are given below.
Weighting function acts as a multiplier to modify the original interpolation functions and has to vanish on any Dirichlet boundary . As implicit function is not unique for the definition, weighting function of different forms can be used, as discussed in [24]. Here, two construction techniques are presented. First, a straightforward method is to construct the weighting function by means of the product . Clearly, if any on . Nevertheless, a product of all the implicit functions might lead to a surge of function values which will cause a numerical overflow and then result in poor computing accuracy and robustness. Therefore, each implicit function is preferred to be normalized in advance according to the size of the physical domain. For instance, the implicit function of a circle can be normalized in terms of its radius.
Now, we will choose a basis function web-spline that will be called the web-spline basis function adapted. This basis will respond to any type of nonhomogeneous boundary conditions especially in the case . Let be the displacement solution of problem (15), so we can write with being the solution of (15) in case of homogeneous Dirichlet boundary and is the nonhomogeneous Dirichlet boundary. is written as a linear combination of the web-spline family cited in Definition (5)
It can be rewritten again in the form where denotes the weighting coefficient associated with in which basic properties of are as follows [23].
, with denoting the Kronecker delta function. With this condition, different forms are proposed below for the definition of .
Transfinite interpolation form: the weighting coefficients can be defined by extending the transfinite interpolation that was studied by Rvachev et al. [25] in the CAD community.
The above equation shares a property of symmetry and similarity. is generally supposed to be positive in the physical domain to ensure the nonzero value of the denominator. Evidently, only if and the partition of unity holds with . The value of is used to interpolate normal derivatives prescribed on , which must be 1 greater than the order of the prescribed derivatives.
Finally, the basis that we can propose for the nonhomogeneous boundary condition is as follows:
6. Discretization of the Navier-Lamé Problem Using WEB-Spline Basis
In the Navier-Lamé equation, and are the terms with derivatives of the highest order for the displacement and div-displacement , respectively. Thus, the orders of the differential operators differ by 1. This suggests the rule of thumb. The degree of the basis functions used to approximate the displacement should be one larger than the approximation of the div-displacement. Also, to satisfy the Dirichlet boundary conditions, the displacement basis functions are multiplied by a suitable weight function . Hence, in this article, we choose -linear weighted extended B-spline for displacement approximation along with -Haar wavelet basis function for div-displacement approximation as our linear-constant element and for quadratic-linear element we denote as quadratic WEB-spline and the mean zero linear function (see Figures 1 and 2) as defined above. In the following, the inf-sup condition and, therefore, the well posedness of the discrete Navier-Lamé problem are settled first for the linear-constant element and then for the quadratic-linear element.
6.1. Linear Constant Element
More precisely, we give and , the displacement and div-displacement finite element spaces, respectively, as follows:
A sufficient condition for to satisfy the inf-sup condition is given in [21, 22]. A similar result has been proved in [26] for a different pair .
6.2. Quadratic-Linear Element
For , where is defined as
For constructing the displacement approximation space, we do the following. For , let be the WEB-spline of order 3, which is given by (40) and (48) where is defined as in (51), the tensor product of the scaled translate of the function where the quadrangulation is the collection of all cells such that (those which are fully inside the domain ) and such that (those portions which intersect the boundary ) and is the one dimensional subspace spanned by the function given by
The function is chosen so that the bubble function vanishes on the edges of the cell . is the node corresponding to the cell . The weight function is multiplied to ensure that the bubble function also vanishes on the boundary . The distance function is used near the boundary, where it is free of singularities and blended smoothly with a plateau inside the domain. More precisely, is defined as , , where controls the height of the plateau and the smoothness of the weight function. The plateau facilitates the use of precomputed values while assembling the Galerkin matrix and also avoids the use of high-order quadratures for the integration of the bubble functions supported on the cells which are fully inside the domain.
The displacement approximation space is taken to be .
This pair of discretization spaces satisfies the discrete inf-sup condition.
Lemma 7 (see [3]). A web-basis is a stable basis with respect to the -norm, , where are constants.
We can also bound higher order Sobolev norms in terms of the 2-norm of the coefficients.
Lemma 8 (see [3]). A web-basis satisfies , and is a constant.
Lemma 9 (inverse estimate). Let be the finite element space considered as in (49). Then, there exists a constant such that .
The proof follows from Lemmas 7 and 8.
Lemma 10 (see [21, 22]). Suppose is a family of uniform quadrilateral mesh on the domain and is the -projection operator onto . Then where is a constant which is -independent.
Lemma 11. There exists an operator which satisfies the following properties:
Proof. Let be the usual -projection operator with The first inequality (59) follows from Lemma (10) and the density of in , and the second (62) is an -error estimate. Now we fix a linear mapping such that We may interpret the map as a process with two steps. First, we apply the -projection onto the space of piecewise constant functions. Afterwards, in each cell the constant is replaced by a bubble function with the same integral. In this way, we get More precisely, we define as follows: In view of condition (61), the constant is taken to be . Observe that for and , where are constants, The above inequality is due to Cauchy-Schwarz and the area of the cell . Let us set Verify the two properties in the statement of the lemma. By the virtue of the construction of and we impose that on Because is piecewise linear, we apply Green’s theorem It is easy to prove (57) by applying Green’s theorem and using (56) This proves the assertion (58) with , such as of Lemma (9).
Theorem 12. Let the space be such that . The pair satisfies the discrete inf-sup condition:
Proof. Let be given. By the continuous form of the inf-sup condition (20) and (56) and (57) of Lemma 11, there exist and such that ; then, we find According to (56) and the second condition of (20), we have .
Remark 13. With a similar demonstration as above, the bilinear form (16) on verifies the inf-sup condition cited in Theorem 12.
7. Matrix Form of the Navier-Lamé Problem
After establishing the discrete inf-sup condition, we are ready to prove the existence and uniqueness of the discrete solution. (This part of the section is common to both the linear-constant and quadratic-linear elements, so we call the finite element spaces as and in general.) To find the discrete solution pair , it is enough to find the coefficient vectors of displacement , where , , and div-displacement .
Assuming that , and for . Let us define
The matrix form of the discrete Navier-Lamé equations is written as
By defining the assembled matrices , , , we rewrite the above matrix form as
Decoupling of div-displacement and displacement: from the above matrix form, we write
Premultiplying the first equation above by , followed by , and using the second equation, we get
In the following, we show that is invertible, and hence, the above equation can be solved for .
The positive definiteness follows from the positive
Let us define , with , where of .
The positive definiteness of follows from the positive definiteness of matrix (which is guaranteed by the coercivity condition) and the following inequality which is an equivalent condition to the inf-sup condition (20).
Indeed, one has
Condition (21) gives
The inequality of Cauchy-Schwarz gives
On the other hand,
By taking , , and . Then, we find that
By combining (77), (78), and (81), we obtain as soon as .
This proves the existence of the discrete solution .
8. Numerical Results
We tested the performance of the WEB-spline method for the Navier-Lamé problem (6). We take the centrifugal force and the domain is a rotating steel disc defined by , with , , and the edge of consists of two disjoint parts in the form of two circles such as , .
We give in (42) the Dirichlet matrix on , on , and we give the Neumann matrix on , on , the traction force on .
Figures 5 and 6 show the numerical solution for rotating steel disc with that is fixed at a slightly eccentric axis with radius . We show the domain as well as direction and size (colored) of displacement for each number of cells per coordinate direction.
As mentioned in Table 1, we check for this example the accuracy of the numerical approximations. Since in view of the smoothness of the weighted B-spline basis functions that are continuously differentiable for we can compute the pointwise residuals [27] for the Navier-Lamé boundary value problem (6), and rate . We can substitute the finite element approximation into the partial differential equation Navier-Lamé where for the rotating disc. This would not be possible for standard finite element approximation which usually merely belong to .
We show the relative residual for grid widths . While the residuals decay, estimates of the rates exhibit a regular increasing behavior when is changed from 64 to 256 (e.g., the error pde norm [27] has dropped from to when is changed from 8 to 256) which we attribute to numerical stabilities.
Note that in the case when is quite small, we will refer the readers to see [28–31].
9. Conclusion
Using the equations of linear elasticity as a model problem with the boundary condition , we have described the implementation of finite element methods with weighted B-splines using mixed finite element. In this article, we have shown that our WEB-spline-based quadratic-linear finite elements satisfy the inf-sup condition, which is necessary for the existence and uniqueness of the solution, and we proved the existence of the discrete solution. We have shown the convergence of the numerical solution for the quadratic case. Because of limited regularity, the Navier-Lamé problem will not change by increasing the degree of the WEB-spline.
We have computed the relative errors and rates and have shown that it is of order , which is theoretically correct; these numerical results are attributed to the numerical stabilities of the solution. Moreover, the advantage of this problem with boundary condition is the program level MATLAB; it is enough to make a single program MATLAB and can be reduced to ordinary problems as Dirichlet and Neumann.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.