1.
Introduction
The Stokes equation [1] is a fundamental tool used to describe viscous fluid flow [2] at low Reynolds numbers (Re) [3], which typically indicates laminar flow conditions [4]. Re characterizes the ratio of inertial forces to viscous forces in fluid dynamics. When Re is very small, the characteristic velocity of the flow can be considered to approach zero. In this limit, the quadratic terms involving velocity in the Navier-Stokes equations become negligible. Consequently, the Navier-Stokes equations simplify the Stokes equations, helping to analyze more complicated fluid problems, with a very wide range of applications [5,6]. With the development of computer science, many numerical methods have been developed to solve Stokes problems, such as finite element method (FEM) [7,8,9], finite difference method [10], mixed FEM [11], boundary element method [12,13], and coupling of FEM [14]. Among them, the FEM has gradually become an important numerical computational method for approximating partial differential equations (PDEs) because of its many advantages such as strong program versatility, high accuracy, flexible mesh selection, and ability to deal with complex boundaries and high-order problems.
The FEM, as an important numerical method for solving mathematical [15] and physical problems [16], has been widely applied in the field of engineering mechanics [17]. In 1943, Courant [18] introduced the concept of FEM by using continuous functions on triangular regions to solve the torsion problem of St. Venant. By the mid-1960s, Feng [19,20] had independently established the mathematical foundation of FEM, making it a systematic and widely used numerical method. Since then, the scope of FEM's application has expanded from single structural analysis to various fields such as sound field analysis, flow field analysis, and electromagnetic field analysis. Based on the variational principle and subdivision interpolation, FEM uses interpolation functions in each element to approximate the unknown function in the domain piece by piece.
With the continuous progress of computer science, the FEM has undergone remarkable developments and improvements. Many new computational methods have emerged, including the finite volume method [21], upwind FEM [22], and spectral methods [23]. These new methods not only enrich the technical means of numerical simulation but also play a crucial role in improving computational efficiency and enhancing model accuracy. In the FEM, the selection of finite element basis functions is key, and appropriate test and trial function spaces ensure the accuracy and stability of the solution. Common choices include Lagrange functions [24], Hermite functions [25], Argyris functions [26], and Bernstein functions [27]. In 1979, Shi [28] used cubic B-spline variational methods to solve equilibrium problems in composite elastic structures of plates and beams in regular domains, introducing spline FEM. In the same year, Qin [29] proposed the spline finite point method based on spline functions, beam vibration functions, and energy variation. In 2005, Hughes et al. [30] used spline basis functions for approximate numerical calculations of field variables in physical problems in finite element analysis. In 2007, Bhatti and Bracken [27] proposed applying Bernstein polynomial bases to solve PDEs. Zhu and Kang [31], in 2010, used cubic B-spline quasi-interpolation to numerically solve the Burgers-Fisher equation. Dutykh [32], in 2014, solved the KdV and KdV-BBM equations using B-spline FEM. More recently in 2022, Pranta [33] solved 2D nonlinear parabolic PDEs using bicubic Bernstein polynomial bases. These developments highlight the ongoing evolution and versatility of FEM in addressing complex engineering and scientific challenges.
Lagrange interpolation functions are typically global, offering high accuracy in certain scenarios but potentially leading to numerical instability, especially with high-degree polynomials or complex boundary conditions. Although the Runge phenomenon [34] is less pronounced in FEM due to the integral approximation approach, it can still occur in specific cases, particularly with high-order polynomials or complex boundaries. Conversely, Bernstein polynomial functions have local support, enhancing numerical stability. They facilitate the construction of higher-order test and trial function spaces and are adept at handling complex boundary conditions. Additionally, Bernstein polynomials are non-negative and shape-preserving, making them uniquely suitable for shape-preserving approximations. Therefore, we have chosen Bernstein polynomial basis functions for our FEM implementation to ensure enhanced numerical stability and accuracy. However, it is important to note that Bernstein polynomials also have some limitations. While they offer stability and flexibility, the computational cost can increase significantly with higher polynomial degrees, making them less efficient for large-scale or real-time applications. Moreover, the theoretical foundation for certain specific problems, such as those involving very-high-order polynomials or highly oscillatory functions, may still require further research and development.
Consider the boundary value problem of the 2D Stokes equation,
where Ω⊂R2 is a bounded polygonal domain, u=u(x,y) is the velocity vector, p refers to fluid pressure, ∇⋅u is the divergence of u, and T(u,p)=2νD(u)−pI is the stress tensor. In more details, the deformation tensor can be written as
Hence, the stress tensor can be written as
f describes the external force, g is the velocity on the domain boundary, and ν=ULRe represents the kinematic viscosity of the fluid where U and L represent characteristic speed and characteristic length, respectively. The Stokes equation is a basic equation of fluid mechanics, which simulates the motion of low velocity or viscous fluid [35,36] and has important applications in fluid mechanics [37], geophysics [38], telecommunication technology [39], and aerospace [40], among others [41,42,43]. We use the mixed FEM based on Bernstein polynomial basis to solve Stokes equations, and calculate the errors of the L∞, L2, and H1-semi norms.
The traditional FEM is a versatile numerical technique that can handle both univariate and multivariate equations. However, when applied to systems involving multiple physical quantities, such as the Stokes equation (Eq (1.1)), traditional FEM requires careful consideration to ensure the existence and uniqueness of the solution. The Stokes equation involves a tight coupling between velocity and pressure, which necessitates precise numerical treatment. To guarantee the uniqueness [44] of the solution to the variational problem, the finite element approximation space must satisfy the Lax-Milgram theorem [45]. Additionally, to ensure the stability of the solution, especially for coupled variables, the inf-sup condition [45] must be satisfied. While traditional FEM can theoretically meet these requirements, the selection of appropriate finite element spaces for velocity and pressure is crucial. If the selection is not appropriate, the solution can become unstable and lose accuracy. Therefore, to ensure that the finite element approximation for the Stokes equation is both convergent and stable, we have chosen the mixed FEM. The mixed FEM can better handle the coupling between velocity and pressure by selecting suitable finite element spaces for these variables. This approach more effectively satisfies the inf-sup condition, thereby providing a more stable and accurate solution. Besides, we found that only the gradient term of pressure appeared in the Stokes equation, which cannot guarantee that the solution of pressure is unique. Therefore, in the process of solving, we need to impose additional conditions for pressure. In this article, we fix pressure at one point in the region. Furthermore, the mixed FEM is not limited to the Stokes equation. It can be equally effective in other multivariate systems, such as the velocity-stress formulation of the wave equation [46,47]. By using mixed FEM, we can achieve higher accuracy and stability in solving a wide range of coupled PDEs.
This paper is organized as follows. In Section 2, we first review some basic contents of Bernstein polynomial basis, Bézier curves, and surfaces. In Section 3, we use the mixed FEM based on the Bernstein polynomial basis to derive the discrete scheme of Stokes equation. In Section 4, the error result is obtained by some numerical examples. In Section 5, we summarize the work.
2.
Bernstein basis functions, curves, and surfaces
In this section, we will recall the definitions and properties of Bernstein polynomial bases, Bézier curves, and surfaces.
Definition 1. Bernstein polynomial bases of degree n are defined by
where, (ni)=n!i!(n−i)!,i=0,1,⋯,n. For simplicity, when i<0 or i>n, let Bni(x)=0,x∈[0,1].
Definition 2. Given control points Pi(x,y)∈R2(i=0,1,⋯,n), the Bézier curve of n degrees is defined by
where Bni(x)(i=0,1,⋯,n) is defined as Eq (2.1), and the n-edge polygon obtained by connecting two adjacent control points with straight line segments is called the control polygon.
Bernstein polynomial bases of tensor product type can be obtained by tensor product from Bernstein polynomial bases of one variable.
Definition 3. The tensor product type Bernstein polynomial bases of m×n degree are defined by
Definition 4. For a continuous function f(s,τ) defined on [0,1]×[0,1], the tensor product Bernstein polynomial interpolation operator Bh is defined as
Next, we prove that Bh is a bounded interpolation operator.
Since f is a continuous function, it is bounded on [0,1]×[0,1]. Therefore, there exists a constant M such that |f(s,τ)|≤M for all (s,τ)∈[0,1]×[0,1]. Hence
So we can get that Bh is a bounded interpolation operator.
Since Bh is a bounded interpolation operator, by the Bramble-Hilbert Lemma [45], we can conclude that
3.
Mixed finite element algorithm based on Bernstein polynomial basis for Stokes equations
In this section, we first construct function spaces of the mixed FEM with Bernstein polynomial basis, and the discrete scheme of Stokes equation in Eq (1.1) is derived.
3.1. Galerkin formulation
First of all, consider the subspace H^{1}_0(\Omega) of Sobolev space H^{1}(\Omega) :
Multiplying the first equation of Eq (1.1) by test vector function \mathbf{v}(x, y) \in H^{1}_0(\Omega)\times H^{1}_0(\Omega) and then integrating on \Omega yields
Second, by multiplying the divergence-free equation by a test function q(x, y) , we get
Then, applying Green's identity,
where,
Introducing bilinear form,
Then, the variational formulation of the mixed FEM of Eq (1.1) is to find \mathbf{u}\in H_0^1(\Omega)\times H_0^1(\Omega) and p\in L^2(\Omega) , satisfying the following equation
for any \mathbf{v}\in H_0^1(\Omega)\times H_0^1(\Omega) and q\in L^2(\Omega) , where (\mathbf{f}, \mathbf{v}) = \int_{\Omega}\mathbf{f}\cdot\mathbf{v}\:dxdy.
Then, we consider the discrete form of variational Eq (3.1).
3.2. Finite element discretization
Let \Omega_{h} be a uniform rectangle partition of \Omega , \mathbf{h} = [h_1, h_2] = [\frac{1}{N_1}, \frac{1}{N_2}] is the mesh size, where N_1 and N_2 represent the number of subintervals on the x -axis and y -axis of quasi-uniform subdivision. For each T\in \Omega_{h} , the local finite element space Q(T, m, n) is spanned by Bernstein polynomial basis defined on T , i.e.,
Consider a finite element space U_h(m, n)\subset H^1(\Omega) for the velocity \mathbf{u} and a finite element space W_h(h, l)\subset L^2(\Omega) for the pressure p . Assume that the polynomial space in the construction of U_h contains P_k, k\geq 1 and that of W_h contains P_{k-1} , where,
Define U_{h0} to be the space that consists of the functions of U_h with vanishing boundary values.
Subsequently, the discrete scheme of Eq (3.1) is to find \mathbf{u}_h\in U_h\times U_h and p_h\in W_h , where \mathbf{u}_h = (u_{1h}, u_{2h}) such that
for any \mathbf{v}_h\in U_{h0}\times U_{h0} and q_h\in W_h .
In order to verify if \mathbf{v}_h and q_h satisfy the inf-sup condition, we now define an interpolation \pi_h so that it is a modification of B_h , that is, satisfying \pi_h\mathbf{u} = B_h(\mathbf{u}) . From (2.4), we know \pi_h is bounded. Discrete compatibility is similar to that proved in [48], that is, b(\mathbf{u}-\pi_h\mathbf{u}, q) = 0 . So, \mathbf{v}_h and q_h satisfy the following inf-sup condition:
where \beta > 0 is a constant independent of mesh size h .
In the scalar format, Eq (3.2) is to find u_{1h} \in U_h, u_{2h} \in U_h, and p_h \in W_h such that
for any v_{1h} \in U_h, v_{2h} \in U_h, q_h \in W_h.
Since u_{1h}, \:u_{2h}\in U_h = span\{r_j\}_{j = 1}^{N_b} and p_h\in W_h = span\{w_j\}_{j = 1}^{N_{bp}} , then
for some coefficients u_{1j}, \:u_{2j} (j = 1, \cdots, N_b) , and p_j (j = 1, \cdots, N_{bp}).
Now, we set up a linear algebraic system for u_{1j}, u_{2j} (j = 1, \cdots, N_b) and p_j (j = 1, \cdots, N_{bp}). Then we can solve it to obtain the finite element solution \mathbf{u}_h = (u_{1h}, u_{2h})^t and p_h .
For the first equation in the Eq (3.3), in the first set of test functions, we set \mathbf{v}_h = (r_i, 0)^t , namely, v_{1h} = r_i(i = 1, \cdots, N_b) and v_{2h} = 0 . Then
After that, we let \mathbf{v}_h = (0, r_i)^t , i.e., v_{1h} = 0 and v_{2h} = r_i(i = 1, \cdots, N_b) ,
Lastly, set q_h = w_i(i = 1, \cdots, N_{bp}) in the second equation of the Eq (3.3), getting
Simplify the above three sets of equations, obtaining
Define the stiffness matrix as
where \mathbb{O} = [0]_{i = 1, j = 1}^{N_{bp}, N_{bp}},
Define the load vector
where \vec{0} = [0]_{i = 1}^{N_{bp}} ,
Define the unknown vector
where,
Then, we get a linear system of ordinary differential equations for u_{1j}, u_{2j}(j = 1, \cdots, N_b) and p_j(j = 1, \cdots, N_{bp}) ,
so we can solve system (3.4) and obtain the unknown vector group \vec{X} .
4.
Numerical examples
In this section, we verify the feasibility and effectiveness of this method using several numerical examples. Tensor product Bernstein polynomial bases are used to construct the trial function space and test function space of the mixed FEM, the approximate solutions are solved by MATLAB 2022 b, and the error and convergence order of the exact solution and the finite element solution under L^{\infty} , L^{2} , and H^{1} -semi norms are obtained. The numerical results obtained by solving Stokes equation with bilinear, biquadratic, and bicubic Lagrange basis functions are consistent with those obtained by using Bernstein polynomial basis with corresponding orders. Since using Lagrange basis functions of higher than bicubic order leads to the Runge phenomenon when solving the Stokes equations, we only present the error results of Bernstein polynomial basis.
Example 1. Consider the following two-dimensional stokes equation with Dirichlet boundary in rectangular domain \Omega = [0, 1]\times [0, 1] .
where the exact solution \mathbf{u} = (u_1, u_2)^t is
the exact solution p(x, y) is
and the body force \mathbf{f} = (f_1, f_2)^t is
where we set \nu = 1 .
The domain \Omega is partitioned into uniform rectangles. Here, we use biquadratic, bicubic, and biquartic Bernstein polynomial basis to solve the Stokes Eq (4.1), and calculate the L^{\infty} , L^{2} , and H^{1} -semi norms between the approximate solution and the exact solution. Tables 1 and 2 show the numerical errors for these kinds of basis functions in L^{\infty} , L^{2} , and H^{1} -semi norms; the corresponding convergence orders are shown in Tables 3 and 4. The comparison of errors are shown in Figures 1 and 2.
When solving 2D Stokes equations, with equal mesh sizes, for velocity \mathbf{u} , the numerical accuracy of the bicubic Bernstein polynomial basis is 1 and 2 orders of magnitude higher than that of the biquadratic Bernstein polynomial basis, while the biquartic Bernstein polynomial basis is 4–7 orders of magnitude higher than the bicubic. For pressure p, the numerical accuracy of the biquadratic Bernstein polynomial basis is 3–5 orders of magnitude higher than that of the bilinear Bernstein polynomial basis, and the bicubic Bernstein polynomial basis is 1–5 orders of magnitude higher than the biquadratic.
When solving Eq (4.1) using biquartic Bernstein polynomial basis for velocity \mathbf{u} and bicubic Bernstein polynomial basis for pressure p , we attempted many methods, including adjusting the accuracy setting of MATLAB 2022 b to improve the accuracy and convergence order of the mixed FEM, but the effect was not obvious due to the limitation of computer hardware, so the convergence order could not be computed. In the future, we will continue to explore ways to improve performance.
Example 2. Consider the following Stokes equation
where \Omega = [0, 1]\times[0, 1] , the exact solution \mathbf{u} = (u_1, u_2)^t is
the exact solution p(x, y) is
and the body force \mathbf{f} = (f_1, f_2)^t is
where we set \nu = 1 .
Analogous to Example 1, mixed FEM with bivariate Bernstein polynomial basis are used to solve the above problems. The numerical errors in L^{\infty} , L^{2} , and H^{1} -semi norms are listed in Tables 5 and 6. The corresponding convergence orders are shown in Tables 7 and 8. Figures 3 and 4 display the error image.
It can be observed from the error line in Figures 3 and 4, and the error convergence order in Tables 7 and 8 that when the mesh size is equal, the higher the degree of Bernstein polynomial basis, not only the higher the numerical accuracy of the error norm, but also the higher the error convergence order.
Example 3. Consider the following non-homogenous 2D Stokes equation
where \Omega = [0, 1]\times[0, 1] , the exact solution \mathbf{u} = (u_1, u_2)^t is
the exact solution p(x, y) is
and the body force \mathbf{f} = (f_1, f_2)^t is
where we set \nu = 1 .
We continue to use the above method to solve Stokes Eq (4.3). The numerical errors are shown in Tables 9 and 10.
As can be seen from Tables 9 and 10, the Bernstein basis function shows good convergence under the three norms as the grid size decreases, which further verifies its advantages in numerical stability. In particular, the cubic or higher-degree Lagrange interpolation shows unstable oscillations, while the Bernstein basis function can provide a stable and consistent solution, while maintaining good geometric properties and flexible boundary condition processing capabilities. The numerical stability and global approximation characteristics of Bernstein polynomial make the results more reliable than Lagrange interpolation.
The Stokes equations are primarily used to describe fluid flow phenomena at low Re, where the inertial forces are significantly smaller compared to the viscous forces and can thus be neglected. This results in a flow that is smooth and orderly. Through the three numerical experiments presented above, we observe that as the mesh size decreases, the errors also gradually diminish. This indicates that our numerical solutions are progressively approaching the true laminar flow state. This trend demonstrates the effectiveness and accuracy of our numerical method in handling low Re fluid dynamics problems.
In this study, we use tensor product Bernstein polynomial basis function and Lagrange basis function to solve Stokes equation and verify the basis functions of different orders in detail. The results show that the solutions obtained by using bicubic or low-order Lagrange basis functions are basically the same as those obtained by using Bernstein polynomial basis functions in numerical accuracy and convergence order, with slight differences only after the decimal point of some p values, but the performance of Bernstein basis functions is slightly better overall. This shows that the performance of the two basis functions is equivalent in the case of lower order, but Bernstein basis function shows better stability in detail processing.
However, when the biquartic Lagrange basis function is used to solve the problem, the situation has changes significantly. In Example 1, the error result of biquartic Lagrange basis function is not as good as that of biquartic Lagrange basis function. In Examples 2 and 3, the solution of biquartic Lagrange basis function appears an obvious oscillation phenomenon, which leads to unreliable numerical results. This phenomenon shows that with the increase of the order of the basis function, and with Lagrange basis function being prone to numerical instability when dealing with complex problems, especially in the case of high order, this instability will be aggravated. In contrast, Bernstein basis functions show excellent numerical stability and higher accuracy in high-order cases. By using high-order Bernstein polynomial basis, we not only effectively alleviate the oscillation problem caused by high-order Lagrange basis function, but also generate more stable and accurate numerical solutions. Therefore, in order to ensure the reliability and stability of numerical results, we only show the solution results using Bernstein polynomial basis functions.
5.
Conclusions
We review the Bernstein polynomial basis and use it to construct the mixed finite element function space. Then, the Galerkin mixed FEM based on the bivariate Bernstein polynomial basis is used to solve the 2D Stokes equation, and the L^{\infty} , L^{2} , and H^{1} -semi norms of the error and convergence order between the exact solution and the finite element solution are calculated. At the same time, compared with the Lagrange basis function, the numerical accuracy and convergence order of solving Stokes equation with bicubic and below Lagrange interpolation polynomial basis and Bernstein polynomial basis are the same. High-order Lagrange interpolation function is often limited by Runge's phenomenon, so we use high-order Bernstein polynomial basis to effectively overcome this problem and obtain significantly better numerical results.
Author contributions
Lanyin Sun: Writing-review & editing, methodology, funding acquisition, conceptualization, visualization, data curation; Siya Wen: Writing-review & editing, writing-original draft, software. All authors have read and approved the final version of the manuscript for publication.
Use of Generative-AI tools declaration
The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.
Acknowledgments
The authors are very grateful to the editor and anonymous referees for their valuable comments and constructive suggestions, which helped to improve the paper significantly. This work is partly supported by Program for Science Technology Innovation Talents in Universities of Henan Province (No. 22HASTIT021), the Science and Technology Project of Henan Province (No. 212102210394).
Conflict of interest
The authors declare that there is no conflict of interest.
DownLoad:
