-
When the horizontal velocity v = (u, v) is represented by (ψ, χ), we have v = k × ∇ψ + ∇χ, where k is the unit vector in the vertical direction and ∇ is the gradient operator in the two-dimensional space of x ≡ (x, y). As mentioned in the introduction section, (ψ, χ) contains no harmonic component in an unbounded or periodic domain, so (ψ, χ) can be uniquely defined by uniquely relating k × ∇ψ and ∇χ to the rotational and divergent parts of v, respectively. In a limited domain, (ψ, χ) contain a harmonic component and thus cannot be uniquely defined and determined due to the nonunique partition and attribution of the harmonic component to (ψ, χ). In this case, the following relationships hold:
where ζ = k • ∇×v =
$\partial $ xv –$\partial $ yu is the vorticity, α = ∇ • v =$\partial $ xu +$\partial $ yv is the divergence, D denotes the limited domain, S denotes the boundary of D,$\partial $ n (or$\partial $ s) denotes the boundary-normal (or boundary-tangential) component of ∇, and vn (or vs) is the boundary-normal (or boundary-tangential) component of v on S.Applying the Gauss and Stokes's theorems to Eqs. (1) and (2) leads to the following solvability conditions:
where ∫S ( )dl denotes the line-integration of ( ) along the closed loop of S and ∫D ( )dx denotes the area-integration of ( ) over D. These two solvability conditions are accurately satisfied in discrete forms by the original v fields used in this paper in three mesoscale domains (as shown later in Table 2) but not accurately by the reconstructed velocity fields (as shown later in Table 3).
Domain Method SCC RRD D1 SOR-based 0.9392 0.2820 M-spectral 0.7836 0.5479 Integral 0.9870 0.1156 Integral-spectral 0.9867 0.1151 Integral-SOR 0.9958 0.0720 D2 SOR-based 0.8644 0.2641 M-spectral 0.8661 0.3680 Integral 0.9809 0.1115 Integral-spectral 0.9822 0.1032 Integral-SOR 0.9932 0.0648 D3 SOR-based 0.8889 0.5024 M-spectral 0.8753 0.5786 Integral 0.9740 0.2240 Integral-spectral 0.9633 0.2698 Integral-SOR 0.9909 0.1288 Table 1. The SCC and RRD value between the original v and reconstructed vc from each listed method for each domain. See Eqs. (9a, b) for definitions of SCC and RRD. Bold fonts indicate those values where both SCC and RRD satisfy the accuracy adequacy criterion given by Eq. (10).
Domain ∫S vsds ∫D ζdx ∫S vndl ∫Dαdx D1 –10.5308 –10.5308 –10.8613 –10.8613 D2 2.6230 2.6230 –6.7034 –6.7034 D3 –1.9023 –1.9023 –13.8323 –13.8323 Table 2. Values of the left-hand-side and right-hand-side terms of solvability conditions in Eqs. (3a) and (3b) computed by using the original v in each domain are listed in each row (units: 10–5 m2 s–1) .
Domain Method ∫Svsdl ∫Svndl D1 SOR-based 0.9852 0.5884 M-spectral 1.7570 1.0915 Integral 1.1236 1.0047 Integral-SOR 0.9973 1.0047 Integral-spectral 0.9868 1.0757 D2 SOR-based 0.8368 0.4332 M-spectral 0.5385 1.4400 Integral 1.7196 0.8844 Integral-SOR 0.9247 0.8844 Integral-spectral 0.8302 0.9090 D3 SOR-based 0.7653 0.4765 M-spectral 0.7808 1.1484 Integral 0.6659 0.9457 Integral-SOR 0.9812 0.9457 Integral-spectral 0.8176 0.9693 Table 3. Normalized values of the left-hand-side terms of solvability conditions in Eqs. (3a) and (3b) computed with (vs, vn) given by vc in each domain from each method are listed in each row. The values of ∫Svsdl from the five methods listed for each domain are normalized by the value of ∫Svsdl listed for the same domain in Table 2. The values of ∫Svsdl from the five methods listed for each domain are normalized by the value of ∫Svndl listed for the same domain in Table 2.
The above two solvability conditions can be expressed in terms of (ψ, χ). In this case, each solvability condition can be split into two parts as shown below. Since (ψ, χ) are single-valued continuous functions along the closed loop of S, the first part of the solvability condition in Eq. (3a) or (3b) in terms of (ψ, χ) satisfies:
The remaining second part of the solvability condition can be then obtained by substituting Eqs. (2a) and (4a) [or (2b) and (4b)] into Eq. (3a) [or (3b)], which gives
-
In the integral method, the internally induced (ψ, χ), denoted by (ψin, χin), is given by
where r = |x' − x| and (2
$\pi $ )–1 lnr is the free-space Greenꞌs function for Poisson equation in an unbounded domain but here it is applied to ($\zeta $ ,$\alpha $ ) in D. Three discretization schemes were developed in CX11 for computing the integral in Eq. (6), and the first scheme that used a staggered grid for ($\zeta $ ,$\alpha $ ) is used in this paper. Ideally, (ψin, χin) should exactly satisfy those solvability conditions that are further partitioned from Eqs. (4) and (5), that is,Practically, the computed (ψin, χin) cannot accurately satisfy the second solvability condition in Eq. (7b) due mainly to discretization errors in computing ∂n(ψin, χin), although the first solvability condition in Eq. (7a) is exactly satisfied. The errors involved are usually small and negligible. Still, they can increase significantly as the original v becomes strong with sharp variations across or along S (as mentioned in the introduction and seen later from Tables 3 and 4). Note that Green's function in Eq. (6) involves no boundary condition, and the derived integral form of (ψin, χin) is simple and unaffected by any complications in the boundary condition. Owing to the simplicity of this integral form and the use of a staggered grid in discretizing (ζ, α), the computed (ψin, χin) is more accurate and less sensitive to increased and sharpened spatial variations of original v in D than the solutions computed by the other two methods. This is an important strength of the integral method that should be adopted for the purpose of first obtaining (ψin, χin) in a hybrid approach.
Domain ∫Svs,indl ∫Svn,indl D1 0.9725 1.0047 D2 0.9421 0.8844 D3 0.7188 0.9457 Table 4. As in Table 3, but computed with (vs,in, vn,in) given by the internally-induced component of vc in each domain obtained from the first step of the integral method listed in each row.
The remaining externally induced harmonic component of (ψ, χ), denoted by (ψe, χe), is computed in the second step by simply setting χe = 0. With χe = 0, (ψ, χ) can be uniquely defined by relating k × ∇ψ (or ∇χ) to the nondivergent (or divergent) part of v, where the nondivergent part consists of the rotational part and harmonic part of v. With the above simplification, four discretization schemes were developed in CX11 for computing ψe and the Cauchy-integral scheme (shown in section 4 of CX11) is used in this paper. This scheme computes the complex velocity potential defined by ωe = φe – iψe where i is the imaginary unit, while the boundary value of ψe (or φe) is obtained by integrating
$\partial $ sψe = vn,in – vn (or$\partial $ sφe = vs – vs,in) along S as shown in (3.9) [or (3.10)] of Xu et al. (2011), where vn,in (or vs,in) is the boundary-normal (or boundary-tangential) component of vin ≡ k × ∇ψin + ∇χin on S. Since ∇2(ψe, φe) = (0, 0) and ∂s(ψe, φe) =$\partial$ n(–φe, ψe) on S according to the Cauchy-Riemann conditions, the paired solvability conditions for (ψe, φe) can be derived from the partitioned solvability conditions in Eqs. (7a, b) in the following form:Note from Eqs. (2a, b) that vs,in =
$\partial $ nψin +$\partial $ sχin and vn,in =$\partial $ nχin –$\partial $ sψin on S, so the solvability conditions in Eq. (8) can be exactly satisfied only if (ψin, χin) exactly satisfies the solvability conditions in Eqs. (7a, b). As explained earlier, for (ψin, χin) computed in the first step, the paired solvability conditions in Eq. (7a) are satisfied exactly, but the paired solvability conditions in Eq. (7b) are not; therefore, the paired solvability conditions in Eq. (8) are also not exactly satisfied. Consequently, the boundary value of ψe (or φe) obtained by integrating ∂nχin – vn (or vs – ∂nψin) along the boundary loop undergoes a discontinuous jump when the integration returns to the beginning point. Although such a jump can be spread out and somewhat diluted by redistributing it evenly over the entire boundary loop, the paired solvability conditions in Eq. (7b) are still not exactly satisfied, and their induced errors in the computed boundary values of (ψe, φe) remain essentially intact. While more accurate than the other three schemes developed in CX11 for computing ψe (with χe = 0), the Cauchy-integral scheme is quite sensitive to small errors in the computed boundary values of (ψe, χe) that cannot accurately satisfy the paired solvability conditions in Eq. (8), because the scheme is derived from an analytical integral and this integral requires the paired solvability conditions in Eq. (8) to be precisely satisfied. This appears to be a weakness for the integral method. Still, it can be avoided if ψe (with χe = 0) can be computed by another method with an improved accuracy after (ψin, χin) is obtained by the integral method via a hybrid approach. -
The SOR-based method (Sangster, 1960) solves the coupled system of Eqs. (1)–(2) with either ψ or χ freely prescribed on S. For example, if ψ is chosen to be freely prescribed on S, then ψ can be solved first from Eq. (1a) in D, which gives
$\partial{\boldsymbol{}} $ nψ along S. Thus,$\partial $ sχ can be obtained along S from Eq. (2a), and integrating$\partial $ sχ along S gives the boundary value of χ that allows χ be finally solved from Eq. (1b) in D. In this case, the solvability condition in Eq. (4b) is satisfied exactly, but the solvability condition in Eq. (5a) is satisfied only approximately due to discretization errors in solving for ψ and computing$\partial $ nψ. Thus, when the boundary value of$\partial $ sχ is obtained from Eq. (2a), its boundary loop integral is not exactly zero. In other words, the solvability condition in Eq. (4a) is not exactly satisfied, so the boundary value of χ obtained by integrating$\partial $ sχ along S undergoes a discontinuous variation when the loop integration returns to the beginning point.Similarly but alternatively, if χ is chosen to be freely prescribed on S, with χ solved first from Eq. (1b) in D, then the solvability condition in Eq. (4a) is satisfied exactly, but the solvability condition in Eq. (5b) is satisfied only approximately due to discretization errors in solving for χ and computing
$\partial $ nχ. Thus, when the boundary value of$\partial $ sψ is obtained from Eq. (2b), its boundary loop integration is not exactly zero. In this case, the solvability condition in Eq. (4b) is not exactly satisfied and the boundary value of ψ, obtained by integrating$\partial $ sψ along S, undergoes a discontinuous jump when the loop integration returns to the beginning point.The abovementioned jump can be negligibly small if the original v is sufficiently smooth. Still, it can become large if the original v demonstrates sharp variations not adequately resolved by the analysis grid, especially across or along the domain boundary (as shown later in Table 3). Although this jump can be thinned out by redistributing it evenly over the entire boundary loop, its implied error in the boundary condition remains essentially intact.
In this paper, the SOR is implemented with (ψ, χ) and (ζ, α) discretized on the same grid. The discretization errors of ∇2(ψ, χ) and (ζ, α) are large in areas of complex flow, and these large errors tend to accumulate as the iterative process of SOR goes sequentially point-by-point through the areas of complex flow. The SOR-based method is thus relatively sensitive to (and its accuracy can be significantly reduced by) increased spatial variations of original v and thus further sharpened variations of (ζ, α) that were computed from original v inside the domain D (as shown later in section 4). On the other hand, if the SOR-based method is used to solve for (ψe, χe) with a vanished (ζ, α) inside D, then the discretization errors of ∇2(ψ, χ) and (ζ, α) will become small or much smaller everywhere inside D, so (ψe, χe) can be solved with improved accuracy. This implies that the SOR-based method should have an advantage if it is used to solve for the externally induced (ψe, χe) after the internally induced (ψin, χin) is obtained by the integral method via a hybrid approach. Such a hybrid approach will be considered in section 2.4.
-
As mentioned in the introduction, the spectral method of CK92a,b has two versions, called the S-version and C-version, respectively. In either version, (ψ, χ) is divided into two parts, called the inner part and harmonic part and denoted by (ψin,n, χin,n) and (ψh,a, χh,a), respectively. The inner part is well defined in the S-version by the solution of ∇2(ψin,n, χin,n) ≡ (ζ, α) in D with a zero Dirichlet boundary value but is ill-defined in the C-version by the solution of ∇2(ψin,n, χin,n) ≡ (ζ, α) in D with a zero Neumann boundary value [see Eq. (2.7) of CK92b] without considering the two solvability conditions in Eqs. (5a,b). The ill-defined (ψin,n, χin,n) in the C-version has not caused significant inconsistency problems in previous applications to synoptic-scale flows in large-scale limited domains in which ∫D ζdx and ∫D αdx nearly vanish. Thus, the two solvability conditions in Eqs. (5a,b) can be approximately satisfied with a zero Neumann boundary value. Perhaps because of this, the inconsistency problem of the C-version was not noted or addressed when the S-version was refined (Boyd et al., 2013). This inconsistency problem can become serious for applications of the C-version method to complex flows with increased ∫D ζdx and/or ∫D αdx (as seen later from Table 2). In this case, the ill-defined (ψin,n, χin,n) in the C-version must be properly redefined, as shown below.
The inner part in the C-version can be properly redefined if it is further divided into two sub-parts, denoted by (ψin,s, χin,s) and (ψin,c, χin,c), respectively, while (ψin,s, χin,s) is well-defined by the solution of ∇2(ψin,s, χin,s) = (∫D ζdx, ∫D
αdx) in D with zero Dirichlet boundary value and (ψin,c, χin,c) is well defined by the solution of ∇2(ψin,c, χin,c) = (ζ – ∫D ζdx, α – ∫D αdx) in D with a zero Neumann boundary value. The spectral formulations derived for the S-versionꞌs (or C-versionꞌs) inner part in CK92a (or CK92b) can be then used to obtain (ψin,s, χin,s) [or (ψin,c, χin,c)] consistency with the solvability conditions in Eqs. (5a,b). The redefined inner part is thus given by (ψin,n, χin,n) = (ψin,s, χin,s) + (ψin,c, χin,c) in the first step. After this, the harmonic part, (ψh,a, χh,a), can be obtained in the second step using the spectral formulations derived in CK92b for C-version harmonic part. The C-version spectral method is modified in this way and implemented in this paper. According to CK92b, the C-version spectral method is more accurate and efficient than the S-version. The modified C-version spectral method, called the M-spectral method for short, is more accurate than the two original versions of the spectral method of CK92a, b especially for mesoscale domains, so the M-spectral method is used in this paper. Like the SOR-based method, the second part of the M-spectral method is also not very sensitive to increased variations of the original v across or along S. However, the first part of the M-spectral method is more sensitive than the SOR-based method to increased spatial variations of the original v inside D. The increased sensitivity is due mainly to the use of a zero Dirichlet (or Neumann) boundary value in defining and solving for (ψin,s, χin,s) [or (ψin,c, χin,c)] that can cause sharp near-boundary variations in the solution poorly resolved by the analysis grid, especially when (∫D ζdx, ∫D αdx) is not small or (ζ, α) has sharp near-boundary variations. This is because (ψin,s, χin,s) [or (ψin,c, χin,c)] has a zero Dirichlet (or Neumann) boundary value and thus must have sharp near-boundary variations in the boundary-normal (or boundary-parallel) direction so that ∇2(ψin,s, χin,s) [or ∇2(ψin,c, χin,c)] can have the same sharp near-boundary variations as (ζ, α) to satisfy the Poisson equation. Poorly resolved near-boundary variations in (ψin,n, χin,n) [= (ψin,s, χin,s) + (ψin,c, χin,c)] can generate additional errors in the boundary conditions for (ψh,a, χh,a) [see Eqs. (2.11)–(2.15) of CK92b], and these additional errors can propagate and spread throughout the domain when (ψh,a, χh,a) is solved iteratively. The above-explained additional errors and their propagation make the M-spectral method less accurate than the SOR-based method for the difficult cases considered in this paper (as shown later in Table 1). However, since the additional errors are generated from (ψin,n, χin,n) in the M-spectral method, the second part of the method for obtaining (ψh,a, χh,a) should still have an advantage if it is used not in the M-spectral method but in a hybrid method to obtain the externally induced (ψe, χe) after the internally induced (ψin,n, χin,n) is obtained by the integral method; this hybrid approach will also be considered in section 2.4.
-
Two hybrid methods are designed to reduce the increased loss of accuracy caused by increased complexities of the original v field, not only in D but also across or along S. The mathematical foundations for the hybrid methods are the same as the adopted three methods mentioned above. In these two hybrid methods, (ψ, χ) is divided into (ψin, χin) and (ψe, χe) for the internally and externally induced components, respectively, in the same way as in the integral method. To retain the strength and related advantages explained in section 2.1 for the integral method, the internally induced components, (ψin, χin), should be calculated in the first step in the same way as that in the integral method. The SOR-based method can be used adaptively for computing the externally induced ψe (with χe = 0) in the second step as described in section 2.2, while the M-spectral method for computing (ψe, χe) as described in section 2.3. This approach can improve the accuracy of (ψe, χe) computed in the second step because (as explained in sections 2.2 and 2.3), these two methods are less sensitive than the integral method to increased variations of original v across or along the domain boundary and their caused errors in solvability conditions and related boundary condition errors. The aforementioned hybrid methods are called the integral-SOR and integral-spectral methods, respectively. The detailed computational procedures of the two hybrid methods are shown in Fig. 1. Test results from these two hybrid methods in computing (ψ, χ) from original v fields with complex flow patterns in mesoscale domains will be shown in section 4 and compared with the results from three previous methods.
Domain | Method | SCC | RRD |
D1 | SOR-based | 0.9392 | 0.2820 |
M-spectral | 0.7836 | 0.5479 | |
Integral | 0.9870 | 0.1156 | |
Integral-spectral | 0.9867 | 0.1151 | |
Integral-SOR | 0.9958 | 0.0720 | |
D2 | SOR-based | 0.8644 | 0.2641 |
M-spectral | 0.8661 | 0.3680 | |
Integral | 0.9809 | 0.1115 | |
Integral-spectral | 0.9822 | 0.1032 | |
Integral-SOR | 0.9932 | 0.0648 | |
D3 | SOR-based | 0.8889 | 0.5024 |
M-spectral | 0.8753 | 0.5786 | |
Integral | 0.9740 | 0.2240 | |
Integral-spectral | 0.9633 | 0.2698 | |
Integral-SOR | 0.9909 | 0.1288 |