Numerical solutions of a gradient-elastic Kirchhoﬀ plate model on convex and concave geometries using isogeometric analysis

In this work, we study numerical solutions of a gradient-elastic Kirchhoﬀ plate model on convex and concave geometries. For a convex plate, we ﬁrst show the well-posedness of the model. Then, we split the sixth-order partial diﬀerential equation (PDE) into a system of three second-order PDEs. The solution of the resulting system coincides with that of the original PDE. This is veriﬁed with convergence studies performed by solving the sixth-order PDE directly ( direct method ) using isogeometric analysis (IGA), and the system of second-order PDEs ( split method ) using both IGA and C 0 ﬁnite elements. Next, we study a concave pie-shaped plate which has one re-entrant point. The well-posedness of the model on the concave domain is proved. Numerical solutions obtained using the split method diﬀer signiﬁcantly from that of the direct method. The split method may even lead to nonphysical solutions. We conclude that for gradient-elastic Kirchhoﬀ plates with concave corners, it is necessary to use the direct method with IGA.


Introduction
Classical continuum mechanics has been successfully applied to macro-scale problems in civil [33], mechanical [22], material [20], structural engineering [10] and more [29,30]. Many of these problems exhibit multi-scale features and have broad engineering applications [12,32]. To study such problems, it is necessary to adopt a continuum theory which includes multi-scale effects. However, the standard theory of continuum mechanics does not consider the internal length scale of the material at micro-or nano-scale. Enrichment or modification of the classical theory becomes important to study problems across multiple scales. The theory of strain gradient elasticity [15,34] considers the effects of micro-structures using gradients of strains, and requires the identification of more material parameters. Simplified models are developed to reduce the number of higher-order terms as well as material coefficients [4]. These models can capture the size effects of the material [1,9,46]. In this work, we are interested in numerical solutions of a simplified one-parameter gradient-elastic Kirchhoff plate model developed in [2].
The gradient-elastic Kirchhoff plate model is a sixth-order partial differential equation (PDE) in terms of the transverse deflection. If the sixth-order PDEs is solved directly (direct method ), H 3 -conforming finite element methods (FEM) become necessary. Isogeometric analysis (IGA) [23], which offers basis functions with high-order inter-element continuity, is an ideal candidate for solving third-and higher-order PDEs [18,19,44], including the gradient-elastic Kirchhoff plate model [21,38,42]. There are other numerical methods [5,31,45] based on the direct method. Another approach which is widely used to solve thirdand higher-order PDEs is the split method [3,6,7,13,14,25,43]. By rewriting any thirdor higher-order PDE into a system of first or second-order PDEs, the use of classical C 0 FEM becomes possible. Reformulating the biharmonic equation into two Poisson equations is an example of the split method. In this work, we will use both methods (direct and split) along with IGA and C 0 FEM to solve the gradient-elastic Kirchhoff model.
Past works [5,31,38,42,45] on gradient-elastic Kirchhoff plate models are limited to convex geometries. Geometries with finitely many concave corners are important as they appear in many engineering applications. Can one expect numerical solutions obtained from the direct and split methods on convex and concave geometries to be the same? For the biharmonic equation, solutions obtained using the two methods differ on domains with concave corners [17,37,41,44]. We will borrow the ideas and tools from analyzing the biharmoinc equation to study how the direct and split methods perform on the gradientelastic Kirchhoff plate model on both convex and concave geometries.
In this work, we study numerical solutions of the gradient-elastic Kirchhoff plate model with both convex and concave geometries using direct and split methods. We first show the well-posedness of the model when the plate is convex. Then, we split the sixth-order PDE into a system of three second-order PDEs, which enables the use of classical C 0 FEM. For concave pie-shaped geometries, the existence of a unique solution of the model is also provided. We solve the sixth-order PDE directly using IGA and the system of secondorder PDEs using both IGA and C 0 FEM. Convergence studies of direct and split methods are performed and the observed convergence rates are in good agreement with [38]. Next, numerical examples on a concave pie-shaped domain are presented. Numerical solutions obtained using the split method differ significantly from that of the direct method. The split method leads to nonphysical solutions for certain loads. Thus, for plates with concave corners, it is necessary to use the direct method with IGA. This paper is organized as follows. In Section 2, the gradient-elastic Kirchhoff plate model along with the direct and split formulations is presented. Section 3 discusses discretization methods. Convergence studies on a convex domain and numerical examples on an L-shaped domain are shown in Section 4. Finally, conclusions are made in Section 5.

Gradient-elastic Kirchhoff plate model
In this section, we first present the gradient-elastic Kirchhoff plate model. Then we discuss two methods to obtain a solution.

Model equation
We briefly present the gradient-elastic Kirchhoff plate model, and we refer readers to [39] for a detailed derivation. Let the plate occupy a three-dimensional domain Ω × (−t/2, t/2), where the mid-surface of the plate is denoted as Ω ⊂ R 2 , and 0 < t diam(Ω) is the thickness of the plate. We assume Ω is open and bounded with straight boundaries. In addition, t is assumed to be constant.
Let w : Ω → R be the transverse deflection of the mid-surface. A standard argument of dimension reduction techniques gives the gradient-elastic Kirchhoff plate model equation [38,39] as Here, E and ν are the Young's modulus and Poisson ratio of the material, g is a constant that represents the length scale of the micro-structure of the material and is also called the gradient parameter, and f is the transverse body load scaled by the thickness of the plate. In order to have a well-posed boundary value problem, we need to provide eq. (1) with boundary conditions. Boundary conditions for strain-gradient elastic plates are non-trivial and there is a number of options [4,35]. In this work, we are mainly interested in numerical solutions of the gradient-elastic Kirchhoff plate model on both convex and concave domains using the direct and split methods. For simplicity, we supply eq. (1) with the following boundary conditions and formulate the model as wheref = f /D, w, M g nn , G g nn are given data, and n is the unit vector in outward normal direction to ∂Ω.
1. For a rectangular plate, when g = 0, the gradient-elastic model in eq. (2) is equivalent to the classical simply supported Kirchhoff plate model [39].
2. The boundary value problem in eq. (2) is well-posed, and can be split into a system of equations which has a unique solution.
3. Boundary conditions of eqs. (2c) and (2d) can be imposed naturally in the weak forms of the direct and split methods to be presented later.
By differentiating convex and concave geometries, we discuss next two methods, namely, direct and split methods, to solve eq. (2), and show the well-posedness of the problem on convex geometries and the existence of a solution on a special type of concave geometries.

Direct method on convex geometries
We introduce more notations that are used throughout the paper. The standard L 2 inner product in Ω for scalar-or vector-valued quantities is denoted as (·, ·), and ·, · is the duality pairing on the boundary, ∂Ω. L 2 and H s (s > 0) denote the classical Lebesgue and Sobolev spaces, respectively. Without loss of generality, we further assume w = M g nn = G g nn = 0 for the rest of this section.
It is shown in [39] that eq. (1) equipped with fully clamped boundary conditions is wellposed. We present next the well-posedness of the gradient-elastic Kirchhoff plate model with the boundary conditions given in eq. (2).

Split method on convex geometries
The direct method discussed in the previous section requires w ∈ H 3 , which prohibits the use of classical C 0 FEM. We introduce a split method [3,6,7,13,14,25,41,43] commonly used for third and higher-order PDEs in this section.
From classical results for elliptic boundary value problems [8], one can show that eq. (7) is well-posed by successively solving the system of equations in (10).
, the solutions of eq. (2) and eq. (7) coincide. Using numerical examples, we will verify that numerical solutions of eqs. (2) and (7) converge to the same solution.
Remark 2.5. As shown in eq. (5), the stability constant depends on g. The coercivity constant in eq. (10b) also depends on g. We will explore how the solution behaves with different g using numerical examples.

Concave geometries
For less regular domains such as domains with finitely many concave corners, one can not expect solutions of eqs. (2) and (7) to be the same. Without loss of generality, we consider a pie-shaped domain with only one re-entrant point as shown in fig. 1. In this section, we assume ω ∈ (π, 2π), and replace Ω with Ω ω . It is convenient to express Ω ω in polar coordinates, namely, where the transformation from polar coordinates x = (r, θ) ∈ Ω ω to Cartesian coordinates Because we do not expect solutions of direct method eq. (3) and split method eq. (9) to be the same on Ω ω , we differentiate them as w d and w s , respectively. We next present the uniqueness of w s followed by the existence of w d .

Direct method
On concave domains with re-entrant points, we cannot use the technique adopted in proposition 2.2. By following [17,36,37], we instead show in this section that for a givenf , there exists a unique solution in H 3 (Ω ω ) ∩ H 1 0 (Ω ω ) that satisfies eq. (3). We first introduce a series of functions We note that u ω,j ∈ L 2 (Ω ω ) but u ω,j / ∈ H 1 (Ω ω ). It is easy to check that u w,j satisfies the following, −∆u ω,j = 0, in Ω ω , Next, we show the well-posedness of the gradient-elastic Kirchhoff plate model eq. (2) on the pie-shaped domain Ω ω .
It is easy to see from eq. (19) that w d is different from the solution of the split method w s . The H 1 solution w s fails to capture additional terms as demonstrated in eq. (19).

Discretization
We use the Galerkin method along with the two methods discussed in Sections 2.2.1 and 2.2.2 to solve the gradient-elastic Kirchhoff plate model. In this section, we describe the numerical methods used to discretize the weak forms obtained using the direct and split methods.

Spline basis
We briefly discuss the spline basis used in IGA, and refer to [23] for more details. Setting the weights of all control points to 1, we can reduce NURBS to B-splines. We will use tensor product basis functions. In this paper, we use bases with C p−1 continuous derivatives for B-splines of order p ≥ 1.

IGA for direct method
As discussed in Section 2.2.1, it is necessary to use H 3 -conforming functional spaces for the direct method. The weak form of eq. (2) is given as: find w h ∈ H 3 (Ω), and w h = w on ∂Ω, such that Note that for the rest of the paper, we relax the assumption of the boundary data, namely, M g nn and G g nn do not have to be zero. As discussed in Remark 2.1, eqs. (2c) and (2d) are imposed naturally in eq. (25).
IGA uses spline basis functions (of order p ≥ 1) which facilitates the construction of H p -conforming approximation spaces, and has been used to discretize the gradient-elastic Kirchhoff plate model in [38]. We will use IGA (p ≥ 3) to discretize eq. (25).

IGA and C 0 finite elements for the split method
The split method proposed in Section 2.2.2 only requires the functional spaces to be H 1 -conforming. The weak form of eq. (7) is given by: The above equation confirms that eqs. (2c) and (2d) can be imposed naturally in the weak form. For the split method it is sufficient to use classical C 0 finite elements, which are H 1conforming. For the ease of implementation, we use IGA of order p ≥ 1 to solve eq. (26). We remark that IGA and C 0 finite elements coincide when p = 1.

Numerical examples
In this section, we present numerical examples using the direct and split methods discussed in the previous section to solve the gradient-elastic Kirchhoff plate model. We first perform convergence studies on a convex domain. Then, we study the pie-shaped plate.
Thus, the boundary data become w = 0, M g nn = 0, and G g nn = −2g 2 π 3 [cos(πx) sin(πy), cos(πy) sin(πx)] · n, which are used in eqs. (25) and (26). In this case, Ω is convex andf ∈ L 2 (Ω), thus as in Remark 2.4 we expect the solutions using direct (Section 3.2) and split (Section 3.3) methods coincide. Convergence studies are performed using four uniform meshes. The number of elements in each direction is denoted as N = 2, 4, 8, and 16. On each mesh, we fix g = 0.01 and solve eq. (26) equipped with eqs. (27) and (28) for p = 1, 2, . . . , 5. Numerical solutions are compared against the manufactured solution. Convergence results of w, u 1 , and u 2 in L 2 and H 1 norms are reported in figs. 2 to 4. In addition, eq. (25) together with eqs. (27) and (28) is solved for p = 3, 4, 5. Convergence results of w in L 2 and H 1 norms using the direct method are included in fig. 2 for comparison with the split method. Convergence results of w in H 2 and H 3 norms using the direct method can be found in Appendix A.
As shown in figs. 2a and 2b, convergence rates of w in L 2 and H 1 norms are optimal, namely, for p-th (p ≤ 4) order B-splines, the convergence rates in L 2 and H 1 norms are of order p + 1 and p, respectively. Figures 3 and 4 indicate that optimal convergence rates of the additional variables u 1 and u 2 in L 2 and H 1 norms are also achieved for p ≤ 4. When p = 5, convergence rates of w, u 1 , and u 2 deteriorate in L 2 norm, see figs. 2a, 3a and 4a, but they are optimal in H 1 norm, see figs. 2b, 3b and 4b. As presented in fig. 2a, convergence rates in L 2 norm using the direct method suffer from some deterioration, which has been observed in [38]. Convergence rates in H 1 , H 2 and H 3 are optimal as shown in Appendix A. Therefore, we have verified that on convex domains, numerical solutions for the gradient-elastic Kirchhoff plate model using the direct and split methods converge to the same solution.  Next, we study the convergence behaviour using the split method with linear splines (p = 1) for fixed g = 0.01, 0.05, 0.14, and 0.20. As discussed in Remark 2.5, the stability constants of the direct and split methods depend on g. Thus, convergence rates with gdependence are expected for both methods. It is observed in [38] that convergence rates of w in L 2 and H 1 norms using the direct method change slightly with g. On the contrary, figs. 5a and 5b show that the convergence behaviour of w in L 2 and H 1 norms using the split method are independent of the gradient parameter g. This independence upon g is also true for u 1 and u 2 , and more convergence results for different g are presented in Appendix B.

Concave pie-shaped domain
In this section, we consider the pie-shaped plate with ω = 10π/9 and 3π/2 as shown in fig. 1. We use polar parametrization to represent the pie-shaped geometry [27]. For ω = 10π/9, the control points used to generate the geometry with cubic splines are shown in fig. 6a. The mesh, consisting of 256×802 elements, is demonstrated in fig. 6b. The control points and mesh (256×905 elements) for ω = 3π/2 are similar to fig. 6, and thus are omitted for conciseness. We remark that the corresponding basis functions are H 3 -conforming on Ω ω \O. The inverse of the geometric mapping is singular at O due to repeated control points [11], such as e 1−20 which are shown in fig. 1. Nevertheless, the use of Gaussian quadrature points avoid the evaluation of the basis function at the singular point O. It is shown in [24] that polar parametrization on a sphere yields optimal convergence rates in H 3 norm for the solution of a sixth-order PDE. Thus, we still use the resulting basis functions for the direct method.
(a) Control points (b) Mesh Figure 6: (a) A fraction of control points used to generate the pie-shaped geometry for ω = 10π/9. Points e1, ..., e20 overlap with each other. In the simulation, 256 control points in the tangential direction and 5 control points in the radial direction are used to generate the geometry. (b) Mesh demonstration for ω = 10π/9. In the simulation, the mesh has 256 × 802 elements.
We take g = 0.01 and w = M g nn = G g nn = 0 for all numerical experiments in this section. A nonuniform body load acting on the pie-shaped plate is considered, where d(x) = x 2 + (y − 0.05) 2 , and x = (x, y). The nonuniform load,f in eq. (29) which is centered at F c in fig. 6b, is nonnegative and decays quickly to zero. We remark that the support off is fully resolved by the mesh in fig. 6b. First of all, we solve eq. (25) on Ω ω directly, and denote the solution as w d,h . Then, eq. (26) is solved on the pie-shaped geometry so as to obtain w s . Finally, following the procedure presented section 2.3.2 and utilizing w s , we obtain the H 3 solution w d in eq. (19). The resulting transverse deflection is shown in figs. 7a and 8a for w d,h , in figs. 7b and 8b for w d , and in figs. 7c and 8c for w s . Note that the values of w d that we will report are not exact as we can only evaluate eq. (19) computationally. The computational process involves the approximation of w s , the computation of the coefficients d j in eq. (24), and the computation of GG g u ω,j . This entails the numerical approximation of the diffusion and reaction-diffusion problems with singular source terms. Thus, we do not expect a perfect agreement between w d and w d,h ; but we do expect w d and w d,h to be in good agreement and w s to be fundamentally different. This is indeed what the results will show. Comparing fig. 7a against fig. 7b, and fig. 8a against  fig. 8b, we observe that the contour lines of w d,h and w d are similar everywhere on Ω ω especially near the singular point x = 0. The good agreement of w d,h and w d in contour lines justifies the use of the repeated control points at O. However, the maximum values of transverse displacements, w d,h and w d , are slightly different and show an error of ∼ 16% in fig. 7, and of ∼ 10% in fig. 8.
While w d and w d,h are in good agreement, it is easy to see from figs. 7 and 8 that the H 1 solution w s differs significantly from w d,h and w d . Especially, the contour lines near the concave corner x = 0 of w s is smoother than those of w d,h and w d . In fig. 7, the maximum value of the displacement using the split method shows an error of ∼ 186% with respect to w d and an error of ∼ 233% with respect to w d,h . The difference between w s and w d or w d,h is even more dramatic in fig. 8. The maximum transverse displacement from the split method shows an error of ∼ 340% with respect to w d and an error of ∼ 393% with respect to w d,h . This can be explained by that compared to w d , the H 1 solution w s excludes additional terms induced by the series of functions u ω,j in eq. (15) which are singular near the re-entrant point. Therefore, this shows that H 1 solution of the split method cannot capture important features of the reference solution near the geometric singular point. (c) Split method, w s ∈ H 1 Figure 7: Numerical approximation to w obtained using the direct method in Section 3.2, reference solution in Section 2.3.2, and split method in Section 3.3 with ω = 10π/9, and g = 0.01. Note that the deflection changes sign for the direct method and reference solution, but it is strictly positive for the split method.
In addition, as shown in figs. 7a and 7b, both the IGA direct solution w d,h and the reference solution w d give a sign-changing transverse deflection near the concave corner. Meanwhile, fig. 7c indicates that w s is strictly positive everywhere in Ω ω . The sign-changing effect is more evident in fig. 8 when w = 3/2π. This can be explained using the maximum principle [16]. As stated in Remark 2.7, the elliptic operator, −∆w + cw, c ≥ 0, has maximum principle. Therefore, a positive loadf , such as eq. (29), leads to positive u 2 in eq. (13). As a result, u 1 in eq. (13) is positive. Then, w s obtained using eq. (14) becomes positive. However, the sign of w d ∈ H 3 in eq. (19) depends on bothf and the geometric singularity at the concave corners [37], namely, the angle ω.
These results show the superiority of the direct method with IGA over the split method. The reason is twofold. First, the contour lines of the solution obtained using direct method is similar to that of the reference solution; while the solution from the split method is very different from the reference solution. Second, for a positive body load, the reference solution and the solution from the direct method might change signs near the concave corner; while the split method result in positive solution.

Conclusion
In this work, we have studied numerical solutions of a gradient-elastic Kirchhoff plate model. First, we assumed that the plate geometry is convex. We showed the well-posedness of the model in H 3 . Then, we split the sixth-order PDE into a system of three secondorder PDEs, and the resulting solution is only in H 1 . The well-posedness of the system of equations is trivial, and the solutions obtained using the direct and split methods coincide. Next, we consider a pie-shaped domain with one re-entrant point. For any given data in L 2 , the uniqueness of an H 3 solution of the model is proved by utilizing the H 1 solution which is obtained by the system of second-order PDEs. The difference between the H 3 solution and H 1 solution, is presented.
The fact that the H 3 and H 1 solutions coincide on convex geometry is verified numerically by conducting convergence studies using a manufactured solution. Optimal convergence rates are obtained using the direct method along with IGA and the split method with IGA and C 0 FEM. In addition, for the range of the gradient parameters studied in this work, we observe an independence between the performance of the split method and the gradient parameter. However, it is found in [38] that convergence rates in L 2 and H 1 norms using the direct method depend moderately on the gradient parameter.
Next, numerical examples on a concave pie-shaped domain, which is subjected to a nonuniform body load, were studied. The transverse deflection obtained using the direct method agrees well with the reference solution. The H 1 solution, which is obtained using the split method, differs significantly from the reference solution in H 3 . The deflection contour lines near the concave corner are also dissimilar. More importantly, for a positive nonuniform body load, the reference solution and the solution from the direct method both change sign near the concave corner while the H 1 solution from the split method preserves the positivity. For plates with geometric singularities, it is hence necessary to use the direct method with IGA.
This work opens up a number of possibilities for future work. We have only considered a boundary condition that converges to the classical simply supported plate model as the gradient parameter, g, goes to zero. It would be interesting to include more types of boundary conditions such as, singly and doubly clamped, double simply supported, and free boundaries. The boundary of the plate is limited to straight boundaries. Extending the current work to curved boundaries would significantly increase our understanding of the model and the algorithms. For the pie-shaped geometry, even though the reference solution and the solution from the direct method agree well, they are not the same. The difference measured in H 3 norm becomes even more evident. The discrepancy might be caused by the singular parametrization of the geometry and by the singularity of the basis function at the singular point. Thus, it is necessary to develop basis functions that are H 3 -conforming everywhere for such concave geometry.