3.4.1 Literature review
The Galerkin finite element model has been proven eminently successful in problems of solid/structural mechanics and in other situations such as the heat conduction which is governed by diffusion-type equation. However, the computational difficulty is perceived by the Galerkin finite element model in fluid mechanics. The solutions of convection
these oscillations become even worse when high Peclet number problem is considered.
The above difficulty was blamed on the central difference approximation generated by the Galerkin model and its inability to account for the upwind effect of the convection term. The descretization of the convection term by Galerkin model produces a set of equations that involve the adjacent nodes, and it results in node-to-node oscillations. That is, the best approximation property in the energy norm of the Galerkin model, which is the basis for success in symmetric operator, is lost when convection dominates the diffusion [14].
In principle, these oscillations can be suppressed by a successive refinement of the mesh to make the mesh Peclet number to be less than 2. However, the necessary degree of the mesh refinement is always economically impracticable. This computational difficulty and perceived shortcomings of the Galerkin model have motivated the development of the so-called Petrov-Galerkin models, in which the trial function and weighting function are chosen from the different functional spaces.
Inspired by the upwind operator in finite difference context [5], Christie [30] and Hein-rich [31, 32] developed the early upwind finite element models for the steady-state one-and two-dimensional convection-diffusion equation, respectively. The weighting func-tions they proposed are equal to the weighting function of Galerkin model plus a high order polynomial in order to place more weight on the upwind side. Later, Hughes [33], in a different manner, modified the numerical quadrature rule for the convection term to achieve the upwind effect. Unfortunately, the success of these early Petrov-Galerkin models was only limited to one-dimensional problem because they produced in general excessive crosswind diffusion which reduces the solution accuracy in multi-dimensional problem. This problem is particularly significant in cases when mesh lines and flow di-rections are not in good alignment. In addition, these models are only first-order accurate and its extension to the unsteady cases or for the case of non-zero source term would not give a consistent formulation.
3.4.2 Streamline Upwind Petrov-Galerkin model
The drawbacks of the early Petrov-Galerkin models [30–33] revealed that the flow-oriented convection discretization scheme should be a better strategy to enhance convective sta-bility without sacrificing solution accuracy. In the beginning of the 1980s, Hughes and Brook pointed out that to enhance convective stability and reduce the crosswind diffusion, it is necessary to add a proper amount of artificial damping along the direction of primary flow. The introduced stabilization term was added to the Galerkin weighting function and it is called the streamline upwind Petrov-Galerkin (SUPG) model [34]. It is also known as a perturbation method since the added stabilization term (or called the artificial damping) can be viewed as a perturbation to the Galerkin weighting function. That is, the weighting function of SUPG model is given below :
wSUPG = wGalerkin+τSUPGBi (3.10)
Bi= u· ∇ϕ (3.11)
where Bidenotes the biased part. Later, Mizukami and Hughes extended the application range of SUPG to flow problem containing a sharp layer by demanding that the discrete system underlying the streamline upwind scheme ensures the satisfying of the maximum principle [35]. In the presence of a boundary and an internal sharp layers, the predicted SUPG solution quality was deteriorated further. To overcome this shortcoming, Hughes and his colleagues improved the SUPG model by adding the discontinuity-capturing term to the SUPG weighting function so as to produce a smooth solution in the boundary or interior layer problem [36].
In Eq. (3.10), the stabilization parameter τSUPG has the significant effect on the sup-pression of the oscillations and the solution accuracy. The parameter τSUPG in [37] is defined as
τSUPG =δSUPGuH 2|u|2
where H andδ represent the element characteristic length (H = 2h) and the upwind
coefficient, respectively. In this study, the quadratic element is adopted in two- and three-dimensional problem. The upwind coefficientδSUPG will be, therefore, determined at the center and the corner node, respectively. One way to determine the coefficient δSUPG is to impose one-dimensional exactness of the homogeneous and steady-state convection diffusion equation given below [37].
δSUPG(γ) =
4 tanh(γ/2) − sinh(γ) − (6/γ)sinh(γ)tanh(γ/2) at corner node 1
2coth(γ 2)−1
γ at center node
(3.12)
whereγ = |u|H2ε represents the element Peclet number.
In this study, in the different manner, the upwind coefficient will be derived by mini-mizing the wavenumber error of the convection term in next section.
3.4.3 Modified wavenumber error optimizing model
Since the numerical instability is generated from the Galerkine finite element model treat-ment of the convection term, one way to enhance convective stability is to take the dis-persive nature of the investigated convection (or first-order derivative) term into consider-ation.
A scheme for approximating the convection term can be regarded to be optimized if the error between the exact and numerical wavenumbers is minimized. This can be achieved by applying the Fourier transform ˜ϕ(k) = 21π∫−∞∞ ϕ(x)exp(−ikx)dx and its in-verseϕ(x) =∫−∞∞ ϕ(k)exp(ikx)dk, where i =˜ √
−1, to the convection term. In this study, the stabilization parameter in modified wavenumber error (MWE) modelτMWE is chosen as
τMWE=δmMWEH
2|u| (3.13)
Following the work of Tam and Webb [38], derivation of the upwind coefficientδmMWEis to minimize the wavenumber error of the convection term. As mentioned before, the nodal point in a quadratic element can be divided into the center and corner nodes, respectively.
The superscript m in Eq. (3.13) isα or β depending on the nodal classification within a quadratic element.
To minimize the wavenumber error, we take the convection termϕxin one-dimensional equation as an example to determineτMWE. Substituting the trial function approximation ϕ = ∑3j=1ϕjNj and the weighting function wi = Ni+τu∂N∂xi, where Ni(i = 1∼ 3) rep-resent the quadratic interpolation functions, into the weighted residuals statement forϕx, the discretized equation at the center and corner nodes can be derived as follows :
ϕx|center≈ 1 its inverse transformation on each term shown in Eqs. (3.14), the numerical wavenumber k at the center node can be derived as k≈−ih(a1exp(−ikh)+a2+ a3exp(ikh)). The exact wavenumber ˜k is regarded as the right hand side of the above relation, thereby leading to
˜k = −i
h (a1exp(−ikh) + a2+ a3exp(ikh)).
To demand k be close to ˜k, the integral quantity E(k) defined below should be a very small and positive value
E(k) =
∫ π/2
−π/2|kh − ˜kh|2d(kh)
To get the smallest value of E from the above equation, the limiting condition∂E/∂a1is applied to get the coefficient δαMWE. The coefficientδβMWE can be derived similarly. The derived upwinding coefficients are summarized below
δMWE=
Comparing the upwind coefficients in Eq. (3.12) and Eq. (3.16), the main difference is
that the proposed stabilization parameter is a constant and the computational cost can be therefore reduced when evaluating the stabilization parameterτMWE. Another difference is that the parameterτMWEis independent of the convection term. Such a difference becomes significant when solving the nonlinear incompressible Navier-Stokes equations, where the convection velocity field will be updated in each nonlinear iteration. For the sake of clarify, the weighting functions of different finite element model are plotted in Fig. 3.1. It is clearly seen from Fig. 3.1 that the proposed MWE model gives more weight to upstream nodes than the Galerkin and SUPG model.
The extension of the above Petrov-Galerkin formulations to multi-dimensional prob-lem is, however, obtained at the cost of accuracy. The reason is due to the use of con-vectional Petrov-Galerkin model which has a tendency to add a false diffusion to regions normal to the streamline as the angle between the grid line and the flow direction is large.
To improve this situation, the weighting function should accommodate a streamline oper-ator so as to provide a natural mechanism for adding the artificial damping to the weight-ing function along the direction of primary flow. In a quadratic element, the stabilization parameters are assigned at the center and corner nodes to stabilize the discrete system [39]
τMWE= ∑ni=1δiMWEHi
2|u| , (n = 2 or 3) (3.17)
In the above,δMWE is obtained by Eq. (3.16)