Abstract
A novel numerical manifold method was derived from the cubic B-spline basis function. The new interpolation function is characterized by high-order coordination at the boundary of a manifold element. The linear elastic-dynamic equation used to solve the bending vibration of thin plates was derived according to the principle of minimum instantaneous potential energy. The method for the initialization of the dynamic equation and its solution process were provided. Moreover, the analysis showed that the calculated stiffness matrix exhibited favorable performance. Numerical results showed that the generalized degrees of freedom were significantly fewer and that the calculation accuracy was higher for the manifold method than for the conventional finite element method.
1. Introduction
According to the concepts of topological manifold and differentiable manifold and based on the discontinuous dynamics of DDA block systems, the numerical manifold method (NMM) adopts the finite cover technique to set up a unified calculation form for the finite element method, the DDA, and the analytical method [1]. NMM combines the advantages of the finite element method and DDA and can deal with various types of complex geometric and material boundaries. Thus, NMM has an obvious advantage in the simulation of large geometric deformations and joint deformations, as well as dynamic and static crossings [2]. The application research on NMM mainly focuses on certain fields, such as geotechnical engineering [3, 4] and crack propagation [5–8]. However, NMM has recently been successfully applied to fatigue failure [9], fluid calculations [10], and seepage [11]. In theoretical research, most NMMs adopt shape functions in the finite element method as the weight functions for the manifold method and use the local cover functions of high-order polynomials to improve the continuity and interpolation accuracy of elements [12–14]. Notably, a polynomial form based on the global coordinate system will increase the calculated condition number of the stiffness matrix or even cause a solution failure. For this reason, Lin et al. [15] improved the form of local cover functions to enhance the computational efficiency and interpolation accuracy of elements. However, the number of generalized degrees of freedom remains large. The mesh-less numerical manifold method is more flexible in the selection of mathematical meshes and can avoid problems resulting from mesh generation [16]. However, plenty of matching points have to be used in the calculation, and matrix inversions have to be intensively performed. Thus, the number of calculations remains large, but the integral accuracy is low, which restricts the application [17, 18]. Recently, Wen and Luo [19] studied the numerical manifold method for polygonal elements, which has a stronger adaptability to physical meshes. The newly developed numerical manifold method based on quadratic B-spline interpolation can further improve the continuity of elements and prevent the linear correlation problem in the conventional numerical manifold method [20, 21].
The discrete finite element solution is the most widely adopted solution for structural dynamics. The interpolation function required by the calculation is based on the static node structure, which is inconsistent with actual problems. The interpolation shape function is related to the element shape and the coordinate system, such that the mesh requirements are relatively high, and the interpolation accuracy is obviously reduced as the structure changes significantly. For the finite elements with mesh structures, the shape functions usually have a low order, which indicates that they cannot meet the requirements of the continuity of displacement by structural deformation. NMM based on finite cover can solve the above problems effectively. The local cover function corresponding to the node displacement of the finite element is no longer relevant to the nodes. If the order of the local cover function is increased, the continuity and interpolation accuracy of the elements can easily be improved. In particular, the manifold method can adopt fixed mathematical meshes that are independent of the physical meshes (material boundary) in setting up the unified interpolation scheme. This method has the advantages of the Lagrangian and Eulerian methods. The manifold method can also avoid problems such as mesh distortion, moving boundary, and migration terms that are difficult to deal with in the finite element and conventional numerical manifold methods [22]. Based on the concept of fixed meshes and considering the deficiency of the existing manifold method in terms of interpolation accuracy, continuity, and solution efficiency, a manifold element with a new cover form and interpolation is constructed. To apply the new element to the solution of structural dynamics, the dynamic theory based on the numerical manifold method is investigated. The numerical manifold scheme for the bending vibration of thin plates is derived, and the final calculation results are comparatively analyzed.
2. Basic Principles of NMM
2.1. Basic Concept of NMM
NMM consists of block dynamics, simplex integration, and a finite cover system. Figure 1 shows the concrete structure [23]. NMM adopts two sets of independent meshes: a mathematical mesh (cover) and a physical mesh (cover). The mathematical mesh is used to define the interpolation accuracy of the solution based on user-defined forms, and the selected mathematical covers can overlap with one another. Consistency with the physical cover is unnecessary, but the mathematical cover has to include the whole material under study. The physical cover is used to describe the geometric properties of the material and the domain of integration, such as the material boundaries, cracks, blocks, and the interfaces of different materials.
In the finite cover system, the intersection of the mathematical cover and the physical mesh is defined as a manifold element; that is, the mathematical cover is subdivided by a physical mesh. The weight function in the manifold method is equivalent to a shape function in the finite element method. If the proper weight function is selected, the weighted summation of the local cover functions can generate the overall cover function of the whole solution domain.
2.2. Generation of Overall Cover Function
For simplicity, the manifold method adopts a mathematical cover with a regular structure or a finite element mesh to generate the desired weight function. The generated weight function should meet the property of element decomposition. Meanwhile, the local cover function corresponding to each mathematical cover is not directly relevant to the nodes. This function can be a constant, linear, and high-order polynomial or take other local series forms. The two-dimensional problem can be taken as an example. The local cover function of any mathematical cover is defined as follows:
The local cover function can be expressed as the linear combination of the basis function with a given -order and the corresponding unknown coefficients and . The matrix form is as follows:
Any manifold element is formed by the intersection of the cover ) subdivided by a physical cover. The weight function corresponding to each cover is used to solve the overall cover function in the domain, which can be represented as follows:
3. Two-Dimensional Manifold Method for Cubic B-Spline Interpolation
3.1. Definition and Property of the Cubic B-Spline Basis Function
The B-spline basis function can be defined in many ways to facilitate computer programming. In this paper, the reverse derivation method is used to define it [24]. is assumed to be a monotonically nondecreasing sequence; that is, and , where is the B-spline node. is the th B-spline curve with the following definition: where is the node interval of th B-spline. During the process of reverse derivation, is required.
B-spline has a local support feature. In other words, when , , when , . All B-spline curves in the B-spline node interval meet the property of the element decomposition [19]; that is,
According to the definition of the B-spline basis function, the mathematical expression of the cubic B-spline curve expressed in piecewise can be derived as
Figure 2 shows the cubic B-spline curves when the B-spline nodes are evenly distributed. Formula (6) can be used to solve all the B-spline curves in the figure. Figure 3 shows the B-spline curves within the B-spline interval or within the element. The mathematical expression is as follows:
3.2. Two-Dimensional Manifold Element of the Cubic B-Spline
A method similar to bilinear interpolation is adopted to construct the B-spline overall interpolation function (as shown in Formula (8)) within the rectangular element domain of and :
In the formula, is the generalized displacement vector of the manifold element , such that
Each term of is , such that
In the same way, cubic B-spline curves of , , , and within can be obtained in piecewise form. The calculated B-spline functions are substituted into Formula (8) to obtain the final overall cover function of the manifold element.
We suppose that the numbers of manifold elements in the directions of , are, respectively, , , and the corresponding element nodes are . Through Formulas (7), a certain number of virtual nodes have to be added to construct the B-spline basis functions. We suppose that the added virtual nodes are , , , ; , , , . The selected nodes only need to have the monotonically nondecreasing property.
Based on the continuity of the B-spline surface, the spline surface formed by uniform B-spline basis functions has a -order continuity [25] of , in the element boundary. The above manifold element adopts the model of the generalized node distribution, which is the same as that on the B-spline surface, such that second-order continuity exists in the element boundary. Figure 3 shows that the B-spline has a desirable property of local support, and its values are small within the interval to ensure the reasonability of the calculated condition number of the stiffness matrix. The number of generalized degrees of freedom required by the calculation is , which is significantly reduced compared with that of the conventional manifold method. Figure 2 shows that the B-spline basis functions have no property function in the nodes. Thus, these functions can reflect the dynamic property of the structure.
4. Manifold Scheme of the Bending Vibration of Simply Supported Thin Plates
4.1. Dynamic Equation of the Manifold Element Based on the Principle of Minimum Instantaneous Potential Energy
Suppose that a thin plate has a thickness of , a density of , an elasticity modulus of , a bending deflection (displacement) of , an actuating force in the unit area on the plate surface of , a normal direction for the simply supported boundary of , and corresponding fixed deflection and bending moments of and . According to the Kirchhoff plate theory [26] and based on the processing method of the boundary condition of the manifold element [1], the principle of the minimum instantaneous potential energy of the numerical manifold of the bending vibration of thin plates can be derived as where is the area of the thin plate, is the simply supported boundary, is the damping coefficient, and and are the stiffness coefficients of the penalty function. is the bending stiffness matrix, and is the curvature and torsion of the elastic surface represented by vectors in the directions of , that meets the following equations:
in Formula (11) can be expressed as the following equation by the displacement function, where is the angle from the direction to the direction:
The manifold element of Formula (8) is substituted into Formula (11) along with Formulas (12) and (14) to obtain the discretization form of the minimum instantaneous potential energy:
The first variation is taken from Formula (15) that is the instantaneous variation that satisfies and . Therefore, the dynamic equation for any manifold element is given by where
To facilitate the integration of the stiffness matrix, the transformation matrix is introduced and can be defined by the following formula, where is the total freedom of degree vector:
The transformation matrix corresponding to Formula (8) is where is the identity matrix.
Formula (18) is used to integrate Formula (16) into the overall dynamic equation where
It can be easily observed that (20) has an analogous form as the overall dynamic equation generated by FEM. Meanwhile, considering that and are symmetric positive matrices, we can solve the natural frequencies of the system by use of the techniques employed in FEM.
4.2. Initialization of the Dynamic Equation
Unlike the finite element method, , , and in the manifold method are no longer relevant to the nodes and thus need to be initialized. If the manifold element to be calculated is and the number of degrees of freedom is , the noncoincident points need to be regenerated in the directions of , to obtain the final nodes , , and . If the manifold element corresponding to node is and the real initial displacement is initialized with the velocity of and , we obtain where
Formulas (23) can be used to obtain the generalized displacement and the generalized velocity , which are substituted into Formula (20) to obtain the generalized initial acceleration:
4.3. Integration Scheme of the Dynamic Equation
When the domain of area integration is a regular rectangular element, Gaussian integration can be adopted. When the domain is an irregular geometric configuration, it can be subdivided into several triangular elements. Then, the simplex integration or Hammer integration is used to obtain the calculation results with sufficient accuracy [27]. When the boundary integration is a regular curve, Gaussian integration can be directly performed. If the boundary integration is a complex curve, it should be subdivided into several, straight-line segments before Gaussian integration is performed.
For time integration, the Wilson- method can be selected in this paper. Formula (20) is then converted into the linear equation T at time . Consider where
After solving , the generalized total displacement, velocity, and acceleration are solved:
Formulas (18) and (8) are used to obtain the function of the displacement, velocity, and acceleration of any manifold element at time :
5. Calculation Examples and Numerical Analysis
5.1. Calculation Example
As shown in Figure 4, the thin plate that is simply supported on the four edges has a length of m, a width of m, and a thickness of m. The transverse dynamic load is Pa, Pa, and . In addition, the elasticity modulus is GPa, the Poisson’s ratio is , the stiffness coefficient of the penalty function is , the material density is kg/m3, the damping coefficient is , and the coefficient of the Wilson- method is . Different time steps are selected for the calculation. As shown in Figure 5, the fixed rectangular mesh is used as the mathematical cover. The shaded area in the figure denotes any calculated manifold element.
With the theory illustrated in Section 4.1, the generalized instantaneous potential energy for this example can be expressed as where the latter four terms are the generalized energy form of boundary conditions.
After a large amount of discretization calculation, the matrices , , , and can be calculated by where and are given in (12) and (13), respectively.
The frequency and the displacement for this case are as follows:
5.2. Numerical Analysis
Table 1 shows the frequency of the manifold method, the accuracy of which will improve along with the increasing density of meshes. The relative error of frequency in the table is defined as . When the mesh is adopted, the first-order frequency error is only 0.0041%, whereas the fifth-order frequency error is 0.24%.
, , and are defined as the numerical displacement, velocity, and acceleration, respectively. In the same way, , , and are the partial derivatives of numerical displacement. The corresponding relative errors have the same definitions as in Table 1. Table 2 shows the velocities and acceleration of the given points at certain time point. The accuracy of the numerical solution is high. Considering that the given interpolation functions has the same continuity and coordination for and , Table 3 only shows part of the partial derivatives of displacement. The numerical solution shows that both stress and strain in the manifold method have good accuracy.
Figures 6, 7, and 8 are, respectively, the displacements, velocities, and acceleration at the points of and for each moment. With the reduction in time step , the accuracy of the numerical solution improves. In particular, with the increase in time steps, smaller step lengths of can significantly improve the accuracy of the numerical solution. Figure 9 shows some partial derivations at the points of and . The figure presents the calculated values of the first several periods. As the meshes become denser, the calculation accuracy significantly improves. Therefore, the numerical solution is consistent with the theoretical solution.
Since the manifold method in this paper adopts an overall interpolation approximate, the integrated error norm has to be introduced to reflect the interpolation accuracy fully. is defined as the integrated error of displacement, where is the number of discrete points equally distributed in the and directions within the thin plate. During calculation, is taken as discrete points within the thin plate. The integrated errors of the other quantities in Tables 4 and 5 are defined by a similar method. Table 4 shows the bilinear Hermite finite element. The mathematical software MATLAB is strongly recommended for the calculation, similar to the manifold method. By comparing Tables 4 and 5, the overall degrees of freedom for the manifold method involving the same number of elements are significantly smaller than those for the finite element method, but its calculation accuracy is higher than that of the latter. Furthermore, Table 6 gives the global errors from the cubic B-spline finite element. By comparing Tables 6 and 5, we can safely conclude the cubic B-spline manifold element produces more accurate results than cubic B-spline finite element. All numerical results confirm the validity of the manifold method in solving structural dynamics problems.
6. Conclusions
In this paper, based on the nonuniform cubic B-spline basis function, a new numerical manifold method was derived. By using the same structural model as the B-spline surface, the second-order coordinate can be obtained at the boundary of the manifold element. The local support property of the B-spline basis function enables the good performance of the stiffness matrix, which facilitates the numerical calculation. The new method was used to solve the bending problems of thin plates. The linear elasticity dynamic equation was derived based on the principle of the minimum instantaneous potential energy. The method for the initialization of the dynamic equation and the time integration method were also derived. The fixed mathematical mesh was used to implement the manifold method. The numerical results agreed well with the theoretical solution. Compared with the finite element method, the overall degrees of freedom of the cover calculated with the new manifold method were significantly reduced. The calculation accuracy was obviously improved, which indicates the superiority of the manifold method in solving structural dynamics problems. The theoretical and application studies in this paper are preliminary, confined to linear elasticity problems, and limited to situations in which only the time integration method is adopted. Further research on nonlinear problems should be conducted in the future to perfect the theory.
Conflict of Interests
The authors declare that they have no conflict of interests.
Acknowledgments
This work is financially supported by the Natural Science Foundation of China (11302272), the Natural Science Foundation of China (11102169), and the Science and Technology Project Affiliated with the Education Department of Chongqing Municipality (KJ120801).