Home Journals Articles Updates About
CN
image
收藏切换
A variationally consistent nodal integration for cubic serendipity finite elements with optimal convergence in explicit transient heat conduction analysis
收藏切换
PDF
Songyang Hou, Zhiwei Lin, Zhenyu Wu, Dongdong Wang*
Theoretical and Applied Mechanics Letters | 2025, 15(6) : 100616
Less
收藏切换
Theoretical and Applied Mechanics Letters | 2025, 15(6): 100616
• Research Article •
A variationally consistent nodal integration for cubic serendipity finite elements with optimal convergence in explicit transient heat conduction analysis
Full
Songyang Hou, Zhiwei Lin, Zhenyu Wu, Dongdong Wang*
Affiliations
  • Fujian Key Laboratory of Digital Simulations for Coastal Civil Engineering, Department of Civil Engineering, Xiamen University, Xiamen 361005, China
Published: 2025-11-01 doi: 10.1016/j.taml.2025.100616
Outline
收藏切换

The 13-node quadrilateral and 39-node hexahedral cubic serendipity elements produce nodally integrated positive-definite lumped heat capacity matrices in higher-order finite element analysis. However, these elements display severe convergence deterioration in explicit transient heat conduction analysis with lumped heat capacity matrices. This convergence decay is due to the violation of variational integration consistency by the standard Galerkin formulation with lumped heat capacity matrices. This issue is resolved by introducing the boundary-enhanced Galerkin weak form that incorporates the elemental boundary contribution in the discrete finite element formulation. Subsequently, it is theoretically proven that a direct nodal integration identically fulfills the variational integration consistency in the context of the boundary-enhanced Galerkin weak form. The proposed variationally consistent nodal integration therefore enables optimal convergence for explicit transient heat conduction analysis with lumped heat capacity matrices. The efficacy of the proposed variationally consistent nodal integration formulation for the 13-node quadrilateral and 39-node hexahedral cubic elements is thoroughly demonstrated via numerical examples.

Cubic serendipity element  /  Heat conduction analysis  /  Nodal integration  /  Lumped heat capacity matrix  /  Variational integration consistency  /  Optimal convergence
Songyang Hou, Zhiwei Lin, Zhenyu Wu, Dongdong Wang. A variationally consistent nodal integration for cubic serendipity finite elements with optimal convergence in explicit transient heat conduction analysis[J]. Theoretical and Applied Mechanics Letters, 2025 , 15 (6) : 100616 - . DOI: 10.1016/j.taml.2025.100616
Full Text
Less
1.　Introduction
Less
The serendipity finite elements [1–5] are designed to minimize the number of internal nodes and thus are popular for higher-order element formulation since they have fewer nodes than the same-order tensorproduct Lagrangian elements, for instance, the quadrilateral and hexahedral cubic serendipity elements considered in this study. Moreover, in transient finite element analysis [6–9], diagonal lumped mass matrices are usually employed with explicit time steps because of their trivial inversion cost in contrast to non-diagonal consistent mass matrices. However, particular attention is required for the mass lumping of serendipity elements because typical mass lumping schemes such as the row-sum method [1–3,10–12] and the nodal integration technique [13–16] often generate negative diagonal entries for these elements, rendering them unsuitable for transient analysis [17–19]. Although the widely used diagonal scaling HRZ (Hinton-Rock-Zienkiewicz) method [20] guarantees the positive definiteness of lumped mass matrices, HRZ diagonal lumped mass matrices actually do not exhibit optimal convergence [21,22]. On the other hand, the nodal quadrature approach takes the element nodes as the integration sample points and can yield optimally convergent transient analysis results if the integration precision of the nodal quadrature is sufficient [14,15]. Therefore, a special arrangement of non-uniform nodal locations is usually required for higher-order serendipity elements to improve the integration accuracy [15]. However, this process is often accompanied by negative lumped mass components. Although the heat capacity and conductivity matrices are the formal terms used in heat conduction analysis, for convenience of development, the mass and stiffness matrices are interchangeably used herein without confusion. Moreover, for conciseness, the capacity and conductivity matrices are also utilized in the subsequent discussion to replace the heat capacity and conductivity matrices.
The focus of this study is to develop optimally convergent quadrilateral and hexahedral cubic serendipity elements for explicit transient heat conduction analysis. The quadrilateral and hexahedral cubic serendipity elements originally consisted of 12 nodes for the twodimensional (2D) case and 32 nodes for the three-dimensional (3D) case. Nevertheless, these two elements cannot yield desirable positive-definite lumped heat capacity matrices via the nodal integration technique. This obstacle can be removed by introducing a few auxiliary internal or surface nodes, as has often been used in triangular elements [23–26], tetrahedral elements [27,28], and pyramid elements [29] to construct positive-definite lumped mass or capacity matrices. For the 2D 12-node quadrilateral element, one centroid node was inserted to formulate a 13-node quadrilateral element [1,30]. In this work, a 3D 39-node hexahedral element is further proposed by adding one body centroid node and six surface centroid nodes to the 32-node hexahedral element. The use of the nodal quadrature rules for these 13-node and 39-node cubic serendipity elements naturally yields positive-definite lumped capacity matrices. However, a numerical implementation of the 2D 13-node and 3D 39-node elements with lumped capacity matrices in the standard Galerkin setting still fails to attain the optimal convergence rates in explicit transient heat conduction analysis, which are 4 and 3 regarding L2 and H1 errors for cubic elements. As shown in this study, the fundamental reason behind these numerical observations is that the combination of the nodally integrated lumped capacity matrices with the standard Galerkin formulation does not meet the variational integration consistency, which implies that it cannot pass the corresponding higher-order patch tests.
To fulfill the variational integration consistency by the nodal quadrature, the finite element discrete Galerkin weak form is reformulated herein by incorporating the element boundary contribution into the conductivity matrix, which is referred to as the boundary-enhanced weak form for convenience. This weak form was adopted in [31] to construct lumped mass matrices in structural dynamics for tensor-product Lagrangian elements via the Gauss–Lobatto quadrature. One important fact is that the Gauss–Lobatto quadrature can actually yield optimal convergence in the standard Galerkin formulation; thus, linking the Gauss–Lobatto quadrature with the boundary-enhanced weak form does not appear necessary for the tensor-product Lagrangian elements. In contrast, in this work, the nodal quadrature is inherently merged with the boundary-enhanced weak form to realize variationally consistent nodal integration (VCNI) for cubic serendipity finite elements. More importantly, it is theoretically proven that the proposed VCNI identically satisfies the variational integration consistency and consequently ensures the pass of third-order patch tests and optimal convergence in explicit transient heat conduction analysis.
The rest of this paper is outlined as follows. In Section 2, the cubic serendipity elements, standard Galerkin formulation and temporal discretization are briefly summarized. Moreover, the variational integration consistency for both transient and steady analyses is particularly elaborated. In Section 3, the violation of variational integration consistency by the standard Galerkin formulation with a nodally integrated lumped capacity matrix is disclosed. In Section 4, a VCNI is subsequently proposed on the basis of boundary-enhanced weak form, and it is proven that VCNI identically meets the variational integration consistency criterion and exactly passes third-order patch tests. The effectiveness of VCNI is demonstrated in Section 5 through numerical results. Finally, conclusions are given in Section 6.
2.　Cubic serendipity elements and standard Galerkin formulation
Less
2.1.　 Basics of cubic serendipity elements
Cubic serendipity elements are often used in higher-order finite element analysis. As shown in Figs. 1(a) and 2(a), the conventional quadrilateral and hexahedral cubic serendipity elements have 12 and 32 nodes, where the element parametric domain is and nsd is the number of spatial dimensions. However, their lumped capacity matrices encompass negative entries for the corner nodes, which are marked by the blue solid dots in Figs. 1 and 2. To remove these negative lumped capacity components, 2D 13-node quadrilateral and 3D 39-node hexahedral serendipity elements with uniform nodal distributions are constructed by adding one centroid node to the 12-node element [30], and one body centroid node plus six surface centroid nodes to the 32-node elements. As described in Figs. 1(b) and 2(b), the 2D 13-node quadrilateral and 3D 39-node hexahedral serendipity elements produce positivedefinite lumped capacity matrices via the nodal integration approach, which are subsequently utilized for transient heat conduction analysis in this study.
As shown in Fig. 1(b), the 2D 13-node quadrilateral element has four corner nodes, eight mid-side nodes and one internal node, whose shape functions are given by:
where ξ = (ξ, η) represents the 2D parametric coordinate system.
The 3D 39-node hexahedral element depicted in Fig. 2(b) consists of 8 corner nodes, 24 mid-side nodes, 6 surface centroid nodes and one body centroid node. The shape functions of the 8 corner nodes are as follows:
where ξ = (ξ, η, ζ) denotes the 3D parametric coordinate system. For the 24 mid-side nodes, the shape functions read:
Besides, the shape functions of the 6 surface centroid nodes as well as the body centroid node have the following expressions:
2.2.　 Standard Galerkin formulation with cubic serendipity elements
In this study, transient heat conduction analysis is considered a model problem. For a multidimensional temperature field variable T(x, t) with spatial domain Ω, boundary Γ and time interval [0, 𝒯], the strong form of transient heat conduction analysis can be expressed as:
where Δ = ▽ · ▽, ▽ represents the gradient operator, c is the specific heat capacity, ρ denotes the material density, κ is the thermal conductivity coefficient, and s is the heat supply rate. and stand for the prescribed heat flux and temperature values on the natural boundary ΓN and essential boundary ΓD, respectively. n represents the outward normal vector of ΓN. The overhead dot denotes the time derivative, and T0 symbolizes the initial temperature.
The strong form in Eq. (6) can be transformed into the following weak form:
Following the standard finite element discretization procedure, the spatial domain Ω is partitioned into a set of 13-node or 39-node cubic serendipity elements, i.e., , where Ωe is the element domain, nel is the total number of elements. After the finite element discretization, Eq. (7) can be recast as:
where ñnel denotes the total number of elements pertaining to the natural boundary. By using the shape functions in Section 2.1, the finite element approximants of T, δT, and T in a generic element of Ωe read:
where nen is the number of nodes per element, denotes the approximate nodal temperature. Moreover, Eq. (9) is supplemented by the isoparametric geometry mapping:
Substituting Eq. (9) into Eq. (7) yields:
with
where M, K are the global heat capacity and conductivity matrices, respectively;, d, and f are the global temperature and heat supply vectors; A is the local-global assembly operator; the elemental quantities Me, Ke, , and are given by:
with
When the time effect is omitted, Eq. (11) reduces to the steady case:
In computer implementation, the quantities in Eq. (13) are usually evaluated via numerical integration:
where J = det(J), J is the Jacobian matrix; JS is the boundary Jacobian; and refer to the kth integration sample point and the related weight; and denote the boundary integration sample point and weight; nΩ and nΓ are the total numbers of domain and boundary integration points, respectively.
2.3.　 Temporal discretization with Heun’s method
Since the nodal integration and lumped heat capacity matrices are emphasized in this study, it is natural to employ an explicit timestepping method. The inversion of these nodally integrated lumped heat capacity matrices is computationally trivial, which enables fully explicit and efficient time stepping. Among various explicit schemes, the secondorder accurate Heun’s method [7] effectively balances computational efficiency and accuracy and is thus used herein for temporal discretization. Along this path, the time domain [0, 𝒯] is subdivided into nt time intervals, i.e., [tn, tn+1], n = 0, 1, …, nt – 1, and the time increment is denoted as Δt = tn+1 – tn. The governing equation in Eq. (11) at tn can be expressed as follows:
where dn, vn, and fn denote the nodal temperature, temperature rate, and heat supply vectors at tn, respectively.
According to Heun’s method [7], the solution at tn+1 is predicted by the forward Euler formula:
Subsequently, the solution dn+1 is corrected by the trapezoid rule:
2.4.　 Variational integration consistency
The variational integration consistency has been extensively discussed in meshfree methods [32–38] because the rational meshfree shape functions and overlapping support domains render numerical integration nontrivial. The consistency of variational integration implies that any solution that can be spanned by the basis or shape functions should be exactly reproduced by the corresponding Galerkin numerical solution, which can also be interpreted as the patch tests up to the highest degree of shape functions. Hence, for the time-dependent transient heat conduction analysis, to meet the variational integration consistency, the temperature T(x, t) can be assumed to be the following complete pth degree polynomial:
where α is a positive constant, denotes the set of nonnegative integers, cijk refers to arbitrary constants, and p represents the complete degree of shape functions.
In the discrete sense, the variational integration consistency associated with certain numerical quadrature rules can be straightforwardly verified by setting the nodal values and source term consistent with the given polynomial solutions defined by Eq. (20). Accordingly, the nodal temperature and its time derivative become:
Meanwhile, Eq. (20) corresponds to the following source and boundary flux terms:
Through substituting Eq. (22) into Eq. (16) and employing a specific numerical quadrature rule, M, K, and f are then readily attained. Consequently, the variational integration consistency can be conveniently examined by checking the vanishing of the following residual vector r(t):
where it is noteworthy that all the quantities on the right-hand side are consistent with the assumed field given by Eq. (20). In accordance with Eq. (23), a quadrature rule is said to satisfy the variational integration consistency in the case of r(t) = 0; otherwise, it violates this condition. Moreover, the above procedure degenerates to the assessment of variational integration consistency associated with steady problems by setting α = 0 [36,37].
3.　Violation of variational integration consistency in standard Galerkin formulation with nodally integrated lumped capacity matrix
Less
In this Section, the impact of the nodally integrated lumped capacity matrix on the consistency of variational integration is theoretically elucidated. The elemental heat capacity matrix Me in Eq. (13) is typically non-diagonal when it is integrated exactly. On the other hand, when nodal integration is employed, Me reduces to a lumped matrix due to the Kronecker delta property of shape functions:
where ξa is parametric coordinate of the ath shape function. wa denotes the corresponding integration weight, which can be readily determined by exactly integrating the ath shape function over the parametric element domain:
Substituting Eq. (1) into Eq. (25) yields the following integration weights for the 13-node quadrilateral serendipity element:
Similarly, introducing Eqs. (2)–(5) into Eq. (25) gives the following weights for the 39-node hexahedral serendipity element:
In subsequent development, the nodal integration rules described by Eqs. (24)–(27) are used to generate the lumped heat capacity matrices for the 13-node quadrilateral and 39-node hexahedral elements. For brevity, “GFI” and “GNI” are adopted to represent the standard Galerkin formulations using full Gauss and nodal quadrature rules for the numerical integration of the conductivity matrix Ke and the heat supply vector fe, where the lumped capacity matrix is employed for both circumstances. Notably, regarding the boundary integration in Eq. (13) for , GFI and GNI signify the use of 4-point Gauss and 4-point Newton–Cotes quadrature rules [14]. To examine whether GFI and GNI satisfy the variational integration consistency, a parallelogram element in Fig. 3 is considered. The reason for choosing such a generic parallelogram element is that it gives a constant Jacobian matrix, which is also the basis for higher-order elements to pass the patch tests exactly [39,40].
For 2D transient heat conduction analysis, T(x, t) in Eq. (20) is the complete cubic polynomial:
Without loss of generality and for clarity, the term of c30x3 exp(–αt) in Eq. (28) is taken as an example, and it corresponds to the following element nodal temperature and its time derivative:
Moreover, from Eqs. (22) and (28), and s(x, t) are obtained as:
To evaluate the residual in Eq. (23), M, K, and f of this parallelogram element need to be explicitly quantified. From Eq. (10), the Jacobian matrix of the 13-node element in Fig. 3 can be easily obtained as:
By substituting Eqs. (26) and (32) into Eq. (24), the nodally integrated lumped capacity matrix M of the parallelogram domain is obtained as:
Moreover, the conductivity matrix K and heat supply vector f will be elucidated next on the basis of different integration schemes.
3.1.　 Violation of variational integration consistency for GFI
In accordance with Eqs. (16) and (32), the corresponding K and f can be computed via the full Gauss quadrature rule, which refers to the 4-point Gauss integration in each spatial direction. For brevity, the resulting components of the first row for K are listed as follows:
with
Similarly, introducing Eq. (31) into Eq. (16) together with the full Gauss integration scheme yields the first component f1(t) of f(t):
Based upon Eqs. (29), (30), and (33)–(36), the first component r1(t) of the residual vector r(t) becomes
where M1 denotes the first row of the lumped capacity matrix M. Since r1(t) in Eq. (37) already fails to vanish, it can be concluded that GFI does not meet the variational integration consistency for transient heat conduction analysis. Consequently, GFI cannot achieve optimal convergence in explicit transient heat conduction analysis.
In the meantime, the time-dependent transient residual in Eq. (37) reduces to a time-independent steady residual by setting α = 0:
A repeat of the above procedure further gives:
Furthermore, Eqs. (38) and (39) hold for all the monomial terms in Eq. (28). Therefore, GFI does fulfill variational integration consistency for steady analysis, although it violates this condition in transient analysis when the nodally integrated lumped capacity matrix is used.
3.2.　 Violation of variational integration consistency for GNI
Following an almost identical procedure as that in Section 3.1, in the case of GNI where nodal integration is used to compute K, the first row of K, namely, K1, is obtained as:
with
Introducing the nodal quadrature of Eq. (31) into Eq. (16) yields f1(t) = 0. As a result, the transient residual of GNI becomes:
When α = 0, Eq. (42) gives the steady residual of GNI:
Equations (42) and (43) immediately show that GNI for the 13-node quadrilateral serendipity element is not variationally consistent in both transient and steady scenarios. Therefore, in the context of the standard Galerkin formulation with the 13-node quadrilateral serendipity element, neither GFI nor GNI meets the variational integration consistency requirement for transient heat conduction analysis. In addition, it is straightforward to verify that this conclusion is congruously applicable for the 3D 39-node serendipity element. A direct consequence of the failure to satisfy the variational integration consistency is that the optimal convergence of numerical solutions is lost in the transient heat conduction analysis.
4.　Variationally consistent nodal integration for cubic serendipity elements
Less
4.1.　 Boundary-enhanced weak form for cubic serendipity elements
In the previous Section, it was shown that both GFI and GNI, which are based on the standard Galerkin weak form, violate the variational integration consistency for transient heat conduction analysis. To resolve this issue, the Galerkin weak form is reformulated to specifically incorporate the element boundary contribution in the case of nodal integration. To this end, substituting Eq. (9) into Eq. (8) gives the finite element approximation of the standard Galerkin weak form:
The first term on the left-hand side of Eq. (44) can be subsequently rephrased as follows:
where ne is the outward normal on the element boundary Γe. Equation (45) indicates that each element introduces an additional flux term involving an integration over the element boundary Γe, thus, this formulation is referred to as the boundary-enhanced weak form for convenience. Obviously, the boundary-enhanced weak form only changes the conductivity matrix. By plugging Eq. (9) into Eq. (45), the element conductivity matrix Ke now reads:
with
It is noted that Eqs. (46) and (47) involve the second-order derivatives of the shape functions. When an exact integration is used, the conductivity matrix in Eq. (46) is symmetric and identical to that in Eq. (13) arising from the standard Galerkin formulation. On the other hand, these two conductivity matrices are different when nodal integration is employed. In this case, the proposed VCNI yields a nonsymmetric conductivity matrix.
4.2.　 VCNI
The integration rules for the 2D 13-node quadrilateral and 3D 39-node hexahedral serendipity elements are illustrated in Figs. 4 and 5, respectively. In the 2D case, the element domain quadrature adopts the 13-node nodal integration in Section 3, and the element boundary adopts the 4-point Newton–Cotes integration rule [14]. Moreover, for the 3D 39-node hexahedral element, domain integration employs the element nodes as sample points, and the corresponding weights are illustrated in Fig. 5. For the element boundary Γe, the nodal integration in Eq. (26) is utilized because these nodal locations exactly match the surface boundary nodes of the 39-node hexahedral element.
On the basis of the nodal integration rules in Fig. 4 and following an identical procedure in Section 3, the following quantities for the 2D 13-node quadrilateral element in Fig. 3 can be obtained:
A combination of Eqs. (29), (30), (33), (48), and (49) then gives:
Similarly, it is trivial to verify that Eq. (50) holds for all the residual entries with respect to every term in Eq. (28):
Meanwhile, it is readily confirmed that Eq. (51) is true for all the monomial terms in Eq. (28). In addition, because Eq. (51) is independent of the parameter α, the relationship of ra = 0 still holds in the steady analysis that corresponds to α = 0.
After a slightly lengthy but straightforward derivation, the vanishing residual in Eq. (51) is valid for the 39-node hexahedral serendipity element via the nodal integration rules in Fig. 5. Therefore, in the context of the boundary-enhanced Galerkin weak form, the nodal integration schemes in Figs. 4 and 5 guarantee the consistency of variational integration and ensure the passing of higher-order patch tests and optimal solution convergence. Consequently, the nodal integration schemes in Figs. 4 and 5 are termed VCNI.
5.　Numerical examples
Less
In the subsequent numerical discussions, consistent with the previous nomenclature, VCNI signifies variationally consistent nodal integration in conjunction with the boundary-enhanced weak formulation, and GFI and GNI denote the standard Galerkin approaches that use the full Gauss and nodal quadrature to compute the heat conductivity matrix and supply vector. In all the cases, nodal integration is employed to evaluate the heat capacity matrix. Moreover, the following L2 and H1 error norms are adopted to assess the convergence of these methods:
Without loss of generality, dimensionless geometry and material parameters are used throughout the numerical examples.
5.1.　 Patch tests
To examine whether GFI, GNI, and VCNI satisfy the variational integration consistency, as shown in Fig. 6, the patch tests within the parallelogram (2D) and parallelepiped (3D) domains are considered, where the geometry parameter L and all material properties are set to unity values. In the steady-state analysis, the analytical temperature solutions are assumed as follows:
The boundary conditions are imposed in accordance with the analytical solutions provided in Eq. (53). For the 2D parallelogram problem discretized by four 13-node quadrilateral elements, the natural boundary conditions are imposed on the left and right sides, whereas the essential boundary conditions are prescribed on the bottom and top sides. For the 3D parallelepiped problem meshed by eight 39-node hexahedral elements, the natural boundary conditions are applied on the left and right surfaces, and the remaining surfaces are subjected to the essential boundary conditions. The steady patch test results are presented in Table 1 regarding L2 and H1 errors, and the temperature error contours are presented in Figs. 7 and 8, respectively. These results conclusively demonstrate that both GFI and the proposed VCNI successfully pass the 2D and 3D steady patch tests, whereas GNI fails the patch tests with noticeable errors.
Next, the transient patch tests are investigated by considering the following analytical solutions:
To exclude the errors introduced by the time integration, a sufficiently small time step Δt = 1 × 10–7 is adopted for both 2D and 3D computations. The corresponding transient analysis results at t = 1 are tabulated in Table 2, and the temperature error contours at t = 1 are simultaneously depicted in Figs. 9 and 10, which reveal that among the three methods of GFI, GNI and VCNI, only the proposed VCNI can pass the transient patch tests since it precisely fulfills the variational integration consistency requirement for explicit transient heat conduction analysis.
5.2.　 2D l-shaped domain problem
An l-shaped domain heat conduction problem is described in Fig. 11, where the geometry parameter is L = 2. The related material parameters are the material density ρ = 2700, specific heat capacity c = 900 and thermal conductivity coefficient κ = 205. For this heat conduction problem, the following manufactured analytical solution is considered:
To assess the algorithmic convergence, three progressively refined meshes with 13-node quadrilateral elements are illustrated in Fig. 12, and a time step of Δt = 0.1 is employed for the time integration based upon Eqs. (18) and (19). The convergence results of GFI, GNI, and VCNI at t = 30 are plotted in Fig. 13, in which suboptimal convergence rates of 3 and 2 for GFI and GNI are observed for L2 and H1 errors, respectively. In contrast, the optimal convergence rates of 4 and 3 in accordance with L2 and H1 errors are produced as expected by the proposed VCNI owing to its satisfaction with the consistency of variational integration. Moreover, the temperature distributions at t = 30 in Fig. 14 show that VCNI results agree very well with the analytical solutions in Fig. 11, and the corresponding solution errors are much smaller than those generated by GFI and GNI.
5.3.　 3D quarter hollow-cylinder problem
The final example is heat conduction analysis in a 3D quarter hollow cylinder. As shown in Fig. 15, the dimensionless geometry and material parameters are as follows: inner and outer radii ri = 1 and ro = 2, height H = 1, material density ρ = 7870, specific heat capacity c = 455 and thermal conductivity coefficient κ = 80. In the explicit transient heat conduction analysis, the following manufactured analytical solution is employed:
The finite element discretizations of this 3D problem using 39-node hexahedral elements are given in Fig. 16. In the finite element computation, as illustrated in Fig. 15, the essential boundary conditions are enforced on the top and bottom faces, and the remaining boundary surfaces are imposed by the natural boundary conditions. All the prescribed boundary conditions are evaluated based on the manufactured analytical solution in Eq. (56). The time step used in Heun’s method for this 3D transient problem is Δt = 1 × 10–2. The corresponding convergence results produced by different methods at t = 10 are plotted in Fig. 17, which testifies that the proposed VCNI guarantees optimal convergence rates of 4 and 3 for L2 and H1 errors. In contrast, the results of both GFI and GNI display severe convergence as well as accuracy deterioration. Moreover, the temperature solutions and the related errors at t = 10 in Fig. 18 further confirm the accuracy superiority of VCNI over GFI and GNI.
6.　Conclusions
Less
A VCNI was developed to achieve optimal convergence in the explicit transient heat conduction analysis via cubic serendipity finite elements, where the variational integration consistency was particularly elaborated. The cubic serendipity elements here specifically refer to the 2D 13-node quadrilateral element and the 3D 39-node hexahedral element, which are constructed by inserting a small number of auxiliary internal or surface nodes into the conventional 12-node and 32-node serendipity elements. Both 2D 13-node quadrilateral and 3D 39-node hexahedral elements yield desirable positive-definite lumped heat capacity matrices via nodal integration, which are preferable for explicit finite element computations, for example, Heun’s method. It was theoretically found that the standard Galerkin formulation with nodally integrated lumped capacity matrices loses the optimal solution convergence rates of 4 and 3 for L2 and H1 errors, and the underlying reason is the violation of variational integration consistency. Alternatively, the nodal integration for cubic serendipity finite elements was specifically implemented in the boundary-enhanced weak form that takes into account the element boundary contribution, as leads to the proposed VCNI. It was then proven that the proposed VCNI perfectly fulfills the variational integration consistency and then ensures optimally convergent numerical solutions for explicit transient heat conduction analysis. The optimal convergence and accuracy superiority of the proposed VCNI were consistently demonstrated by the numerical results.
References
Less
[1]
Szabó B.A., Babuška I., Finite Element Analysis, John Wiley & Sons, Hoboken, 1991.
[2]
Hughes T.J.R., The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover Publications, Mineola, 2000.
[3]
Zienkiewicz O.C., Taylor R.L., Zhu J.Z., The Finite Element Method: Its Basis and Fundamentals (7th Edition), Elsevier, Singapore, 2015.
[4]
Arnold D.N., Awanou G., The serendipity family of finite elements, Found. Comput. Math. 11 (2011) 337-344.
[5]
Arbogast T., Wang C., Direct serendipity finite elements on cuboidal hexahedra, Comput. Methods Appl. Mech. Eng. 433 (2025) 117500.
[6]
Duczek S., Gravenkamp H., Critical assessment of different mass lumping schemes for higher order serendipity finite elements, Comput. Methods Appl. Mech. Eng. 350 (2019) 836-897.
[7]
Ferziger J.H., Computational Methods for Fluid Dynamics (4th Edition), Springer, Switzerland, 2020.
[8]
Oñate E., de Pouplana I., Zárate F., Explicit time integration scheme with large time steps for first order transient problems using finite increment calculus, Comput. Methods Appl. Mech. Eng. 402 (2022) 115332.
[9]
Jiang X., Gan J., Teng S., Simulation of anchor chain based on lumped mass method, Theoret. Appl. Mech. Lett. 13 (2023) 100460.
[10]
Li X., Zhang H., Wang D., A lumped mass finite element formulation with consistent nodal quadrature for improved frequency analysis of wave equations, Acta Mech. Sin. 38 (2022) 521388.
[11]
Li X., Wang D., Xu X., Sun Z., A nodal spacing study on the frequency convergence characteristics of structural free vibration analysis by lumped mass lagrangian finite elements, Eng. Comput. 38 (2022) 5519-5540.
[12]
Yang Y., Zheng H., Sivaselvan M.V., A rigorous and unified mass lumping scheme for higher-order elements, Comput. Methods Appl. Mech. Eng. 319 (2017) 491-514.
[13]
Fried I., Malkus D.S., Finite element mass matrix lumping by numerical integration with no convergence rate loss, Int. J. Solids. Struct. 11 (1975) 461-466.
[14]
Stroud A.H., Approximate Calculation of Multiple Integrals, Prentice-Hall, Englewood Cliffs, 1971.
[15]
Li Y., Liang R., Wang D., On convergence rate of finite element eigenvalue analysis with mass lumping by nodal quadrature, Comput. Mech. 8 (1991) 249-256.
[16]
Duczek S., Gravenkamp H., Mass lumping techniques in the spectral element method: on the equivalence of the row-sum, nodal quadrature, and diagonal scaling methods, Comput. Methods Appl. Mech. Eng. 353 (2019) 516-569.
[17]
Malkus D.S., Plesha M.E., Zero and negative masses in finite element vibration and transient analysis, Comput. Methods Appl. Mech. Eng. 59 (1986) 281-306.
[18]
Malkus D.S., Plesha M.E., Liu M.R., Reversed stability conditions in transient finite element analysis, Comput. Methods Appl. Mech. Eng. 68 (1988) 97-114.
[19]
Voet Y., Sande E., Buffa A., A mathematical theory for mass lumping and its generalization with applications to isogeometric analysis, Comput. Methods Appl. Mech. Eng. 410 (2023) 116033.
[20]
Hinton E., Rock T., Zienkiewicz O.C., A note on mass lumping and related processes in the finite element method, Earthq. Eng. Struct. Dyn. 4 (1976) 245-249.
[21]
Hou S., Li X., Wang D., Lin Z., A mid-node mass lumping scheme for accurate structural vibration analysis with serendipity finite elements, Int. J. Appl. Mech. 13 (2021) 2150013.
[22.]
Li Hou X., Lin Z., Wang D., An aspect ratio dependent lumped mass formulation for serendipity finite elements with severe side-length discrepancy, Comput. Mech. 74 (2024) 819-847.
[23]
Cohen G., Joly P., Roberts J.E., Tordjman N., Higher order triangular finite elements with mass lumping for the wave equation, SIAM. J. Numer. Anal. 38 (2001) 2047-2078.
[24]
Mulder W.A., Higher-order mass-lumped finite elements for the wave equation, J. Comput. Acoust. 9 (2001) 671-680.
[25]
Mulder W.A., More continuous mass-lumped triangular finite elements, J. Sci. Comput. 92 (2022) 38.
[26]
Fan H., Huang D., Wang G., A new family of high-order triangular weight functions for mass lumping in vibration analysis, Int. J. Numer. Methods Eng. 124 (2023) 3796-3833.
[27]
Danielson K.T., Fifteen node tetrahedral elements for explicit methods in nonlinear solid dynamics, Comput. Methods Appl. Mech. Eng. 272 (2014) 160-180.
[28]
Geevers S., Mulder W.A., Van Der Vegt J.J.W., New higher-order mass-lumped tetrahedral elements for wave propagation modelling, SIAM J. Sci. Comput. 40 (2018) A2830-A2857.
[29]
Browning R.S., Danielson K.T., Littlefield D.L., Second-order pyramid element formulations suitable for lumped-mass explicit methods in nonlinear solid mechanics, Comput. Methods Appl. Mech. Eng. 405 (2023) 115854.
[30]
Babuška I., Strouboulis T., Upadhyay C.S., Gangaraj S.K., Computer-based proof of the existence of superconvergence points in the finite element method; superconvergence of the derivatives in finite element solutions of Laplace’s, Poisson’s, and the elasticity equations, Numer. Methods Partial. Differ. Equ 12 (1996) 347-392.
[31]
Schillinger D., Evans J.A., Frischmann F., Hiemstra R.R., Hsu M.C., Hughes T.J.R., A collocated Cº finite element method: reduced quadrature perspective, cost comparison with standard finite elements, and explicit structural dynamics, Int. J. Numer. Methods Eng. 102 (2014) 576-631.
[32]
Chen J.S., Wu C.T., Yoon S., You Y., A stabilized conforming nodal integration for Galerkin mesh-free methods, Int. J. Numer. Methods Eng. 50 (2001) 435-466.
[33]
Wang D., Sun Y., An efficient Galerkin meshfree formulation for shear deformable beam under finite deformation, Theor. Appl. Mech. Lett. 1 (2011) 051010.
[34]
Chen J.S., Hillman M., Rüter M., An arbitrary order variationally consistent integration for Galerkin meshfree methods, Int. J. Numer. Methods Eng. 95 (2013) 387-418.
[35]
Duan Q., Gao X., Wang B., Li X., Zhang H., Belytschko T., Shao Y., Consistent element-free Galerkin method, Int. J. Numer. Methods Eng. 99 (2014) 79-101.
[36]
Wang D., Wu J., An inherently consistent reproducing kernel gradient smoothing framework toward efficient Galerkin meshfree formulation with explicit quadrature, Comput. Methods Appl. Mech. Eng. 349 (2019) 628-672.
[37]
Wu J., Wang D., An accuracy analysis of Galerkin meshfree methods accounting for numerical integration, Comput. Methods Appl. Mech. Eng. 375 (2021) 113631.
[38]
Ying J., Wang D., Deng L., Lin Z., An immersed boundary fast meshfree integration methodology with consistent weight learning, Comput. Methods Appl. Mech. Eng. 429 (2024) 117121.
[39]
Taylor R.L., Simo J.C., Zienkiewicz O.C., Chan C.H., The patch test-A condition for assessing FEM convergence, Int. J. Numer. Methods Eng. 22 (1986) 39-62.
[40]
Zienkiewicz O.C., Taylor R.L., The finite element patch test revisited-A computer test for convergence, validation and error estimates, Comput. Methods Appl. Mech. Eng. 149 (1997) 223-254.
Appendix
Less
Year 2025 volume 15 Issue 6
PDF
56
34
Cite this Article
BibTeX
Article Info
doi: 10.1016/j.taml.2025.100616
  • Receive Date：2025-06-18
  • Online Date：2026-03-24
  • Published：2025-11-01
Article Data
Affiliations
History
  • Received：2025-06-18
  • Revised：2025-09-07
  • Accepted：2025-09-07
Affiliations
    Fujian Key Laboratory of Digital Simulations for Coastal Civil Engineering, Department of Civil Engineering, Xiamen University, Xiamen 361005, China

Corresponding:

*

E-mail address: (D. Wang).
References
Share
https://castjournals.cast.org.cn/joweb/taml/EN/10.1016/j.taml.2025.100616
Share to
QR

Scan QR to access full text

Cite this article
BibTeX
Citations
表12种不同金属材料的力学参数
科
Family		 属数
Number of
genus		 种数
Number of
species		 占总种数比例
Percentage of
total species (%)		 属
Genus		 种数
Number of
species		 占总种数比例
Percentage of total
species (%)
鹅膏菌科Amanitaceae 2 11 5.26 鹅膏菌属 Amanita 10 4.78
小菇科 Mycenaceae 2 12 5.74 丝盖伞属 Inocybe 5 2.39
多孔菌科 Polyporaceae 8 14 6.70 蜡蘑属 Laccaria 5 2.39
红菇科 Russulaceae 3 23 11.00 小皮伞属 Marasmius 6 2.87
小菇属 Mycena 11 5.26
光柄菇属 Pluteus 5 2.39
红菇属 Russula 17 8.13
栓菌属 Trametes 5 2.39
关闭全屏
  • BibTeX
  • EndNote
  • RefWorks
  • TxT
Links
Articles
Latest Articles
Most Read
Collections
Updates
Events
News
Multimedia
About
About Us
Contact
No. 86 Xueyuan South Road, Haidian District, Beijing
100081
010-62199257
qkjq@cast.org.cn

Copyright © 2025 China Association for Science and Technology. All rights reserved. For all open access content, the relevant licensing terms apply.
Sponsored by the Office of the Leading Group for Cybersecurity and Informatization of CAST, and supported by Science and Technology Review Publishing House