R-Function and variation method for bending problem of clamped thin plate with complex shape

Solving ordinary thin plate bending problem in engineering, only a few analytical solutions with simple boundary shapes have been proposed. When using numerical methods (e.g. the variational method) to solve the problem, the trial functions can be found only it exhibits a simple boundary shape. The R-functions can be applied to solve the problem with complex boundary shapes. In the paper, the R-function theory is combined with the variational method to study the thin plate bending problem with the complex boundary shape. The paper employs the R-function theory to express the complex area as the implicit function, so it is easily to build the trial function of the complex shape thin plate, which satisfies with the complex boundary conditions. The variational principle and the R-function theory are introduced, and the variational equation of thin plate bending problem is derived. The feasibility and correctness of this method are verified by five numerical examples of rectangular, I-shaped, T-shaped, U-shaped, and L-shaped thin plates, and the results of this method are compared with that of other literatures and ANSYS finite element method (FEM). The results of the method show a good agreement with the calculation results of literatures and FEM.


Introduction
Thin plate refers to a flat plate with the thickness significantly smaller than the size of the plate surface. The problem of thin plate bending is recognized as a critical problem in the elastic theory. On the whole, the car body floor, the engine cylinder, the gearbox housing, the floor of building structure, the bridge deck, the aircraft panel, and the missile panel in aerospace engineering belong to thin plate bending problem. 1 Given Kirchhoff assumption, the governing equation of elastic thin plate bending problem is the biharmonic equation. To work out this type of higher order partial differential equation boundary value problem, a few regular plates with simple geometric shapes and boundary conditions can be obtained analytical solutions. So most engineering practical problems should be solved by using the finite element method (FEM), the boundary element method and other numerical methods. 2 Two main methods are available for solving the mentioned problems, namely, the numerical method and the analytical method. Analytical methods include the superposition method, the quasi-Green function method, the trigonometric series method, the hierarchical wavelet method, the extended separationof-variable method, the optimized separation-ofvariable method, the finite integral transform method and the up-to-date symplectic superposition method. Numerical methods cover the boundary element method, the R-function method, the difference method, the finite element method, the weighted residual method, the high-order meshless method, the spline finite point method, the three-hinge-line method, and the wave-based method. In a wide range of methods, the deflection functions should be selected in advance, and the selection of the mentioned functions is arbitrary without any definite rule. In several novel analysis methods, the bending function is not required to be selected in advance (e.g. the symplectic superposition method and the finite integral transformation). Zhang et al. proposed a two-dimensional generalized finite integral transform method for new analytic bending solutions of orthotropic rectangular thin foundation plates. 3 A Hamiltonian system-based variational principle via the Lagrangian multiplier method was proposed to formulate the thin plate buckling in the symplectic space. 4 Hu et al. presented a first attempt to extend an up-to-date symplectic superposition method to linear buckling of side-cracked rectangular thin plates. 5 An et al. developed a double finite integral transform method for analytical bending solutions of non-Levytype cylindrical shell panels without a free edge that were not obtained by classical semi-inverse methods. 6 Numerous engineering calculation problems (e.g. arbitrary geometry, complex boundary conditions, diverse, and uneven material properties) can be obtained directly from the mathematical model by numerical methods. 7 The bending, vibration, and buckling of orthotropic rectangular thin plate under various boundary conditions by employing Fourier series and superposition principle with MATLAB programming were studied by Li 8 . Li et al. adopted the quasi-Green's function method to study vibration analysis of simply supported polygonal thin plates on Winkler foundation. 9-12 Huang obtained the general solution of the differential equation of the bending problem of anisotropic rectangular thin plates by exploiting the trigonometric series solution as complex numbers and the double sinusoidal series solution, and solved the bending problems of various boundaries under arbitrary loads. 13 Yu applied a hierarchical wavelet homotopy methodology for studying the nonlinear bending of anisotropic thin plates. 14 Xing and Wang 15 adopted a method to solve closed-form analytical solutions of the free vibration problems of orthotropic rectangular thin plates with arbitrary homogeneous boundary conditions. The optimized separation-of-variable method for the free vibration of orthotropic rectangular thin plates was investigated by Xing and Wang. 16 The finite integral transform method was developed to conduct the bending analysis of thin plates with the combination of simply supported, clamped, and free boundary conditions. 17 Hu et al. 18 developed an advanced symplectic superposition method for several novel analytic free vibration solutions of side-cracked rectangular thin plates. Feng 19 employed the boundary element method to solve the bending problem of thin plate with small deflection and using variational principle, and derived the boundary integral equation of the bending vibration of the base plate and thin plate on Winkler elastic foundation. Dong 20 used the boundary element method to study the bending problem of elastic thin plate in inhomogeneous infinite domain. The boundary element method was adopted to reduce the dimension of the problem and simplify the calculation. 21 Chen and Xiao 22 studied the numerical characteristics of the preconditioned GMRES algorithm in solving the boundary element method equations of elastic mechanics problems. 22 Kurpa and Mazur 23 proposed a method of R-function theory for studying parametric vibration of complex shape plates with irregular shapes and different boundary conditions. Kurpa et al. applied the R-function theory for the free vibration of shallow shells, combined the spline approximation method, and the Ritz method to analyze the stress-strain state of the laminated shallow shell and the vibration of the laminated shallow shell with complex notch. [24][25][26] The variational-difference method was used to study thin plate bending problem. 27 The classical theory of the thin plate bending problem according to the generalized variational principle of three types of variables in 3D elastic theory was analyzed by Zhang et al. 28 The bending problems of Kirchhoff plate and Winkler plate under different transverse loads effectively was researched by generalized difference method. 29 The approximate calculation method for rectangular thin plate by using the principle of minimum potential energy and calculating the deflection of thin plate under uniform load was investigated by Dong and Lv. 30 The bending vibration problems of rectangular thin plates under simply supported boundary and fixed boundary conditions was studied by conducting the finite element analysis. 31 A method for solving the bending problems of thin plates with arbitrary geometric boundaries and supporting conditions was studied by combining the weighted residual method. 32 The application of small deflection theory of elastic thin plate in engineering design was conducted by Xu et al. 33 Analytical solutions for bending, buckling, and vibration analyses of thick rectangular plates with various boundary conditions were presented using two variable refined plate theory. 34 Based on the parametrized mixed variational principle, a high-performance finite element model has been developed for bending and vibration analysis of thick plates. 35 A novel mixed-field theory with relatively low number of unknown variables was introduced for bending and vibration analysis of multi-layered composite plates. 36 A four-node quadrilateral partial mixed plate element with low degrees of freedom was developed to conduct the static and free vibration analysis of functionally graded material plates rested on Winkler-Pasternak elastic foundations. 37 A spline finite point method was developed to study the nonlinear bending behavior of functionally graded material plates with different thickness in thermal environments. 38 Wang et al. presented a mesh free Galerkin method to study SG thin beams and plates to satisfy the continuity and convergence requirement. Besides, moving least square or reproducing kernel shape functions were employed with the cubic approximation bases. 39 The three-hinge-line method was adopted to study the crushing deformation of sheet metals. 40 Liu et al. 41 applied finite element method and wave-based method simulations to study data on the flexural vibration of thin plates.
The Ritz method which is a common energy method has been extensively used in mechanical research because of its simple solving process. 42 The R-function theory can be used to describe complex regions as implicit functions. The key is to express the boundary of the problem by standardized equation v = 0, and the region of the problem is expressed by inequality v.0. 23 The introduction of the R-function theory can easily build the deflection function w of the irregular thin plate to meet the boundary conditions, and then the specific deflection function of the concave-convex irregular thin plate can be obtained by combining the variational method.
It is required to develop efficient and accurate numerical analysis methods for mechanical problems (e.g. bending, vibration, and buckling problem of thin plate) to ensure the design and safety of plate structures in practical engineering. 43 Besides, the governing equation of thin plate bending problem is a fourth order partial differential equation with deflection as independent variable. Developing numerical methods for the mentioned problem is of great scientific significance for solving high-order partial differential equations. 44 The paper introduces the calculation principle and gives numerical examples of rectangular, I-shaped, Lshaped, T-shaped, and U-shaped thin plates. After deducing the variational equation, the results are calculated by MATLAB program, and the feasibility of the method is verified by comparing with the results of other literatures and ANSYS finite element method (FEM) calculation results which show a good agreement.

Calculation of small deflection bending of thin plate by variational method
To the small deflection bending problem of thin plates, according to the calculation assumption, the deformation components e z , g zx , g zy are not considered, thus the expression of deformation potential energy is simplified to the expression as follows The stress component and strain component in the mentioned equation are expressed by deflection w as follows Substituting equations (2)-(7) into equation (1), it yields For equation (8) each item does not change with z, integral z, from À d = 2 to d = 2 , and using the equation The flexural stiffness D in the equation, in the case of a plate with variable thickness, will be a function of x and y, since d will be a function of x and y.
In the thin plate of equal thickness, D is a constant, and equation (8) can be rewritten as where According to Green's theorem, there are ðð According the equation (12), it yields The line integral on the right is along the boundary of the thin plate.
All the boundaries of the thin plate are fixed edges. Thus, regardless of the boundary shape, there is a firstorder partial derivative on the boundary equal to zero.
The boundary conditions of the fixed edge are According to the calculation assumption and geometric equation in the small deflection bending problem of thin plate, the displacement component u and v can be expressed by deflection w, which need not be taken as the basic unknown function. Thus, there's only one w fundamentally unknown function. Now set the w expression as follows Where, C m are m individual undetermined coefficients independent of each other, w m is the setting function satisfying the displacement boundary condition of thin plate. In this way, regardless of the value of C m , the deflection w shown in equation (15) can always satisfy the displacement boundary condition. Note that the variation of deflection w is only realized by the variation of coefficient C m . As for the set function w, it changes only with coordinates and is completely independent of the above variation.
To determine the coefficients C m , the expression is adopted as follows To the bending problem of thin plates, both the body force and the surface force are attributed to the load q, and there is dS = dxdy on the panel, so equation (11) is written as The C m can be obtained from the m individual linear equations, thus the deflection w can be obtained by equation (15), and the internal force of the thin plate can be calculated.

R-function theory
According to R-function theory, [23][24][25][26] the v 0 describing a complex region can be described by Boolean operations _ a and^a of each subdomain, which are corresponding to the union set and intersection set. The function v l describing each subdomain satisfies the first-order normalization equation, and v 0 also satisfies the first-order normalization equation.
Definition of the first-order normalized equation as follows Let X and Y satisfy the definition of the first-order normalized equation, then the following Boolean operations X _ a Y and X^aY are also the first-order normalized equations Where À 1 ł a ł 1.
For the thin plate bending problem, the R-function is introduced into equation (15). This paper discusses the situation of the clamped edge, according to the properties of the R-function, the obtained v 0 is squared to satisfy the boundary conditions of the clamped edge.

Derivation of variation equations
The expression of deflection is Thus, the boundary conditions are expressed by v 2 0 . The expression of deformation potential energy of thin plate with clamped boundary conditions as follows Namely Substituting equation (18) into equation (19), we can get Take the derivative of C i (i = 1, 2, 3 Á Á Á Á Á Á m À 1, m), we can obtain Substitute the R-function expression into equation (17), the follow expression can be obtained Take the derivative of C i (i = 1, 2, 3 Á Á Á Á Á Á m À 1, m), we can get Substituting equations (22)-(25) into equation (27), we can get The coefficient equation system can be obtained by derivation as follow Where A = a ij À Á m 3 m , Xia and li Here i = 1, 2, Á Á Á Á Á Á , m; j = 1, 2, Á Á Á Á Á Á m, O is the area contained in the plane of thin plate (x-y) in the Cartesian coordinate system, which C m can be obtained from equation (17), and C m can be substituted into the equation (20) to obtain the deflection expression.

Example 1 Rectangular thin plate
Rectangular thin plates shown in Figure 1 with side lengths of 2a and 2b, where a = b = 0:5, the elastic modulus E of the material is 100GPa, the Poisson's ratio n is 0.3, the thickness d of the thin plate is 0:05m, and the uniform load is q = 10 6 Pa.
According to the R-function theory, 23-26 v 0 can be taken as follows Setting the deflection w to be an even function of x and y, so w m take 1, x 2 , y 2 , x 4 , y 4 , x 2 y 2 , x 6 , y 6 , x 2 y 4 , and x 4 y 2 .
Combining R-function with the trial function of variational method, in the process of calculating integral, it will appear that the denominator of integral function is zero on the boundary of integral. In this case, it is impossible to calculate the integral of integral function. Thus, the rectangular area is divided to n 3 n numbers by gridding.
When number of the trial function is m = 10, deflection of the center point of thin plate is listed in Table 1 Table 2, and m takes different values. It can be seen from Tables 1 and 2 that the deflection of thin plate calculated by the method in this paper is convergent and effective, and it is a good agreement with the analytic solution 45 and the FEM solution. These prove the feasibility and correctness of this method. When the number of grids is 40 3 40, m take different values, the deflection curve of the center line y = 0 of the plate is shown in Figure 2 which show the convergence of the method.
In FEM modeling, SHELL181 is selected as the element type of the thin plate, and the edge length of each element of the rectangular thin plate is 0.05 m. The rectangular thin plate is divided into a total of 400 square   grids. When m = 10, the results of the paper and FEM are presented in Figure 3 for comparison, and which shows a good agreement.

Example 2 I-shaped thin plate
A I-shaped thin plate is shown in Figure 4. The elastic modulus E of the material is 100Gpa, the Poisson's ratio n is 0.3, the thickness d of the thin plate is 0:05m, and the uniform load is q = 5 3 10 5 pa. We set a = b = 1, c = d = 0:5.
According to the R-function theory, [23][24][25][26] Where v 1 = a 2 Àx 2 2a ø 0, v 2 = b 2 Ày 2 2b ø 0, v 3 = c 2 Àx 2 2c ø 0 and v 4 = y 2 Àd 2 2d .0. w m take 1, x 2 , y 2 , x 4 , y 4 , x 2 y 2 , x 6 , y 6 , x 2 y 4 , and x 4 y 2 . When programming with MATLAB, thin plate is divided into three small rectangles, the first rectangular part is divided into 4n 3 n, the second rectangle is divided into 2n 3 2n, and the third rectangle is divided into 4n 3 n. When m = 10, and n takes different values, the results are listed in Table 3. When n = 9, and m takes different values, the results of center point are listed in Table 4, and the deflection curve of the center line y = 0 of the thin plate is calculated which shown in Figure 5 which show the convergence of the method. In FEM model, SHELL181 is taken as the element type of the thin plate, and the respective element of the Ishaped thin plate takes the edge length of 0.05 m. The I-shaped thin plate is divided into a total of 1200 square grids. When n = 9 and m = 10, the results are    compared with the FEM results which shows a good agreement in Figure 6, and that proves the convergence and correctness of the method in this paper.

Example 3 T-shaped thin plate
A T-shaped thin plate is shown in Figure 7. The elastic modulus E of the material is 100Gpa, the Poisson's ratio n is 0.3, the thickness d of the thin plate is 0:005m, and the uniform load q = 10 3 , and set a = b = 0:5, c = d = 0:25.
According to the R-function theory. [23][24][25][26] v 0 can be taken Where v 1 = a 2 Àx 2 2a ø 0,v 2 = b 2 Ày 2 2b ø0,v 3 = c 2 Àx 2 2c ø 0 and v 4 = d À y ð Þ.0. w m take 1, y, x 2 , y 2 , y 3 , x 2 y, x 4 , y 4 , x 2 y 2 , and y 5 . When programming with MATLAB, the T-shaped thin plate is divided into three rectangles. The first rectangular part is divided into 2n 3 n grids, the second rectangular part is divided into 4n 3 2n, and the third rectangular part is divided into 4n 3 n. When m = 10,n takes different values, the results are listed in Table 5. When n = 9, m takes different values, the results of center point are listed in Table 6, and the deflection curve of center line y = 0 of thin plate is calculated as shown in Figure 8 which show the convergence of the method. In FEM model, SHELL181 is taken as the element type of the thin plate, and the element of the T-shaped thin plate is the edge length of 0.05 m. The T-shaped thin plate is divided into a total of 350 rectangular grids. When n = 9 and m = 10, the results of the paper and FEM are compared which shows a good agreement in Figure 9, and that proves the convergence and correctness of this method.

Example 4 U-shaped thin plate
A U-shaped thin plate is shown in Figure 10. The elastic modulus E of the material is 100Gpa, the Poisson's ratio n is 0.3, the thickness d of the thin plate is 0:005m, and the uniform load is q = 10 3 , and set a = b = 1, c = d = 0:5.
According to the R-function theory, [23][24][25][26] Þø 0 and v 4 = y 2 Àd 2 2d ø 0. w m take 1, y, x 2 , y 2 , y 3 , x 2 y, x 4 , y 4 , x 2 y 2 , and y 5 . When programming with MATLAB, the U-shaped thin plate is divided into three rectangles. The first rectangular part is divided into 4n 3 n grids, the second rectangular part is divided into 3n 3 2n, and the third rectangular part is divided into 4n 3 n. When m = 10, n takes different values, the results of center point are listed in Table 7. When n = 9, m takes different values, the results of center point are listed in Table 8, and the deflection curve of center line y = 0 of thin plate is calculated as shown in Figure 11 which show the convergence of the method. In FEM model, SHELL181 is taken as the element type of the thin plate, and the element of the U-shaped thin plate is the edge length of 0.05 m. The U-shaped thin plate is divided into a total of 1400 rectangular grids. When n = 9 and m = 10, the results of the paper and FEM are compared which shows a good agreement in Figure 12, and that proves the convergence and correctness of this method.

Example 5 L-shaped thin plate
A L-shaped thin plate is shown in Figure 13. The elastic modulus E of the material is 100Gpa, the Poisson's ratio    Figure 11. The curve of w at y = 0.
n is 0.3, the thickness d of the thin plate is 0:005m, and the uniform load is q = 5 3 10 2 pa, and set a = b = 0:5, c = d = 0:125. According to the R-function theory, [23][24][25][26] w m take 1, x, y, x 2 , xy, y 2 , x 3 , x 2 y, xy 2 , and y 3 . When programming with MATLAB, the L-shaped thin plate is divided into two rectangles. The first rectangular part is divided into 7n 3 n grids, the second rectangular part is divided into 8n 3 7n. When m = 10,n takes different values, the results of center point are listed in Table 9. When n = 9, m takes different values, the results of center point are listed in Table 10, and the deflection curve of center line y = 0 of thin plate is calculated as shown in Figure 14 which show the convergence of the method. In FEM model, SHELL181 is taken as the element type of the thin plate, and the element of the L-shaped thin plate is the edge length of 0.025 m. The L-shaped thin plate is divided into a total of 1575 rectangular grids. When n = 9 and m = 10, the results of the paper and FEM are compared which shows a good agreement in Figure 15, and that proves the convergence and correctness of this method.
In the calculation of example 2, example 3, example 4, and example 5, the final results have some errors with the results calculated by the FEM, which may be caused by three factors as follows.
Firstly, the trial function selected in this paper may meet the additional non-existence boundary conditions    in the non-boundary area and impact the calculation results of this research. For instance, the trial function x 2 will constantly maintain zero on the line x = 0, which reveals that the trial function will maintain zero deflection on the line x = 0, whereas this obviously does not comply with the conditions under the consideration of this research. Secondly, in this paper, the number of the trial functions may be used in the process of research is insufficient, the data of example 2 can indicate that when m is small, the results will cause a large error compared with the results of FEM. Thus, to use the variational method for obtaining the correct results, a maximal number of trial functions should be taken.
Thirdly, when the R-function is integral, the denominator is equated to zero on the boundary, which reveals that the R-function cannot be solved by the traditional symbolic integral method. It must be solved by the method of adding the grid after meshing consistent with the FEM, and this also means that the density of the grid will affect the correctness of the results. Accordingly, the R-function and the variational method should be employed to solve the problem.

Conclusions
The variational method, as an energy method, is capable of effectively calculating the deflection variation of the thin plate bending problem. The disadvantage of the variational method is that when the shape of the thin plate is very complex, the advantage of the R-function should be used. The advantage of the R-function is that the implicit function can be used to represent the complex boundary shape. Taking the advantages of the R-function, the combination of the R-function and the function of the variational method can be applied easily to calculate the deflection variation of the thin plate bending problem with complex boundary conditions.
In this paper, the effect of the combination of the R-function and the variational method is illustrated by the theoretical derivation and relevant examples. The advantage of R-function is that it can simply represent the complex boundary shape. For the thin plate often used in engineering, the advantages of R-function are reflected. As one of the energy methods, the variational method can easily calculate the deflection variation of the thin plate. Combining the variational method with the R-function theory to solve the problem with complex boundary conditions, the thin plate bending problem with complex shape can be easily calculated.
In this paper, the examples are presented to illustrate that the combination of the R-function and the variational method can be adopted to calculate the bending problem of thin plates with complex shapes. It is suggested that the combination of the R-function and the variational method can feasibly calculate the bending problem of thin plates with complex shapes. Moreover, the correlation between the results and the trial function of variational method and the number of integral grids is listed, thereby it carry forward the R-function on solving the bending problem of thin plates.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors gratefully acknowledge the financial support provided by the Science and Technology Scheme of Guangzhou City (no.201904010141).

Data availability statement
The data used to support the findings of this study are available from the corresponding author upon request.