- Email: [email protected]
Boundary element methods for transient convective diffusion. Part II: 2D implementation
Comput. Methods Appl. Mech. Engrg. 192 (2003) 4299–4312 www.elsevier.com/locate/cma
Boundary element methods for transient convective diffusion. Part II: 2D implementation M.M. Grigoriev, G.F. Dargush
*
Department of Civil Engineering, State University of New York at Buffalo, 212 Ketter Hall, North Campus, Amherst, Buffalo, NY 14260, USA Received 22 November 2002; received in revised form 16 April 2003; accepted 21 May 2003
Abstract A boundary element method for transient convective diffusion phenomena presented in Part I of the paper is extended to two dimensional problems. We introduce a series representation for the transient convective kernel and perform a time integration for the double integrals to evaluate coefficients of the time-discrete boundary integral equation. The time-integrated kernels are evaluated for the linear, quadratic and quartic time interpolation functions utilized in the paper. Then, linear, quadratic and quartic boundary elements as well as bi-linear, bi-quadratic and biquartic volume cells are introduced to ensure proper resolution in space for the two-dimensional formulation. Due to the singular nature of the transient convective diffusion kernels, integration of the kernels over the boundary elements and volume cells requires a considerable effort to maintain a desired level of accuracy. We define influence domains due to time-integrated and time-delayed kernels arising for the surface and volume integrals, respectively. Note that the kernel influences are extremely localized due to the convective nature of the kernels, thus, the surface and volume integrations are performed only within these domains of influence. The localization of the kernels becomes more prominent as the Peclet number of the flow increases. Due to increasing sparsity of the global matrix, iterative solvers become the primary choice for the convective diffusion problems. 2003 Elsevier B.V. All rights reserved. Keywords: Boundary element methods; Unsteady convective diffusion; Higher-order methods
1. Introduction In Part I of this paper [1], we have presented a boundary element formulation for transient convective heat diffusion valid for problems with arbitrary dimensionality. We have introduced linear, quadratic, and quartic time interpolation functions and evaluated time-integrated convective kernels for the one-dimensional formulation. A full set of closed form time integrals have been presented for the first time. Then, we have defined the domains of kernel influences and obtained the boundaries of these domains. Since the
*
Corresponding author. Tel.: +1-716-645-2114; fax: +1-716-645-3711. E-mail address: [email protected]ffalo.edu (G.F. Dargush).
0045-7825/$ - see front matter 2003 Elsevier B.V. All rights reserved. doi:10.1016/S0045-7825(03)00389-X
4300
M.M. Grigoriev, G.F. Dargush / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4299–4312
contribution due to time-integrated and time-delayed kernels are localized within the influence domains, the segmentation of the volume cells and subsequent volume integration has been performed only locally. Part II extends the higher-order BEM to the two-dimensional transient convective diffusion problems. Since the closed form time integrals are hardly possible for two-dimensions, we introduce a series representation for the time-dependent kernels and provide time-integrated formulae for linear, quadratic and quartic time interpolation functions. Although the implicit time marching formulation utilized in this paper does not impose any restrictions on the time step size, the applicability of the numerical algorithm is restricted to the following time step sizes: Dt 6 Dtc ¼
4 ; V 2 Pe
ð1Þ
where V is the characteristic velocity of the flow. The use of time steps Dt less than the critical time step size Dtc ensures fast convergence of the series introduced for the time-dependent convective diffusion kernels in two dimensions. Although we demonstrate that larger time steps Dt P Dtc are possible, these computations require more terms in the series to attain the desired level of accuracy. Consequently, the boundary element method becomes less efficient and even fails to converge if the time steps are too large. Note that the critical time step size Dtc in (1) is independent of the boundary element mesh discretization. To ensure proper resolution in time and space, we consider linear, quadratic and quartic boundary elements to represent spatial variation of boundary temperatures and normal fluxes. Accordingly, bi-linear, bi-quadratic and bi-quartic volume cells are introduced to evaluate volume integrals arising in the timerecurring formulation. Due to the singular nature of the transient convective diffusion kernels, integration of the kernels over the boundary elements and volume cells requires a considerable effort to maintain a desired level of accuracy. We propose an adaptive approach for the subsegmentation of boundary elements and volume cells that provides us with very intelligent and numerically efficient algorithms for surface and volume integration. We should emphasize the importance of the adaptive approach for highly convective flows with the strong manifestation of the kernel singularity. The convective kernels are significant only at the vicinity of the pseudo-source point that is convected by the flow. The kernels decay rapidly away from the pseudo-source point and integration of the kernels is not needed beyond a certain range. Note that the area of the kernel influence at high Peclet numbers becomes confined to a tiny oval-shaped contour which is significantly less than the area of the smallest volume cells considered in the paper. Therefore, the subsegmentation and subsequent integration is performed only over boundary elements and volume cells influenced by the time-dependent kernels. Thus, the efficiency of the algorithm increases dramatically as the flow becomes more convective. In this paper, we extend the higher-order boundary element formulation presented in Part I to transient convective diffusion problems in two dimensions. The governing equation, integral formulation, and temporal discretization are presented in Part I of the paper. Here, Section 2 outlines the evaluation of the time-integrated kernels arising in the surface and volume integrals. Spatial discretization together with the boundary element and volume cell integrations is outlined in Section 3 and Appendix A. Section 4 discusses the assembly of the global matrix, while Section 5 presents some conclusions. Numerical results for both one- and two-dimensional problems are considered in Part III.
2. Time-integrated kernels 2.1. Introduction We should emphasize that the integral formulations and temporal discretizations presented in Sections 2 and 3 of Part I, respectively, are valid for one- and two-dimensional BEM formulations. As noted in Part I
M.M. Grigoriev, G.F. Dargush / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4299–4312
4301
of the paper, the difference between these formulations begins when evaluating time-integrated kernels. For the one-dimensional formulation, we have obtained the time integrals in the closed forms. Due to the time multiplier in the transient convective kernel, the closed form time integration is readily performed only for odd-dimensional problems, i.e., in one and three dimensions. Unfortunately, the closed form integration of the kernels over a finite time increment has not been possible for the two dimensional formulations. We should note that the analytical time integration is mandatory in order to develop efficient and accurate time-dependent BEM formulations. Thus, we introduce a series representation for the transient convective kernel gðy; tÞ and then perform the time integration over a finite time interval term by term. Although the series converge rapidly for relatively small time step sizes Dt 6 Dtc , where the critical time step size is given by (1), the radius of convergence is independent of the boundary element mesh discretization. 2.2. Evaluation of P(a)(y; t) and R(a)(y; t) We note that the evaluation of the time-integrated coefficients GðaÞ ðy; kDtÞ and F ðaÞ ðy; kDtÞ outlined in Appendix I of Part I is also valid for the two-dimensional formulation considered in this paper. However, evaluation of the following time integrals Z t ðaÞ P ðy; tÞ ¼ sa gðy; sÞ ds ð2Þ 0
and ðaÞ
R ðy; tÞ ¼
Z
t
sa fn ðy; sÞ ds
ð3Þ
0
for a ¼ 0; 1; 2; 3; 4 depends on the problem dimension. Since the straightforward closed form integration of (2) or (3) is not possible for any a ¼ 0; 1; 2; 3; 4, we substitute the kernel gðy; sÞ presented in Part I into integral (2): rV cos # Pe Z t Pe exp r2 Pe sV 2 Pe 2 P ðaÞ ðy; tÞ ¼ sa1 exp ds: ð4Þ 4p 4s 4 0 In (4), # is the angle between the radius vector yi and the velocity of the media Vi . Unfortunately, the closedform integration of (4) is possible only for some particular cases, e.g., t ! 1 (see [2,3] for details). For the general case, we introduce the following series representation [4] k V 2 Pe sk 1 X sV 2 Pe 4 : ð5Þ exp ¼ 4 k! k¼0 Note that the series (5) is valid for all
V 2 Pe 4
P 0, however, the convergence is fast only when
2
V Pe 6 1: 4 Substitution of (5) into (4) and its subsequent integration leads to the following rV cos # k 2 Pe X Pe exp 1 1 V 2 Pe r Pe 2 ðaÞ kþa P ðy; tÞ ¼ t Ekþaþ1 ; 4p k! 4 4t k¼0
ð6Þ
ð7Þ
where a ¼ 1, 0, 1, 2, 3, 4 and En ðÞ is the exponential integral [4]. The evaluation of the integral P ðaÞ ðy; tÞ requires evaluation of a series of the exponential integrals En ðzÞ for very large ranges of z P 0 and orders n > 0. The basic recurrence relation
4302
M.M. Grigoriev, G.F. Dargush / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4299–4312
1 Enþ1 ðzÞ ¼ ½expðzÞ zEn ðzÞ n
ð8Þ
for integer n > 0;
given in [4] is stable only for small parameters z < 1. In this paper, the exponential integrals En ðzÞ for n > 0 are evaluated using a very efficient and accurate algorithm provided by Gautschi [5]. Finally, the following equation RðaÞ ðy; tÞ ¼
yi ni ða1Þ Vi ni ðaÞ P P ðy; tÞ for a ¼ 0; 1; 2; 3; 4 ðy; tÞ 2 2
ð9Þ
can be used to evaluate integral (3).
3. Spatial discretization 3.1. Boundary elements and volume cells We divide the boundary into a number of boundary elements Cl , while the domain X is segmented into volume cells Xe . Linear, quadratic and quartic boundary elements are considered to represent a spatial variation of the temperatures T ðkÞ ðxÞ ¼ Na ðgÞ ðaÞ T ðkÞ
ð10Þ
and normal fluxes qðkÞ ðxÞ ¼ Na ðgÞ ðaÞ qðkÞ
ð11Þ
over boundary elements at any time level k. In (10) and (11), g ¼ gðxÞ represents an intrinsic coordinate, and a summation over the repeated index a is assumed. Linear, quadratic and quartic shape functions Na ðgÞ are defined as usual and may be found in the literature [6–9]. Accordingly, bi-linear, bi-quadratic and bi-quartic interpolation functions are utilized to represent a variation of the temperature over the volume cells Xe T ð0Þ ðxÞ ¼ Mb ðg1 ; g2 Þ ðbÞ T ð0Þ ;
ð12Þ
where shape functions Mb ðg1 ; g2 Þ are also defined as usual (see [10] for details). Similarly, we imply a summation over the index b. The typical boundary element mesh utilized for the two-dimensional problems considered in this paper is shown in Fig. 1. Once the boundary element (10), (11) and volume cell (12) approximations are introduced, we discretize the boundary integral equation (14) in Part I to the following form: Z Z X ðbÞ ðaÞ cðnÞT ðkÞ ðnÞ þ qn Nb ðgÞGðaÞ ðy; kDtÞ dCðxÞ þ ðbÞ T ðaÞ Nb ðgÞF a ðy; kDtÞ dCðxÞ þ ðbÞ T ðaÞ
Z
l¼1
Cl
Cl
Cl
X Z ðbÞ ð0Þ a Nb ðgÞun ðxÞG ðy; kDtÞ dCðxÞ ¼ T Mb ðg1 ; g2 Þgðy; kDtÞ dXðxÞ: e¼1
ð13Þ
Xe
Here, the summation over repeated indices a ¼ 0; 14 ; 12 ; 34 ; 1 is implied. Note that the discrete boundary integral equation (13) requires evaluation of the integrals over the boundary elements Cl and volume cells Xe . It is well-known that the accuracy of the boundary element method greatly depends on the integration accuracy over the boundary elements and volume cells. The issue of the accurate integration becomes more acute at higher Peclet number when the kernels are extremely localized with very large gradients.
M.M. Grigoriev, G.F. Dargush / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4299–4312
4303
Fig. 1. Mesh of 4 · 4 regular rectangles; quartic (five-noded) boundary elements (open circles) and bi-quartic (twenty-five noded) volume cells (pluses).
3.2. Evaluation of surface integrals As follows from Appendix I in Part I, evaluation of the surface integrals involved in the discrete boundary integral equation (13) leads to the following integrals over the boundary element Cl : Z ðaÞ Ub ðn; tÞ ¼ P ðaÞ ðy; tÞNb ðgÞ dCðxÞ ð14Þ Cl
and ðaÞ Wb ðn; tÞ
¼
Z
RðaÞ ðy; tÞNb ðgÞ dCðxÞ
ð15Þ
Cl
for a ¼ 0; 1; 2; 3; 4 and b ¼ 1; 2; 3; 4; 5. We notice that the time-integrated terms P ðaÞ ðy; tÞ and RðaÞ ðy; tÞ are significant only at the proximity of the source point n. Similar to the one-dimensional formulation, we need to establish the domain of influence for P ðaÞ ðy; tÞ, RðaÞ ðy; tÞ and perform integration only within that domain. In Section 2, we have obtained the following series representation for P ðaÞ ðy; tÞ: P ðaÞ ðy; tÞ ¼
e
rV cos # Pe 2
4p
1 Pe X ðaÞ Sk ;
for a ¼ 1; 0; 1; 2; 3; 4;
ð16Þ
k¼0
where k 2 1 V 2 Pe r Pe kþa ¼ t Ekþaþ1 : k! 4 4t 2 In (17), Ekþaþ1 r 4tPe is the exponential integral [4]. Since exponential integral [10] ðaÞ Sk
ð17Þ tV 2 Pe 4
6 1, and using the following property of the
4304
M.M. Grigoriev, G.F. Dargush / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4299–4312
Ekþaþ2
r2 Pe 4t
we get S ðaÞ kþ1 ðaÞ 6 1 S
< Ekþaþ1
r2 Pe ; 4t
for all k P 0 and
for k P 0
and
a ¼ 1; 0; 1; 2; 3; 4
ð18Þ
a ¼ 1; 0; 1; 2; 3; 4:
ð19Þ
k
Moreover, we can show that for all a ¼ 1; 0; 1; 2; 3; 4 S ðaÞ kþ1 lim ¼ 0: k!1 S ðaÞ k
ð20Þ
Now that we showed the convergence of the series (16), we need an appropriate estimator for the domain of the influence to perform the surface integration. Obviously, a good candidate for estimating the P ðaÞ ðy; tÞ term would be the first term in the series (16) for a ¼ 0 since it follows from (20) that ð0Þ
ðaÞ
S 0 P S0
for a ¼ 1; 2; 3; 4:
ð21Þ
ðaÞ
ðaÞ
ðaÞ
Using (9) to express R ðy; tÞ in terms of P ðy; tÞ, we notice that the domain of influence for R ðy; tÞ can be ð1Þ estimated using the first term of the series representation for P ð1Þ ðy; tÞ, namely S0 . Taking into account ð1Þ the series representation (16), we can show that the estimator for the P ðy; tÞ term is stronger than that for P ð0Þ ðy; tÞ. Finally, the domain of influence Rs ðyÞ for the time-integrated terms P ðaÞ ðy; tÞ and RðaÞ ðy; tÞ can be estimated by the first term in the series for P ð1Þ ðy; tÞ as follows: Pe rV cos # r2 Pe Pe rs ðyÞ ¼ 2 exp : ð22Þ pr 2 4t The domain of influence is then bounded by the curve rs ðyÞ ¼ es ;
ð23Þ
where es represents a small tolerance for surface integrals, e.g., es ¼ 1030 , to ensure the contribution of the time-integrated kernels P a ðy; tÞ and RðaÞ ðy; tÞ from outside the domain Rs ðyÞ is negligible. Thus, the integration over the boundary elements (14) and (15) is performed only within the domain Rs ðyÞ associated with the source point n. No integration is needed if boundary element Cl lies outside the domain Rs ðyÞ, i.e., Cl ðxÞ \ Rs ðx nÞ ¼ 0. Notice that when the collocation point n does not lie on the boundary element Cl , the time-integrated kernels P ðaÞ ðy; tÞ and RðaÞ ðy; tÞ involve no singularity and hence can be evaluated using a standard Gaussian quadrature. However, to perform the integration in a more efficient way, we could either introduce Ns contours of pre-set levels ðiÞ rðiÞ s P es ;
i ¼ 1; 2; . . . ; Ns
and
rsðNs Þ ¼ es
ð24Þ
and segment the boundary element Cl via these contours. Or, alternatively, the contours can be evaluated directly from the outer-level contour via the following non-uniform scaling rðiÞ ð#Þ ¼ Hi rðNs Þ ð#Þ;
i ¼ 1; 2; . . . ; Ns1 ;
ð25Þ
sÞ . In the where 0 < Hi < 1, and rðNs Þ ð#Þ is the distance from the source point n to the outer-level contour rðN s former approach, difficulties arise in specifying the preset level contours to attain a desired clustering of the segments for different combinations of t and Pe as limy!0 rs ðyÞ ! 1. Thus, we utilize the latter approach throughout the study. Note that the parameter Hi is much less sensitive to the Peclet number than the preset levels in the former approach.
M.M. Grigoriev, G.F. Dargush / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4299–4312
4305
For both approaches considered above, the non-uniform clustering of the segments allows very efficient and accurate integration over the non-singular boundary element Cl influenced by the transient kernel (Fig. 2). Variable order Gaussian quadrature is utilized for every segment of interest to ensure a high level of accuracy. However, when the collocation point n lies on the straight boundary element Cl , the evaluation of the boundary integrals (14) and (15) involves a singularity as y ! 0. Analysis of terms in (9), (16), and (17) shows that the contribution due to kernel P ð1Þ ðy; tÞ vanishes in (9) since yi ni ¼ 0 for the singular boundary element, while the time-integrated kernels P ð0Þ ðy; tÞ lead to log-type singularities. All other terms result in non-singular integrals. Moreover, a boundary element integral involving P ð0Þ ðy; tÞ is singular only when the shape function Nb ðgÞ coincides with the local number of the collocation node c. Otherwise, the integral is non-singular due to the following property of the shape function 1 if c ¼ b; Nb ðgðy ðcÞ ÞÞ ¼ ð26Þ 0 if c 6¼ b: The evaluation of the non-singular integrals is similar to the algorithm described above. ð0Þ In order to remove the singularity from the integral Ub ðn; tÞ, we represent it as follows: Z Z Z P ð0Þ ðy ðbÞ ; tÞNb ðgÞ dCðxÞ ¼ P ð0Þ ðy ðbÞ ; tÞðNb ðgÞ 1Þ dCðxÞ þ P ð0Þ ðy ðbÞ ; tÞ dCðxÞ: Cl
Cl
ð27Þ
Cl
Notice that the first integral on the right-hand side of (27) is no longer singular and can be integrated using the approach outlined above. The second integral in (27) can be split into singular and non-singular parts: "Z 2 Z yV Pe r Pe ð0Þ ðbÞ i2 i Pe P ðy ; tÞ dCðxÞ ¼ e E1 dCðxÞ 4p Cl 4t Cl # k Z 2 1 X yV tk V 2 Pe r Pe i2 i Pe þ e Ekþ1 dCðxÞ : ð28Þ 4 4t k! Cl k¼1
Fig. 2. Schematics of the boundary element segmentation using scaled contours of the time-integrated kernels.
4306
M.M. Grigoriev, G.F. Dargush / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4299–4312
Fig. 3. Scheme for log-singular integration over the boundary element.
Now the log-singular integral on the right-hand-side of (28) can be represented in the following form using the notations introduced in Fig. 3: pffiffiffiffi
2 Z yV r Pe ez ð1 eÞz 4t i2 i Pe ; qk þ J ; qk : ð29Þ J e E1 dCðxÞ ¼ 4t k k Pe Cl pffiffiffiffiffiffi 2 V1 In (29), k ¼ tPe, q ¼ n1 V2 n , z ¼ LV2 Pe, and V Z a J ða; bÞ ¼ ebx E1 ðx2 Þ dx: ð30Þ 0
Unfortunately, no closed-form integral for (30) is available in the published literature except for the particular case of a ¼ b. In Appendix A we present an algorithm to numerically evaluate integral (30) for all possible cases of a > 0 and 1 6 b 6 1. Meanwhile, the non-singular integral in (28) is evaluated using the integration algorithm described earlier in this section. 3.3. Evaluation of volume integrals Although integration over the volume cells involves no singularity, it requires careful attention at small time steps for moderate and high Peclet numbers due to the nature of the kernel functions. Fig. 4 shows an upstream propagation of the convective kernel gðy; tÞ for the source located at the point of origin. Notice that the bell-shaped response is symmetric with respect to the axis of symmetry translating upstream with the speed V . When the Peclet number of the flow is relatively small, the kernel decays quickly due to large diffusion. For higher Peclet numbers, the response becomes more localized around the axis of symmetry since the convection dominates over the diffusion. Thus, it is mandatory to develop an integration algorithm that tracks the propagation of the transient convective kernel. Obviously, a straightforward integration of the kernel over volume cells would be extremely inefficient, especially for high Peclet numbers. Similar to the one-dimensional formulation, we introduce the domain of influence Rv ðyÞ for the transient convective kernel at the current time level. The contribution of the kernel gðy; tÞ is significant only within the domain of influence, whereas it is negligible outside Rv ðyÞ. The domain of influence is bounded by the contour satisfying the following gðy; tÞ ¼ ev :
ð31Þ
M.M. Grigoriev, G.F. Dargush / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4299–4312
4307
Fig. 4. Surface plots of the transient convective g-kernel for Pe ¼ 1000, V1 ¼ 1, V2 ¼ 0.
Here, ev ¼ 1030 is a small parameter to ensure the kernel is significant only within the domain of influence, and the transient convective kernel is scaled to provide 0 < gðy; tÞ 6 1 " # gðy; tÞ ðx1 n1 þ V1 tÞ2 þ ðx2 n2 þ V2 tÞ2 ¼ exp Pe : ð32Þ gðy; tÞ ¼ gðVt; tÞ 4t Introducing a frame of reference ni ¼ ni Vi t
ð33Þ
translating with the kernel, we show that contours of constant values of gðy; tÞ are concentric circles of radius rv ! 4t 1 2 rv ¼ ln : ð34Þ Pe gðy; tÞ Obviously, the center of circles coincides with the axis of symmetry. Fig. 5 shows the propagation of the circular domains of the kernel influence upstream from the source point ni . We should emphasize that the integration over a volume cell Xe is performed if any part of the volume cell is inside the domain of influence Rv ðyÞ: Xe \ Rv ðyÞ 6¼ 0: Otherwise, no integration over Xe is needed. If the volume cell is influenced by the transient convective kernel, we utilize a segmentation template shown in Fig. 6. We introduce a grid in polar coordinates ðr; hÞ so that the point of origin coincides with the axis of symmetry. Angular segmentation is uniform with the increment Dh. For the radial segmentation, we introduce Nv pre-set levels for the scaled kernel gi ðy; tÞ i ðy; tÞ < 1; eðiÞ v 6g
i ¼ 1; 2; . . . ; Nv :
ð35Þ
Here, the radius of each circular level gi ðy; tÞ is evaluated directly from (35). Note that the outer-level circle gNs ðy; tÞ is the boundary of the kernel influence. Generally, the volume cell Xe is not aligned with the segmentation grid and, thus, some of the grid segments generated in the polar coordinate system ðr; hÞ require a volume integration only over the parts
4308
M.M. Grigoriev, G.F. Dargush / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4299–4312
Fig. 5. Circular domains of kernel influence propagate upstream from the source point.
Fig. 6. Schematics of the volume cell segmentation influenced by the kernel. Circular domain of influence is represented by polygons.
that belong to the volume cell. For the sake of simplicity in evaluating the intersections between influence domain and volume cells, the circular domains are replaced with polygons as shown in Fig. 6. For relatively small increments Dh, the accuracy of the integration is not impacted by this approximation. Note that the volume integrals are evaluated using a standard Gaussian quadrature. Similar to the boundary element integration, we utilize a variable order of the Gaussian quadrature for every volume segment of interest to ensure a high level of accuracy.
M.M. Grigoriev, G.F. Dargush / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4299–4312
4309
4. Assembly The assembly procedure for the BEM approach is standard with no additional complications arising from the use of the higher-order time functions. In contrast to assembly procedure for time-dependent heat diffusion presented by Grigoriev and Dargush [11], the transient convective diffusion leads to a sparse global matrix for both one- and two-dimensional formulations. Since the global matrix for one-dimensional problems is very compact, the choice of the solution algorithm has not been a critical issue there. However, the size of the global matrix for the two dimensions can be as large as 100,000 for the problems considered in this study. For the problems of practical importance, the number of degrees of freedom could easily go beyond a million. We should emphasize that the sparsity of the global matrix increases dramatically with an increase of the Peclet number. The non-zero elements of the matrix are clustered around the main diagonal. Off-diagonal elements decay rapidly or are virtually zero for the entire matrix. Therefore, iterative algorithms would be of the primary choice to solve the global set of equations. In this paper, we utilize the preconditioned biconjugate gradient method for sparse linear systems presented by Press et al. [12]. We should note that the use of the main diagonal as a preconditioner results in a stable, accurate and efficient solution algorithm even for very high Peclet numbers.
5. Conclusions The higher-order boundary element methods for time-dependent convective diffusion phenomena presented in Part I of the paper are extended to the two-dimensional problems. First, the transient convective kernel gðy; tÞ is represented as an infinite series. Although the series converges for all time step sizes, the applicability of the numerical algorithm is constrained by inequality (1). We emphasize that the critical time step size Dtc following from (1) is not subject to the mesh discretization contrary to the well-known constraint on local mesh-based Peclet number Pec 6 2 when no upwinding is utilized. Note that the convective kernels used in this paper provide an analytical upwinding for any Peclet number. To relieve the restriction or, rather, recommendation on the time step size of the algorithm, another form of the series representation for the kernels is needed. However, the implementation of a suitable alternative goes beyond the scope of this paper. Here the series of the convective kernel is integrated term-by-term over finite time intervals for all time interpolation functions. To facilitate accurate evaluation of the discrete integral equation coefficients over both singular and non-singular boundary elements, we have developed and implemented an adaptive algorithm that accounts for the complex nature of the time-dependent convective kernels. The algorithm tracks the domains influenced by the kernels entering the surface and volume integrals. Accordingly, it provides an automatic segmentation of the influence domains to allow efficient and accurate numerical integrations of the kernels over higher-order boundary elements and volume cells. The numerical results for the transient convective heat diffusion are presented and discussed in Part III.
Acknowledgements The work described in this paper was partially supported by the National Science Foundation under Grant CMS 9700387. The authors gratefully acknowledge this support. Furthermore, the authors are grateful to the reviewer for many valuable comments that helped improve the paper.
4310
M.M. Grigoriev, G.F. Dargush / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4299–4312
Appendix A. Evaluation of integral J(a; b) A.1. Small parameter a 6 1 Here we present an algorithm to numerically evaluate integral (30) given by Z a J ða; bÞ ¼ ebx E1 ðx2 Þ dx:
ðA:1Þ
0
Note that an accurate evaluation of integral (A.1) is required for any combinations of a > 0 and 1 6 b 6 1. Also, notice that the exponential integral E1 ðx2 Þ leads to a log-type singularity as x ! 0. We should emphasize that the straightforward substitution of the following series representation for the exponential integral [3] E1 ðx2 Þ ¼ c 2 lnðxÞ
1 X ð1Þk x2k kk! k¼1
ðA:2Þ
into (A.1) and a subsequent integration allows us to obtain 1 X c ð1 expðabÞÞ 2 ð1Þ B2k ða; bÞ þ ½c þ ln jbj þ expðabÞ lnðaÞ þ E1 ðabÞ : b b kk! k¼1 k
J ða; bÞ ¼
ðA:3Þ In (A.2) and (A.3), c ¼ 0:57721566. . . is the EulerÕs constant [4], and the integral Z a Bk ða; bÞ ¼ ebx xk dx
ðA:4Þ
0
can be expressed by the following recurrence formula Bk ða; bÞ ¼
1 ½kBk1 ða; bÞ ak expðabÞ b
B2 ða; bÞ ¼
2 expðabÞð2 þ 2ab þ a2 b2 Þ b3
for k > 2:
ðA:5Þ
Here, for k ¼ 2:
Since Bk ða; bÞ are needed only for even k, we modify (A.5) to this form
1 kðk 1Þ k Bk2 ða; bÞ ak1 a þ Bk ða; bÞ ¼ expð abÞ for k P 4: b b b
ðA:6Þ
ðA:7Þ
Unfortunately, the recurrence evaluation of (A.7) is stable only for small parameters a 6 1. For higher a, the use of (A.7) to evaluate the integral J ða; bÞ results in large errors. Thus, this algorithm is limited to a 6 1 only.
A.2. Large parameter 1 < a 6 10 Let us introduce parameter k and recast integral (A.1) to the following form Z a J ða; b; kÞ ¼ expðbxÞE1 ðk2 x2 Þ dx: 0
ðA:8Þ
M.M. Grigoriev, G.F. Dargush / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4299–4312
4311
Differentiation of (A.8) with respect to k gives Z a o 2 J ða; b; kÞ ¼ expðbx k2 x2 Þ dx: ok k 0 Integration of (A.9) is straightforward and leads to: 2 pffiffiffi b p exp o b b 4k2 J ða; b; kÞ ¼ erf erf þ ak : ok 2k 2k k2 We notice that Z 1 o J ða; b; kÞ dk ¼ J ða; bÞ ok 1
ðA:9Þ
ðA:10Þ
ðA:11Þ
since lim J ða; b; kÞ ¼ 0
ðA:12Þ
J ða; b; 1Þ ¼ J ða; bÞ:
ðA:13Þ
k!1
and Integration of (A.10) over the parameter k and the use of (A.11) allow us to recast the desired integral in the following form: # pffiffiffi " Z b 2 2 p ab 2 J ða; bÞ ¼ expðy Þ erf þ y erfðyÞ dy : ðA:14Þ b 2y 0 We should note that the behavior of the integrand in (A.14) is no longer singular. Thus, the Gaussian quadrature can be applied to numerically evaluate the integral in (A.14). Notice that the integral (A.14) can also be used for small 0 < a < 1. For the limiting case of b ! 0, the closed form integration of (A.1) gives: pffiffiffi J ða; 0Þ ¼ aE1 ða2 Þ þ p erfðaÞ: ðA:15Þ A.3. Very large parameter a > 10 For large parameters a > 10, integral (A.1) can be represented in the following form: Z a J ða; bÞ ¼ J ð10; bÞ þ ebx E1 ðx2 Þ dx:
ðA:16Þ
10
It is possible to show that the value of the second integral in (A.16) is negligible for all 1 6 b 6 1. Thus, only the contribution of the first integral J ð10; bÞ is considered for large a.
References [1] M.M. Grigoriev, G.F. Dargush, Boundary element methods for transient convective diffusion. Part I: General formulation and 1D implementation, Comput. Methods Appl. Mech. Engrg. 192 (2003) 4281–4298. [2] A.P. Prudnikov, Yu.A. Brychkov, O.I. Marichev, Integrals and Series, vol. 1: Elementary Functions, Gordon and Breach Science Publishers, New York, 1986.
4312
M.M. Grigoriev, G.F. Dargush / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4299–4312
[3] I.S. Gradshtein, I.M. Ryzhik, Tables of Integrals, Series and Products, Academic Press, 1980. [4] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, Dover, New York, 1964. [5] W. Gautschi, Algorithm 471: Exponential integrals, Commun. ACM 16 (1973) 761–763. [6] C.A. Brebbia, S. Walker, Boundary Element Techniques in Engineering, Newness-Butterworths, London, 1980. [7] P.K. Banerjee, R. Butterfield, Boundary Element Methods in Engineering Science, McGraw-Hill, London, 1981. [8] C.A. Brebbia, J.C.F. Telles, L.C. Wrobel, Boundary Element Techniques, Springer, Berlin, 1984. [9] P.K. Banerjee, The Boundary Element Methods in Engineering, McGraw-Hill, London, New York, 1994. [10] K.J. Bathe, Finite Element Procedures, Prentice Halls, Englewood Cliffs, 1996. [11] M.M. Grigoriev, G.F. Dargush, Higher-order boundary element methods for transient diffusion problems. Part I: Bounded flux formulation, Int. J. Num. Methods Engrg. 55 (2002) 1–40. [12] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical recipes in C. The Art of scientific computing, Cambridge University Press, Cambridge, 1992.
Boundary element methods for transient convective diffusion. Part II: 2D implementation M.M. Grigoriev, G.F. Dargush
*
Department of Civil Engineering, State University of New York at Buffalo, 212 Ketter Hall, North Campus, Amherst, Buffalo, NY 14260, USA Received 22 November 2002; received in revised form 16 April 2003; accepted 21 May 2003
Abstract A boundary element method for transient convective diffusion phenomena presented in Part I of the paper is extended to two dimensional problems. We introduce a series representation for the transient convective kernel and perform a time integration for the double integrals to evaluate coefficients of the time-discrete boundary integral equation. The time-integrated kernels are evaluated for the linear, quadratic and quartic time interpolation functions utilized in the paper. Then, linear, quadratic and quartic boundary elements as well as bi-linear, bi-quadratic and biquartic volume cells are introduced to ensure proper resolution in space for the two-dimensional formulation. Due to the singular nature of the transient convective diffusion kernels, integration of the kernels over the boundary elements and volume cells requires a considerable effort to maintain a desired level of accuracy. We define influence domains due to time-integrated and time-delayed kernels arising for the surface and volume integrals, respectively. Note that the kernel influences are extremely localized due to the convective nature of the kernels, thus, the surface and volume integrations are performed only within these domains of influence. The localization of the kernels becomes more prominent as the Peclet number of the flow increases. Due to increasing sparsity of the global matrix, iterative solvers become the primary choice for the convective diffusion problems. 2003 Elsevier B.V. All rights reserved. Keywords: Boundary element methods; Unsteady convective diffusion; Higher-order methods
1. Introduction In Part I of this paper [1], we have presented a boundary element formulation for transient convective heat diffusion valid for problems with arbitrary dimensionality. We have introduced linear, quadratic, and quartic time interpolation functions and evaluated time-integrated convective kernels for the one-dimensional formulation. A full set of closed form time integrals have been presented for the first time. Then, we have defined the domains of kernel influences and obtained the boundaries of these domains. Since the
*
Corresponding author. Tel.: +1-716-645-2114; fax: +1-716-645-3711. E-mail address: [email protected]ffalo.edu (G.F. Dargush).
0045-7825/$ - see front matter 2003 Elsevier B.V. All rights reserved. doi:10.1016/S0045-7825(03)00389-X
4300
M.M. Grigoriev, G.F. Dargush / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4299–4312
contribution due to time-integrated and time-delayed kernels are localized within the influence domains, the segmentation of the volume cells and subsequent volume integration has been performed only locally. Part II extends the higher-order BEM to the two-dimensional transient convective diffusion problems. Since the closed form time integrals are hardly possible for two-dimensions, we introduce a series representation for the time-dependent kernels and provide time-integrated formulae for linear, quadratic and quartic time interpolation functions. Although the implicit time marching formulation utilized in this paper does not impose any restrictions on the time step size, the applicability of the numerical algorithm is restricted to the following time step sizes: Dt 6 Dtc ¼
4 ; V 2 Pe
ð1Þ
where V is the characteristic velocity of the flow. The use of time steps Dt less than the critical time step size Dtc ensures fast convergence of the series introduced for the time-dependent convective diffusion kernels in two dimensions. Although we demonstrate that larger time steps Dt P Dtc are possible, these computations require more terms in the series to attain the desired level of accuracy. Consequently, the boundary element method becomes less efficient and even fails to converge if the time steps are too large. Note that the critical time step size Dtc in (1) is independent of the boundary element mesh discretization. To ensure proper resolution in time and space, we consider linear, quadratic and quartic boundary elements to represent spatial variation of boundary temperatures and normal fluxes. Accordingly, bi-linear, bi-quadratic and bi-quartic volume cells are introduced to evaluate volume integrals arising in the timerecurring formulation. Due to the singular nature of the transient convective diffusion kernels, integration of the kernels over the boundary elements and volume cells requires a considerable effort to maintain a desired level of accuracy. We propose an adaptive approach for the subsegmentation of boundary elements and volume cells that provides us with very intelligent and numerically efficient algorithms for surface and volume integration. We should emphasize the importance of the adaptive approach for highly convective flows with the strong manifestation of the kernel singularity. The convective kernels are significant only at the vicinity of the pseudo-source point that is convected by the flow. The kernels decay rapidly away from the pseudo-source point and integration of the kernels is not needed beyond a certain range. Note that the area of the kernel influence at high Peclet numbers becomes confined to a tiny oval-shaped contour which is significantly less than the area of the smallest volume cells considered in the paper. Therefore, the subsegmentation and subsequent integration is performed only over boundary elements and volume cells influenced by the time-dependent kernels. Thus, the efficiency of the algorithm increases dramatically as the flow becomes more convective. In this paper, we extend the higher-order boundary element formulation presented in Part I to transient convective diffusion problems in two dimensions. The governing equation, integral formulation, and temporal discretization are presented in Part I of the paper. Here, Section 2 outlines the evaluation of the time-integrated kernels arising in the surface and volume integrals. Spatial discretization together with the boundary element and volume cell integrations is outlined in Section 3 and Appendix A. Section 4 discusses the assembly of the global matrix, while Section 5 presents some conclusions. Numerical results for both one- and two-dimensional problems are considered in Part III.
2. Time-integrated kernels 2.1. Introduction We should emphasize that the integral formulations and temporal discretizations presented in Sections 2 and 3 of Part I, respectively, are valid for one- and two-dimensional BEM formulations. As noted in Part I
M.M. Grigoriev, G.F. Dargush / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4299–4312
4301
of the paper, the difference between these formulations begins when evaluating time-integrated kernels. For the one-dimensional formulation, we have obtained the time integrals in the closed forms. Due to the time multiplier in the transient convective kernel, the closed form time integration is readily performed only for odd-dimensional problems, i.e., in one and three dimensions. Unfortunately, the closed form integration of the kernels over a finite time increment has not been possible for the two dimensional formulations. We should note that the analytical time integration is mandatory in order to develop efficient and accurate time-dependent BEM formulations. Thus, we introduce a series representation for the transient convective kernel gðy; tÞ and then perform the time integration over a finite time interval term by term. Although the series converge rapidly for relatively small time step sizes Dt 6 Dtc , where the critical time step size is given by (1), the radius of convergence is independent of the boundary element mesh discretization. 2.2. Evaluation of P(a)(y; t) and R(a)(y; t) We note that the evaluation of the time-integrated coefficients GðaÞ ðy; kDtÞ and F ðaÞ ðy; kDtÞ outlined in Appendix I of Part I is also valid for the two-dimensional formulation considered in this paper. However, evaluation of the following time integrals Z t ðaÞ P ðy; tÞ ¼ sa gðy; sÞ ds ð2Þ 0
and ðaÞ
R ðy; tÞ ¼
Z
t
sa fn ðy; sÞ ds
ð3Þ
0
for a ¼ 0; 1; 2; 3; 4 depends on the problem dimension. Since the straightforward closed form integration of (2) or (3) is not possible for any a ¼ 0; 1; 2; 3; 4, we substitute the kernel gðy; sÞ presented in Part I into integral (2): rV cos # Pe Z t Pe exp r2 Pe sV 2 Pe 2 P ðaÞ ðy; tÞ ¼ sa1 exp ds: ð4Þ 4p 4s 4 0 In (4), # is the angle between the radius vector yi and the velocity of the media Vi . Unfortunately, the closedform integration of (4) is possible only for some particular cases, e.g., t ! 1 (see [2,3] for details). For the general case, we introduce the following series representation [4] k V 2 Pe sk 1 X sV 2 Pe 4 : ð5Þ exp ¼ 4 k! k¼0 Note that the series (5) is valid for all
V 2 Pe 4
P 0, however, the convergence is fast only when
2
V Pe 6 1: 4 Substitution of (5) into (4) and its subsequent integration leads to the following rV cos # k 2 Pe X Pe exp 1 1 V 2 Pe r Pe 2 ðaÞ kþa P ðy; tÞ ¼ t Ekþaþ1 ; 4p k! 4 4t k¼0
ð6Þ
ð7Þ
where a ¼ 1, 0, 1, 2, 3, 4 and En ðÞ is the exponential integral [4]. The evaluation of the integral P ðaÞ ðy; tÞ requires evaluation of a series of the exponential integrals En ðzÞ for very large ranges of z P 0 and orders n > 0. The basic recurrence relation
4302
M.M. Grigoriev, G.F. Dargush / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4299–4312
1 Enþ1 ðzÞ ¼ ½expðzÞ zEn ðzÞ n
ð8Þ
for integer n > 0;
given in [4] is stable only for small parameters z < 1. In this paper, the exponential integrals En ðzÞ for n > 0 are evaluated using a very efficient and accurate algorithm provided by Gautschi [5]. Finally, the following equation RðaÞ ðy; tÞ ¼
yi ni ða1Þ Vi ni ðaÞ P P ðy; tÞ for a ¼ 0; 1; 2; 3; 4 ðy; tÞ 2 2
ð9Þ
can be used to evaluate integral (3).
3. Spatial discretization 3.1. Boundary elements and volume cells We divide the boundary into a number of boundary elements Cl , while the domain X is segmented into volume cells Xe . Linear, quadratic and quartic boundary elements are considered to represent a spatial variation of the temperatures T ðkÞ ðxÞ ¼ Na ðgÞ ðaÞ T ðkÞ
ð10Þ
and normal fluxes qðkÞ ðxÞ ¼ Na ðgÞ ðaÞ qðkÞ
ð11Þ
over boundary elements at any time level k. In (10) and (11), g ¼ gðxÞ represents an intrinsic coordinate, and a summation over the repeated index a is assumed. Linear, quadratic and quartic shape functions Na ðgÞ are defined as usual and may be found in the literature [6–9]. Accordingly, bi-linear, bi-quadratic and bi-quartic interpolation functions are utilized to represent a variation of the temperature over the volume cells Xe T ð0Þ ðxÞ ¼ Mb ðg1 ; g2 Þ ðbÞ T ð0Þ ;
ð12Þ
where shape functions Mb ðg1 ; g2 Þ are also defined as usual (see [10] for details). Similarly, we imply a summation over the index b. The typical boundary element mesh utilized for the two-dimensional problems considered in this paper is shown in Fig. 1. Once the boundary element (10), (11) and volume cell (12) approximations are introduced, we discretize the boundary integral equation (14) in Part I to the following form: Z Z X ðbÞ ðaÞ cðnÞT ðkÞ ðnÞ þ qn Nb ðgÞGðaÞ ðy; kDtÞ dCðxÞ þ ðbÞ T ðaÞ Nb ðgÞF a ðy; kDtÞ dCðxÞ þ ðbÞ T ðaÞ
Z
l¼1
Cl
Cl
Cl
X Z ðbÞ ð0Þ a Nb ðgÞun ðxÞG ðy; kDtÞ dCðxÞ ¼ T Mb ðg1 ; g2 Þgðy; kDtÞ dXðxÞ: e¼1
ð13Þ
Xe
Here, the summation over repeated indices a ¼ 0; 14 ; 12 ; 34 ; 1 is implied. Note that the discrete boundary integral equation (13) requires evaluation of the integrals over the boundary elements Cl and volume cells Xe . It is well-known that the accuracy of the boundary element method greatly depends on the integration accuracy over the boundary elements and volume cells. The issue of the accurate integration becomes more acute at higher Peclet number when the kernels are extremely localized with very large gradients.
M.M. Grigoriev, G.F. Dargush / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4299–4312
4303
Fig. 1. Mesh of 4 · 4 regular rectangles; quartic (five-noded) boundary elements (open circles) and bi-quartic (twenty-five noded) volume cells (pluses).
3.2. Evaluation of surface integrals As follows from Appendix I in Part I, evaluation of the surface integrals involved in the discrete boundary integral equation (13) leads to the following integrals over the boundary element Cl : Z ðaÞ Ub ðn; tÞ ¼ P ðaÞ ðy; tÞNb ðgÞ dCðxÞ ð14Þ Cl
and ðaÞ Wb ðn; tÞ
¼
Z
RðaÞ ðy; tÞNb ðgÞ dCðxÞ
ð15Þ
Cl
for a ¼ 0; 1; 2; 3; 4 and b ¼ 1; 2; 3; 4; 5. We notice that the time-integrated terms P ðaÞ ðy; tÞ and RðaÞ ðy; tÞ are significant only at the proximity of the source point n. Similar to the one-dimensional formulation, we need to establish the domain of influence for P ðaÞ ðy; tÞ, RðaÞ ðy; tÞ and perform integration only within that domain. In Section 2, we have obtained the following series representation for P ðaÞ ðy; tÞ: P ðaÞ ðy; tÞ ¼
e
rV cos # Pe 2
4p
1 Pe X ðaÞ Sk ;
for a ¼ 1; 0; 1; 2; 3; 4;
ð16Þ
k¼0
where k 2 1 V 2 Pe r Pe kþa ¼ t Ekþaþ1 : k! 4 4t 2 In (17), Ekþaþ1 r 4tPe is the exponential integral [4]. Since exponential integral [10] ðaÞ Sk
ð17Þ tV 2 Pe 4
6 1, and using the following property of the
4304
M.M. Grigoriev, G.F. Dargush / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4299–4312
Ekþaþ2
r2 Pe 4t
we get S ðaÞ kþ1 ðaÞ 6 1 S
< Ekþaþ1
r2 Pe ; 4t
for all k P 0 and
for k P 0
and
a ¼ 1; 0; 1; 2; 3; 4
ð18Þ
a ¼ 1; 0; 1; 2; 3; 4:
ð19Þ
k
Moreover, we can show that for all a ¼ 1; 0; 1; 2; 3; 4 S ðaÞ kþ1 lim ¼ 0: k!1 S ðaÞ k
ð20Þ
Now that we showed the convergence of the series (16), we need an appropriate estimator for the domain of the influence to perform the surface integration. Obviously, a good candidate for estimating the P ðaÞ ðy; tÞ term would be the first term in the series (16) for a ¼ 0 since it follows from (20) that ð0Þ
ðaÞ
S 0 P S0
for a ¼ 1; 2; 3; 4:
ð21Þ
ðaÞ
ðaÞ
ðaÞ
Using (9) to express R ðy; tÞ in terms of P ðy; tÞ, we notice that the domain of influence for R ðy; tÞ can be ð1Þ estimated using the first term of the series representation for P ð1Þ ðy; tÞ, namely S0 . Taking into account ð1Þ the series representation (16), we can show that the estimator for the P ðy; tÞ term is stronger than that for P ð0Þ ðy; tÞ. Finally, the domain of influence Rs ðyÞ for the time-integrated terms P ðaÞ ðy; tÞ and RðaÞ ðy; tÞ can be estimated by the first term in the series for P ð1Þ ðy; tÞ as follows: Pe rV cos # r2 Pe Pe rs ðyÞ ¼ 2 exp : ð22Þ pr 2 4t The domain of influence is then bounded by the curve rs ðyÞ ¼ es ;
ð23Þ
where es represents a small tolerance for surface integrals, e.g., es ¼ 1030 , to ensure the contribution of the time-integrated kernels P a ðy; tÞ and RðaÞ ðy; tÞ from outside the domain Rs ðyÞ is negligible. Thus, the integration over the boundary elements (14) and (15) is performed only within the domain Rs ðyÞ associated with the source point n. No integration is needed if boundary element Cl lies outside the domain Rs ðyÞ, i.e., Cl ðxÞ \ Rs ðx nÞ ¼ 0. Notice that when the collocation point n does not lie on the boundary element Cl , the time-integrated kernels P ðaÞ ðy; tÞ and RðaÞ ðy; tÞ involve no singularity and hence can be evaluated using a standard Gaussian quadrature. However, to perform the integration in a more efficient way, we could either introduce Ns contours of pre-set levels ðiÞ rðiÞ s P es ;
i ¼ 1; 2; . . . ; Ns
and
rsðNs Þ ¼ es
ð24Þ
and segment the boundary element Cl via these contours. Or, alternatively, the contours can be evaluated directly from the outer-level contour via the following non-uniform scaling rðiÞ ð#Þ ¼ Hi rðNs Þ ð#Þ;
i ¼ 1; 2; . . . ; Ns1 ;
ð25Þ
sÞ . In the where 0 < Hi < 1, and rðNs Þ ð#Þ is the distance from the source point n to the outer-level contour rðN s former approach, difficulties arise in specifying the preset level contours to attain a desired clustering of the segments for different combinations of t and Pe as limy!0 rs ðyÞ ! 1. Thus, we utilize the latter approach throughout the study. Note that the parameter Hi is much less sensitive to the Peclet number than the preset levels in the former approach.
M.M. Grigoriev, G.F. Dargush / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4299–4312
4305
For both approaches considered above, the non-uniform clustering of the segments allows very efficient and accurate integration over the non-singular boundary element Cl influenced by the transient kernel (Fig. 2). Variable order Gaussian quadrature is utilized for every segment of interest to ensure a high level of accuracy. However, when the collocation point n lies on the straight boundary element Cl , the evaluation of the boundary integrals (14) and (15) involves a singularity as y ! 0. Analysis of terms in (9), (16), and (17) shows that the contribution due to kernel P ð1Þ ðy; tÞ vanishes in (9) since yi ni ¼ 0 for the singular boundary element, while the time-integrated kernels P ð0Þ ðy; tÞ lead to log-type singularities. All other terms result in non-singular integrals. Moreover, a boundary element integral involving P ð0Þ ðy; tÞ is singular only when the shape function Nb ðgÞ coincides with the local number of the collocation node c. Otherwise, the integral is non-singular due to the following property of the shape function 1 if c ¼ b; Nb ðgðy ðcÞ ÞÞ ¼ ð26Þ 0 if c 6¼ b: The evaluation of the non-singular integrals is similar to the algorithm described above. ð0Þ In order to remove the singularity from the integral Ub ðn; tÞ, we represent it as follows: Z Z Z P ð0Þ ðy ðbÞ ; tÞNb ðgÞ dCðxÞ ¼ P ð0Þ ðy ðbÞ ; tÞðNb ðgÞ 1Þ dCðxÞ þ P ð0Þ ðy ðbÞ ; tÞ dCðxÞ: Cl
Cl
ð27Þ
Cl
Notice that the first integral on the right-hand side of (27) is no longer singular and can be integrated using the approach outlined above. The second integral in (27) can be split into singular and non-singular parts: "Z 2 Z yV Pe r Pe ð0Þ ðbÞ i2 i Pe P ðy ; tÞ dCðxÞ ¼ e E1 dCðxÞ 4p Cl 4t Cl # k Z 2 1 X yV tk V 2 Pe r Pe i2 i Pe þ e Ekþ1 dCðxÞ : ð28Þ 4 4t k! Cl k¼1
Fig. 2. Schematics of the boundary element segmentation using scaled contours of the time-integrated kernels.
4306
M.M. Grigoriev, G.F. Dargush / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4299–4312
Fig. 3. Scheme for log-singular integration over the boundary element.
Now the log-singular integral on the right-hand-side of (28) can be represented in the following form using the notations introduced in Fig. 3: pffiffiffiffi
2 Z yV r Pe ez ð1 eÞz 4t i2 i Pe ; qk þ J ; qk : ð29Þ J e E1 dCðxÞ ¼ 4t k k Pe Cl pffiffiffiffiffiffi 2 V1 In (29), k ¼ tPe, q ¼ n1 V2 n , z ¼ LV2 Pe, and V Z a J ða; bÞ ¼ ebx E1 ðx2 Þ dx: ð30Þ 0
Unfortunately, no closed-form integral for (30) is available in the published literature except for the particular case of a ¼ b. In Appendix A we present an algorithm to numerically evaluate integral (30) for all possible cases of a > 0 and 1 6 b 6 1. Meanwhile, the non-singular integral in (28) is evaluated using the integration algorithm described earlier in this section. 3.3. Evaluation of volume integrals Although integration over the volume cells involves no singularity, it requires careful attention at small time steps for moderate and high Peclet numbers due to the nature of the kernel functions. Fig. 4 shows an upstream propagation of the convective kernel gðy; tÞ for the source located at the point of origin. Notice that the bell-shaped response is symmetric with respect to the axis of symmetry translating upstream with the speed V . When the Peclet number of the flow is relatively small, the kernel decays quickly due to large diffusion. For higher Peclet numbers, the response becomes more localized around the axis of symmetry since the convection dominates over the diffusion. Thus, it is mandatory to develop an integration algorithm that tracks the propagation of the transient convective kernel. Obviously, a straightforward integration of the kernel over volume cells would be extremely inefficient, especially for high Peclet numbers. Similar to the one-dimensional formulation, we introduce the domain of influence Rv ðyÞ for the transient convective kernel at the current time level. The contribution of the kernel gðy; tÞ is significant only within the domain of influence, whereas it is negligible outside Rv ðyÞ. The domain of influence is bounded by the contour satisfying the following gðy; tÞ ¼ ev :
ð31Þ
M.M. Grigoriev, G.F. Dargush / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4299–4312
4307
Fig. 4. Surface plots of the transient convective g-kernel for Pe ¼ 1000, V1 ¼ 1, V2 ¼ 0.
Here, ev ¼ 1030 is a small parameter to ensure the kernel is significant only within the domain of influence, and the transient convective kernel is scaled to provide 0 < gðy; tÞ 6 1 " # gðy; tÞ ðx1 n1 þ V1 tÞ2 þ ðx2 n2 þ V2 tÞ2 ¼ exp Pe : ð32Þ gðy; tÞ ¼ gðVt; tÞ 4t Introducing a frame of reference ni ¼ ni Vi t
ð33Þ
translating with the kernel, we show that contours of constant values of gðy; tÞ are concentric circles of radius rv ! 4t 1 2 rv ¼ ln : ð34Þ Pe gðy; tÞ Obviously, the center of circles coincides with the axis of symmetry. Fig. 5 shows the propagation of the circular domains of the kernel influence upstream from the source point ni . We should emphasize that the integration over a volume cell Xe is performed if any part of the volume cell is inside the domain of influence Rv ðyÞ: Xe \ Rv ðyÞ 6¼ 0: Otherwise, no integration over Xe is needed. If the volume cell is influenced by the transient convective kernel, we utilize a segmentation template shown in Fig. 6. We introduce a grid in polar coordinates ðr; hÞ so that the point of origin coincides with the axis of symmetry. Angular segmentation is uniform with the increment Dh. For the radial segmentation, we introduce Nv pre-set levels for the scaled kernel gi ðy; tÞ i ðy; tÞ < 1; eðiÞ v 6g
i ¼ 1; 2; . . . ; Nv :
ð35Þ
Here, the radius of each circular level gi ðy; tÞ is evaluated directly from (35). Note that the outer-level circle gNs ðy; tÞ is the boundary of the kernel influence. Generally, the volume cell Xe is not aligned with the segmentation grid and, thus, some of the grid segments generated in the polar coordinate system ðr; hÞ require a volume integration only over the parts
4308
M.M. Grigoriev, G.F. Dargush / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4299–4312
Fig. 5. Circular domains of kernel influence propagate upstream from the source point.
Fig. 6. Schematics of the volume cell segmentation influenced by the kernel. Circular domain of influence is represented by polygons.
that belong to the volume cell. For the sake of simplicity in evaluating the intersections between influence domain and volume cells, the circular domains are replaced with polygons as shown in Fig. 6. For relatively small increments Dh, the accuracy of the integration is not impacted by this approximation. Note that the volume integrals are evaluated using a standard Gaussian quadrature. Similar to the boundary element integration, we utilize a variable order of the Gaussian quadrature for every volume segment of interest to ensure a high level of accuracy.
M.M. Grigoriev, G.F. Dargush / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4299–4312
4309
4. Assembly The assembly procedure for the BEM approach is standard with no additional complications arising from the use of the higher-order time functions. In contrast to assembly procedure for time-dependent heat diffusion presented by Grigoriev and Dargush [11], the transient convective diffusion leads to a sparse global matrix for both one- and two-dimensional formulations. Since the global matrix for one-dimensional problems is very compact, the choice of the solution algorithm has not been a critical issue there. However, the size of the global matrix for the two dimensions can be as large as 100,000 for the problems considered in this study. For the problems of practical importance, the number of degrees of freedom could easily go beyond a million. We should emphasize that the sparsity of the global matrix increases dramatically with an increase of the Peclet number. The non-zero elements of the matrix are clustered around the main diagonal. Off-diagonal elements decay rapidly or are virtually zero for the entire matrix. Therefore, iterative algorithms would be of the primary choice to solve the global set of equations. In this paper, we utilize the preconditioned biconjugate gradient method for sparse linear systems presented by Press et al. [12]. We should note that the use of the main diagonal as a preconditioner results in a stable, accurate and efficient solution algorithm even for very high Peclet numbers.
5. Conclusions The higher-order boundary element methods for time-dependent convective diffusion phenomena presented in Part I of the paper are extended to the two-dimensional problems. First, the transient convective kernel gðy; tÞ is represented as an infinite series. Although the series converges for all time step sizes, the applicability of the numerical algorithm is constrained by inequality (1). We emphasize that the critical time step size Dtc following from (1) is not subject to the mesh discretization contrary to the well-known constraint on local mesh-based Peclet number Pec 6 2 when no upwinding is utilized. Note that the convective kernels used in this paper provide an analytical upwinding for any Peclet number. To relieve the restriction or, rather, recommendation on the time step size of the algorithm, another form of the series representation for the kernels is needed. However, the implementation of a suitable alternative goes beyond the scope of this paper. Here the series of the convective kernel is integrated term-by-term over finite time intervals for all time interpolation functions. To facilitate accurate evaluation of the discrete integral equation coefficients over both singular and non-singular boundary elements, we have developed and implemented an adaptive algorithm that accounts for the complex nature of the time-dependent convective kernels. The algorithm tracks the domains influenced by the kernels entering the surface and volume integrals. Accordingly, it provides an automatic segmentation of the influence domains to allow efficient and accurate numerical integrations of the kernels over higher-order boundary elements and volume cells. The numerical results for the transient convective heat diffusion are presented and discussed in Part III.
Acknowledgements The work described in this paper was partially supported by the National Science Foundation under Grant CMS 9700387. The authors gratefully acknowledge this support. Furthermore, the authors are grateful to the reviewer for many valuable comments that helped improve the paper.
4310
M.M. Grigoriev, G.F. Dargush / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4299–4312
Appendix A. Evaluation of integral J(a; b) A.1. Small parameter a 6 1 Here we present an algorithm to numerically evaluate integral (30) given by Z a J ða; bÞ ¼ ebx E1 ðx2 Þ dx:
ðA:1Þ
0
Note that an accurate evaluation of integral (A.1) is required for any combinations of a > 0 and 1 6 b 6 1. Also, notice that the exponential integral E1 ðx2 Þ leads to a log-type singularity as x ! 0. We should emphasize that the straightforward substitution of the following series representation for the exponential integral [3] E1 ðx2 Þ ¼ c 2 lnðxÞ
1 X ð1Þk x2k kk! k¼1
ðA:2Þ
into (A.1) and a subsequent integration allows us to obtain 1 X c ð1 expðabÞÞ 2 ð1Þ B2k ða; bÞ þ ½c þ ln jbj þ expðabÞ lnðaÞ þ E1 ðabÞ : b b kk! k¼1 k
J ða; bÞ ¼
ðA:3Þ In (A.2) and (A.3), c ¼ 0:57721566. . . is the EulerÕs constant [4], and the integral Z a Bk ða; bÞ ¼ ebx xk dx
ðA:4Þ
0
can be expressed by the following recurrence formula Bk ða; bÞ ¼
1 ½kBk1 ða; bÞ ak expðabÞ b
B2 ða; bÞ ¼
2 expðabÞð2 þ 2ab þ a2 b2 Þ b3
for k > 2:
ðA:5Þ
Here, for k ¼ 2:
Since Bk ða; bÞ are needed only for even k, we modify (A.5) to this form
1 kðk 1Þ k Bk2 ða; bÞ ak1 a þ Bk ða; bÞ ¼ expð abÞ for k P 4: b b b
ðA:6Þ
ðA:7Þ
Unfortunately, the recurrence evaluation of (A.7) is stable only for small parameters a 6 1. For higher a, the use of (A.7) to evaluate the integral J ða; bÞ results in large errors. Thus, this algorithm is limited to a 6 1 only.
A.2. Large parameter 1 < a 6 10 Let us introduce parameter k and recast integral (A.1) to the following form Z a J ða; b; kÞ ¼ expðbxÞE1 ðk2 x2 Þ dx: 0
ðA:8Þ
M.M. Grigoriev, G.F. Dargush / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4299–4312
4311
Differentiation of (A.8) with respect to k gives Z a o 2 J ða; b; kÞ ¼ expðbx k2 x2 Þ dx: ok k 0 Integration of (A.9) is straightforward and leads to: 2 pffiffiffi b p exp o b b 4k2 J ða; b; kÞ ¼ erf erf þ ak : ok 2k 2k k2 We notice that Z 1 o J ða; b; kÞ dk ¼ J ða; bÞ ok 1
ðA:9Þ
ðA:10Þ
ðA:11Þ
since lim J ða; b; kÞ ¼ 0
ðA:12Þ
J ða; b; 1Þ ¼ J ða; bÞ:
ðA:13Þ
k!1
and Integration of (A.10) over the parameter k and the use of (A.11) allow us to recast the desired integral in the following form: # pffiffiffi " Z b 2 2 p ab 2 J ða; bÞ ¼ expðy Þ erf þ y erfðyÞ dy : ðA:14Þ b 2y 0 We should note that the behavior of the integrand in (A.14) is no longer singular. Thus, the Gaussian quadrature can be applied to numerically evaluate the integral in (A.14). Notice that the integral (A.14) can also be used for small 0 < a < 1. For the limiting case of b ! 0, the closed form integration of (A.1) gives: pffiffiffi J ða; 0Þ ¼ aE1 ða2 Þ þ p erfðaÞ: ðA:15Þ A.3. Very large parameter a > 10 For large parameters a > 10, integral (A.1) can be represented in the following form: Z a J ða; bÞ ¼ J ð10; bÞ þ ebx E1 ðx2 Þ dx:
ðA:16Þ
10
It is possible to show that the value of the second integral in (A.16) is negligible for all 1 6 b 6 1. Thus, only the contribution of the first integral J ð10; bÞ is considered for large a.
References [1] M.M. Grigoriev, G.F. Dargush, Boundary element methods for transient convective diffusion. Part I: General formulation and 1D implementation, Comput. Methods Appl. Mech. Engrg. 192 (2003) 4281–4298. [2] A.P. Prudnikov, Yu.A. Brychkov, O.I. Marichev, Integrals and Series, vol. 1: Elementary Functions, Gordon and Breach Science Publishers, New York, 1986.
4312
M.M. Grigoriev, G.F. Dargush / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4299–4312
[3] I.S. Gradshtein, I.M. Ryzhik, Tables of Integrals, Series and Products, Academic Press, 1980. [4] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, Dover, New York, 1964. [5] W. Gautschi, Algorithm 471: Exponential integrals, Commun. ACM 16 (1973) 761–763. [6] C.A. Brebbia, S. Walker, Boundary Element Techniques in Engineering, Newness-Butterworths, London, 1980. [7] P.K. Banerjee, R. Butterfield, Boundary Element Methods in Engineering Science, McGraw-Hill, London, 1981. [8] C.A. Brebbia, J.C.F. Telles, L.C. Wrobel, Boundary Element Techniques, Springer, Berlin, 1984. [9] P.K. Banerjee, The Boundary Element Methods in Engineering, McGraw-Hill, London, New York, 1994. [10] K.J. Bathe, Finite Element Procedures, Prentice Halls, Englewood Cliffs, 1996. [11] M.M. Grigoriev, G.F. Dargush, Higher-order boundary element methods for transient diffusion problems. Part I: Bounded flux formulation, Int. J. Num. Methods Engrg. 55 (2002) 1–40. [12] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical recipes in C. The Art of scientific computing, Cambridge University Press, Cambridge, 1992.