Consistent vs. reduced integration penalty methods for incompressible media using several old and new elements

The frequently used reduced integration method for solving incompressible flow problems ‘a la penalty’ is critically examined vis-a-vis the consistent penalty method. For the limited number of quadrilateral and hexahedral elements studied, it is shown that the former method is only equivalent to the latter in certain special cases. In the general case, the consistent penalty method is shown to be more accurate. Finally, we demonstrate significant advantages of a new element, employing biquadratic (2-D) or triquadratic (3-D) velocity and linear pressure over that using the same velocity but employing bilinear (2-D) or trilinear (3-D) pressure approximation.     

Spectral/hp penalty least‐squares finite element formulation for unsteady incompressible flows

In this paper, we present spectral/hp penalty least‐squares finite element formulation for the numerical solution of unsteady incompressible Navier–Stokes equations. Pressure is eliminated from Navier–Stokes equations using penalty method, and finite element model is developed in terms of velocity, vorticity and dilatation. High‐order element expansions are used to construct discrete form. Unlike other penalty finite element formulations, equal‐order Gauss integration is used for both viscous and penalty terms of the coefficient matrix. For time integration, space–time decoupled schemes are implemented. Second‐order accuracy of the time integration scheme is established using the method of manufactured solution. Numerical results are presented for impulsively started lid‐driven cavity flow at Reynolds number of 5000 and transient flow over a backward‐facing step. The effect of penalty parameter on the accuracy is investigated thoroughly in this paper and results are presented for a range of penalty parameter. Present formulation produces very accurate results for even very low penalty parameters (10–50). Copyright © 2008 John Wiley & Sons, Ltd.     

Variable penalty method for finite element analysis of incompressible flow

A new scheme is applied for increasing the accuracy of the penalty finite element method for incompressible flow by systematically varying from element to element the sign and magnitude of the penalty parameter λ, which enters through ?.v + p/λ = 0, an approximation to the incompressibility constraint. Not only is the error in this approximation reduced beyond that achievable with a constant λ, but also digital truncation error is lowered when it is aggravated by large variations in element size, a critical problem when the discretization must resolve thin boundary layers. The magnitude of the penalty parameter can be chosen smaller than when λ is constant, which also reduces digital truncation error; hence a shorter word-length computer is more likely to succeed. Error estimates of the method are reviewed. Boundary conditions which circumvent the hazards of aphysical pressure modes are catalogued for the finite element basis set chosen here. In order to compare performance, the variable penalty method is pitted against the conventional penalty method with constant λ in several Stokes flow case studies.     

A least-squares finite element method for time-dependent incompressible flows with thermal convection

The time-dependent Navier–Stokes equations and the energy balance equation for an incompressible, constant property fluid in the Boussinesq approximation are solved by a least-squares finite element method based on a velocity–pressure–vorticity–temperature–heat-flux ( u –P–ω–T– q ) formulation discretized by backward finite differencing in time. The discretization scheme leads to the minimization of the residual in the l²-norm for each time step. Isoparametric bilinear quadrilateral elements and reduced integration are employed. Three examples, thermally driven cavity flow at Rayleigh numbers up to 10⁶, lid-driven cavity flow at Reynolds numbers up to 10⁴ and flow over a square obstacle at Reynolds number 200, are presented to validate the method.     

Superconvergence of nonconforming finite element penalty scheme for Stokes problem using L 2 projection method

A modified penalty scheme is discussed for solving the Stokes problem with the Crouzeix-Raviart type nonconforming linear triangular finite element. By the L ² projection method, the superconvergence results for the velocity and pressure are obtained with a penalty parameter larger than that of the classical penalty scheme. The numerical experiments are carried out to confirm the theoretical results.     

A new general-purpose least-squares finite element model for steady incompressible low-viscosity laminar flow using isoparametric C1-continuous Hermite elements

This paper deals with a critical evaluation of various finite element models for low-viscosity laminar incompressible flow in geometrically complex domains. These models use Galerkin weighted residuals UVP, continuous penalty, discrete penalty and least-squares procedures. The model evaluations are based on the use of appropriate tensor product Lagrange and simplex quadratic triangular elements and a newly developed isoparametric Hermite element. All of the described models produce very accurate results for horizontal flows. In vertical flow domains, however, two different cases can be recognized. Downward flows, i.e. when the gravitational force is in the direction of the flow, usually do not present any special problem. In contrast, laminar flow of low-viscosity Newtonian fluids where the gravitational force is acting in the direction opposite to the flow presents a difficult case. We show that only by using the least-squares method in conjunction with C¹-continuous Hermite elements can this type of laminar flow be modelled accurately. The problem of smooth isoparametric mapping of C¹ Hermite elements, which is necessary in dealing with geometrically complicated domains, is tackled by means of an auxiliary optimization procedure. We conclude that the least-squares method in combination with isoparmetric Hermite elements offers a new general-purpose modelling technique which can accurately simulate all types of low-viscosity incompressible laminar flow in complex domains.     

Efficient quadrature for a boundary element method to compute three-dimensional Stokes flow

A collocation-type boundary element method based on bilinear B-splines is used for the numerical solution of the Stokes Dirichlet problem in bounded domains D ? R³. The computation of the influence matrix requires the numerical evaluation of weakly singular integrals on the domain boundary if the usual double-layer potential ansatz is chosen. Here mostly standard methods with disjoint grids for collocation and integration are used. We develop a special integration scheme based on triangular co-ordinates near the singularity and show its efficiency compared with the method mentioned above.     

Two-level stabilized finite element method for Stokes eigenvalue problem

A two-level stabilized finite element method for the Stokes eigenvalue problem based on the local Gauss integration is considered.This method involves solving a Stokes eigenvalue problem on a coarse mesh with mesh size H and a Stokes problem on a fine mesh with mesh size h = O(H 2),which can still maintain the asymptotically optimal accuracy.It provides an approximate solution with the convergence rate of the same order as the usual stabilized finite element solution,which involves solving a Stokes eigenvalue problem on a fine mesh with mesh size h.Hence,the two-level stabilized finite element method can save a large amount of computational time.Moreover,numerical tests confirm the theoretical results of the present method.     

An extension of the SIMPLE based discontinuous Galerkin solver to unsteady incompressible flows

In this paper, we present a SIMPLE based algorithm in the context of the discontinuous Galerkin method for unsteady incompressible flows. Time discretization is done fully implicit using backward differentiation formulae (BDF) of varying order from 1 to 4. We show that the original equation for the pressure correction can be modified by using an equivalent operator stemming from the symmetric interior penalty (SIP) method leading to a reduced stencil size. To assess the accuracy as well as the stability and the performance of the scheme, three different test cases are carried out: the Taylor vortex flow, the Orr‐Sommerfeld stability problem for plane Poiseuille flow and the flow past a square cylinder. (1) Simulating the Taylor vortex flow, we verify the temporal accuracy for the different BDF schemes. Using the mixed‐order formulation, a spatial convergence study yields convergence rates of k + 1 and k in the L²‐norm for velocity and pressure, respectively. For the equal‐order formulation, we obtain approximately the same convergence rates, while the absolute error is smaller. (2) The stability of our method is examined by simulating the Orr–Sommerfeld stability problem. Using the mixed‐order formulation and adjusting the penalty parameter of the symmetric interior penalty method for the discretization of the viscous part, we can demonstrate the long‐term stability of the algorithm. Using pressure stabilization the equal‐order formulation is stable without changing the penalty parameter. (3) Finally, the results for the flow past a square cylinder show excellent agreement with numerical reference solutions as well as experiments. Copyright © 2015 John Wiley & Sons, Ltd.     

A sensitivity analysis on the parameter of the GLS method for a second‐gradient theory of incompressible flow

Using a non‐conforming C⁰‐interior penalty method and the Galerkin least‐square approach, we develop a continuous–discontinuous Galerkin finite element method for discretizing fourth‐order incompressible flow problems. The formulation is weakly coercive for spaces that fail to satisfy the inf‐sup condition and consider discontinuous basis functions for the pressure field. We consider the results of a stability analysis through a lemma which indicates that there exists an optimal or quasi‐optimal least‐square stability parameter that depends on the polynomial degree used to interpolate the velocity and pressure fields, and on the geometry of the finite element in the mesh. We provide several numerical experiments illustrating such dependence, as well as the robustness of the method to deal with arbitrary basis functions for velocity and pressure, and the ability to stabilize large pressure gradients. We believe the results provided in this paper contribute for establishing a paradigm for future studies of the parameter of the Galerkin least square method for second‐gradient theory of incompressible flow problems. Copyright © 2015 John Wiley & Sons, Ltd.     

Mixed hybrid penalty finite element method and its application

The penalty and hybrid methods are being much used in dealing with the general incompatible element, With the penalty method convergence can always be assured, but comparatively speaking its accuracy is lower, and the condition number and sparsity are not so good. With the hybrid method, convergence can be assured only when the rank condition is satisfied. So the construction of the element is extremely limited. This paper presents the mixed hybrid penalty element method, which combines the two methods together. And it is proved theoretically that this new method is convergent, and it has the same accuracy, condition number and sparsity as the compatible element. That is to say, they are optimal to each other.Finally, a new triangle element for plate bending with nine freedom degrees is constructed with this method (three degreesof freedom are given on each corner——one displacement and tworotations), the calculating formula of the element stiffness matrix is almost the same as that of the old triangle element for plate bending with nine degrees of freedom But it is converged to true solution with arbitrary irregrlar triangle subdivision. If the true solution u∈H³ with this method the linear and quadratic rates of convergence are obtianed for three bending moments and for the displacement and two rotations respectively.     

An algebraic factorisation scheme for spectral element solution of incompressible flow and scalar transport

A spectral element algorithm for solution of the unsteady incompressible Navier–Stokes and scalar (species/heat) transport equations is developed using the algebraic factorisation scheme. The new algorithm utilises Nth order Gauss–Lobatto–Legendre points for velocity and the scalar, while (N-2)th order Gauss–Legendre points are used for pressure. As a result, the algorithm does not require inter-element continuity for pressure and pressure boundary conditions on solid surfaces. Implementations of the algorithm are performed for conforming and non-conforming grids. The latter is accomplished using both the point-wise matching and integral projection methods, and applied for grids with both polynomial and geometric non-conformities. Code validation cases include the unsteady scalar convection equation, and Kovasznay flow in two- and three-dimensional domains. Using cases with analytical solutions, the algorithm is shown to achieve spectral accuracy in space and second-order accuracy in time. The results for the Boussinesq approximation for buoyancy-driven flows, and the species mixing in a continuous flow micro-mixer are also included as examples of applications that require long-time integration of the scalar transport equations.     

A numerical analysis of a class of generalized Boussinesq‐type equations using continuous/discontinuous FEM

An improved class of Boussinesq systems of an arbitrary order using a wave surface elevation and velocity potential formulation is derived. Dissipative effects and wave generation due to a time‐dependent varying seabed are included. Thus, high‐order source functions are considered. For the reduction of the system order and maintenance of some dispersive characteristics of the higher‐order models, an extra O(μ²ⁿ⁺²) term (n ∈ ?) is included in the velocity potential expansion. We introduce a nonlocal continuous/discontinuous Galerkin FEM with inner penalty terms to calculate the numerical solutions of the improved fourth‐order models. The discretization of the spatial variables is made using continuous P₂ Lagrange elements. A predictor‐corrector scheme with an initialization given by an explicit Runge–Kutta method is also used for the time‐variable integration. Moreover, a CFL‐type condition is deduced for the linear problem with a constant bathymetry. To demonstrate the applicability of the model, we considered several test cases. Improved stability is achieved. Copyright © 2011 John Wiley & Sons, Ltd.     

A Petrov–Galerkin finite element model for one‐dimensional fully non‐linear and weakly dispersive wave propagation

A new finite element method is presented to solve one‐dimensional depth‐integrated equations for fully non‐linear and weakly dispersive waves. For spatial integration, the Petrov–Galerkin weighted residual method is used. The weak forms of the governing equations are arranged in such a way that the shape functions can be piecewise linear, while the weighting functions are piecewise cubic with C²‐continuity. For the time integration an implicit predictor–corrector iterative scheme is employed. Within the framework of linear theory, the accuracy of the scheme is discussed by considering the truncation error at a node. The leading truncation error is fourth‐order in terms of element size. Numerical stability of the scheme is also investigated. If the Courant number is less than 0.5, the scheme is unconditionally stable. By increasing the number of iterations and/or decreasing the element size, the stability characteristics are improved significantly. Both Dirichlet boundary condition (for incident waves) and Neumann boundary condition (for a reflecting wall) are implemented. Several examples are presented to demonstrate the range of applicabilities and the accuracy of the model. Copyright © 2001 John Wiley & Sons, Ltd.     

An exact element method for plane problem

In this paper, based on the step reduction method^[1] and exact analytic method^[2] anew method-exact element method for constructing finite element, is presented. Since the new method doesn ’t need the variational principle, it can be applied to solve non-positive and positive definite partial differential equations with arbitrary variable coefficient. By this method, a quadrilateral noncompatible element with 8 degrees of freedom is derived for the solution of plane problem. Since Jacobi ’s transformation is not applied, the present element may degenerate into a triangle element. It is convenient to use the element in engineering. In this paper, the convergence is proved. Numerical examples are given at the endof this paper, which indicate satisfactory results of stress and displacements can be obtained and have higher numerical precision in nodes.     

Nonstationary pandom vibpation analysis of linear elastic structures with finite element method

At present, the finite element method is an efficient method for analyzing structural dynamic problems. When the physical quantities such as displacements and stresses are resolved in the spectra and the dynamic matrices are obtained in spectral resolving form, the relative equations cannot be solved by the vibration mode resolving method as usual. For solving such problems, a general method is put forward in this paper. The excitations considered with respect to nonstationary processes are as follows: P(t)={P_i(t)},P_i(t)=a_i(t)P_i(t), a_i(t) is a time function already known. We make Fourier transformation for the discretized equations obtained by finite element method, and by utilizing the behaviour of orthogonal increment of spectral quantities in random process^[1], some formulas of relations about the spectra of excitation and response are derived. The cross power spectral denisty matrices of responses can be found by these formulas, then the structrual safety analysis can be made. When a_i(t)=l (i= 1,2,…n), the. method stated in this paper will be reduced to that which is used in the special case of stationary process.     

Least-squares finite element methods for compressible Euler equations

A method based on backward finite differencing in time and a least-squares finite element scheme for first-order systems of partial differential equations in space is applied to the Euler equations for gas dynamics. The scheme minimizes the L²-norm of the residual within each time step. The method naturally generates numerical dissipation proportional to the time step size. An implicit method employing linear elements has been implemented and proves robust. For high-order elements, computed solutions based on the L²-method may have oscillations for calculations at similar time step sizes. To overcome this difficulty, a scheme which minimizes the weighted H¹-norm of the residual is proposed and leads to a successful scheme with high-degree elements. Finally, a conservative least-squares finite element method is also developed. Numerical results for two-dimensional problems are given to demonstrate the shock resolution of the methods and compare different approaches.     

Reduced augmented Lagrangian bi-conjugate gradient method for impact-contact problems

To avoid the numerical oscillation of the penalty method and non-compatibility with explicit operators of conventional Lagrange multiplier methods used in transient contact problems to enforce surface contact conditions, a new approach to enforcing surface contact constraints for the transient nonlinear finite element problems, referred to as “the reduced augmented Lagrangian bi-conjugate gradient method (ALCG)”, is developed in this paper. Based on the nonlinear constrained optimization theory and is compatible with the explicit time integration scheme, this approach can also be used in implicit scheme naturally. The new surface contact constraint method presented has significant advantages over the widely adopted penalty function methods and the conventional Lagrangian multiplier methods. The surface contact constraints are satisfied more accurately for each step by the algorithm, so the oscillation of numerical solution for the explicit scheme is depressed. Through the development of new iteration strategy for solving nonlinear equations, ALCG method improves the computational efficiency greatly. Project supported by State Education Commission Doctoral Foundation and Natural Science Foundation of Liaoning Province.     

On the difficulties and remedies in enforcing the div = 0 condition in the finite element analysis of thermal plumes with strongly temperature-dependent viscosity

A finite element (FE) analysis of experimentally observed creeping thermal plumes in a medium whose viscosity is strongly temperature-dependent is performed. Such plumes are considered to play an important role in numerous geological processes and numerical modelling is often the only option to study their physics. Initial simulations by means of the general-purpose Galerkin finite element package NACHOS-II demonstrated serious deficiencies of the method in modelling plumes with large viscosity contrasts, in spite of several options for the solution (mixed or penalty formulation) and the elements (continuous or discontinous pressure). In agreement with observations from FE simulations of isothermal Stokes flow in other studies, we have isolated the violation of the div = 0 or incompressibility constraint as the major culprit in the failure of the FE method. It is demonstrated that the a posteriori computed discrete divergence (DDIV) can be used as a diagnostic tool to evaluate the reliability of the FE solution and to rank the solution and element options provided in the NACHOS code. On the basis of these considerations, the combination of the mixed method with a Q₂-P₁ (discontinuous pressure) element turns out to be the most suitable for the present plume problem but is still unable to sufficiently enforce the div = 0 condition. With a goal to remedy this detrimental behaviour, several FE modifications and new approaches have been taken. These include: (i) use of a new scaling option for the governing equations which has the effect of equilbrating the stiffness matrices and thus improving their condition; (ii) implementation of several iterative solution techniques such as the iterated penalty and the Uzawa algorithm for the augmented Langrangian to better accommodate the dual role of the pressure; (iii) use of a multistep Newton method to better handle the high non-linearity of the coupled flow/transport problem. Although each of these options (or a combination of them) is able to improve on the quality of FE solution, the most startling amelioration has been gained with option (iii). Use of the latter resulted in very satisfactory modelling of the experimentally observed plumes.     

The application of compatible stress iterative method in dynamic finite element analysis of high velocity impact

There is a common difficulty in elastic-plastic impact codes such as EPIC ^[2,3].NONSAP ^[4],etc. Most of these codes use the simple linear functions usually taken from static problem to represent the displacement components. In such finite element formulation, the stress components are constant in each element and they are discontinuous in any two neighboring elements. Therefore, the bases of using the virtual work principle in such elements are unreliable. In this paper, we introduce a new method, namely, the compatible stress iterative method, to eliminate the above-said difficulty. The calculated examples show that the calculation using the new method in dynamic finite element analysis of high velocity impact is valid and stable, and the element stiffness can be somewhat reduced.This is a part of the author's Ph. D dissertation under the supervision of Professor Chien Wei-zang.