- Email: [email protected]
Localized method of approximate particular solutions with Cole–Hopf transformation for multi-dimensional Burgers equations
Engineering Analysis with Boundary Elements 40 (2014) 78–92
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Localized method of approximate particular solutions with Cole–Hopf transformation for multi-dimensional Burgers equations C.Y. Lin a,b, M.H. Gu c, D.L. Young a,b,n, C.S. Chen d,e a
Department of Civil Engineering, National Taiwan University, Taipei 10617, Taiwan Hydrotech Research Institute, National Taiwan University, Taipei 10617, Taiwan c National Center for Research on Earthquake Engineering, Taipei 10668, Taiwan d Department of Mathematics, The University of Southern Mississippi, Hattiesburg, MS, USA e Department of Engineering Mechanics, Hohai University, Nanjing, Jiangsu 210098, China b
art ic l e i nf o
a b s t r a c t
Article history: Received 23 August 2013 Accepted 28 November 2013 Available online 25 December 2013
The Burgers equations depict propagating wave with quadratic nonlinearity, it can be used to describe nonlinear wave propagation and shock wave, where the nonlinear characteristics cause difficulties for numerical analysis. Although the solution approximation can be executed through iterative methods, direct methods with finite sequence of operation in time can solve the nonlinearity more efficiently. The resolution for nonlinearity of Burgers equations can be resolved by the Cole–Hopf transformation. This article applies the Cole–Hopf transformation to transform the system of Burgers equations into a partial differential equation satisfying the diffusion equation, and uses a combination of finite difference and the localized method of approximate particular solution (FD-LMAPS) for temporal and spatial discretization, respectively. The Burgers equations with behaviors of propagating wave, diffusive N-wave or within multi-dimensional irregular domain have been verified in this paper. Effectiveness of the FD-LMAPS has also been further examined in some experiments, and all the numerical solutions prove that the FD-LMAPS is a promising numerical tool for solving the multi-dimensional Burgers equations. & 2013 Elsevier Ltd. All rights reserved.
Keywords: Burgers equations Cole–Hopf transformation Diffusion equation Meshless methods Irregular domain Localized method of approximate particular solutions
1. Introduction The Burgers equations, are quasi-linear partial differential equations (PDEs) consisting an unsteady term, a nonlinear convection term and a diffusion term. The Burgers equations depict propagating wave with quadratic nonlinearity, it can be used to describe nonlinear wave propagation and shock wave [1]. The Burgers equations have been widely applied for mathematic analysis and problems of fluid mechanics. In this paper we approximate the solutions of the Burgers equations by the localized method of approximate particular solutions (LMAPS) [2] which is a novel developed numerical scheme using radial basis functions (RBFs). In order to overcome the nonlinearity of the Burgers equations, the implementation of the LMAPS to solve the resulting diffusion equation of applying Cole–Hopf transformation is considered in this paper. Basic solution approximation for nonlinear equations usually requires iterative approaches such as the Newton–Ralphson or Picard method [3,4]. However, the convergence of iterative approaches is not
n Corresponding author at: Department of Civil Engineering, National Taiwan University, Taipei 10617, Taiwan. Tel./fax: þ 886 2 2362 6114. E-mail address: [email protected] (D.L. Young).
0955-7997/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.enganabound.2013.11.019
always guaranteed. In contrast to iterative approaches, direct approaches for approximating the solution of the Burgers equations can overcome the nonlinearity through finite time sequence which avoids the issue of convergence. Such direct approaches can be easily found in the literature for solving Burgers equations, for example the well-known Eulerian–Lagrangian method [5], or the direct approach that will be applied in this paper — the Cole–Hopf transformation [6,7]. The Cole–Hopf transformation was introduced by Hopf and Cole [6,7] to deal with the nonlinearity of the viscous Burgers equations. After applying the Cole–Hopf transformation, the system of Burgers equations can be transformed into a single variable PDE satisfying the diffusion equation. The analytical solution of the diffusion equation can be simply found by solving the Cauchy problem in regular domain. Therefore, it provides a method to obtain analytical solutions for the Burgers equations. Although the Cole–Hopf transformation was originally developed for deriving analytical solutions, in this research the Cole–Hopf transformation is only applied to deal with the nonlinearity of Burgers equations for further numerical processes. For the numerical implementation, it is necessary for the essential conditions, i.e. initial and boundary conditions, to be transformed according to the definition of the Cole-Hopf transformation. This implies both the governing equation and the essential
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92
conditions are transformed into the single variable form. Although the Burgers equations have been transformed into single variable form, the number of equations that depict the essential conditions remains the same, which indicates the linear system generated for dealing with multi-dimensional problems will be overdetermined. Hereby a numerical linear solver named MA49 using least-squares QR method [8] with capability for dealing with the sparse over-determined system is applied in this article. After overcoming the nonlinearity by the Cole–Hopf transformation, further solution approximation is required. For numerical analysis, there are well-known and widely applied mesh-dependent numerical methods, such as the finite difference method (FDM) [3,9,10], the finite element method (FEM) [11,12] or the boundary element method (BEM) [13,14]. All these classical numerical methods are popular and effective but require difficult procedures of mesh generation or numerical quadrature. In order to avoid such difficulties, various meshless numerical methods such as the smoothed particle hydrodynamics (SPH) [15], the multiquadrics collocation method (MQ) [16–20], the method of fundamental solutions (MFS) [5,21,22], the method of particular solutions (MPS) [23,24] and the meshless local Petrov–Galerkin method (MLPG) [25,26], to name just a few, have been developed. For global domain-type meshless methods with RBFs, the resultant interpolation matrix of MQ or the MPS is often dense and ill-conditioned which makes the computation inefficient and unstable. These disadvantages cause the global type numerical methods to be restrained for large-scale computational problems. Hereby extension researches on localized formulation of global domain-type meshless numerical methods have been intensively studied among recent literature [2,27–31]. In this paper we focus on the localized formulation to modify the method of approximate particular solutions (MAPS), also called the localized method of approximate particular solutions (LMAPS). The main concept of localization was introduced by Yao, Chen, and Kolibal [2]. The LMAPS is based on the global MAPS with modification of choosing limited number of weighting points within the entire local domain. The breakthrough of the localization technique allows dense matrix of the original global linear system to reduce into a sparse linear system. This will decrease the risks of ill-conditioned linear system and lessen computational memory loading. Therefore the LMAPS as a meshless method has the advantages without requiring any mesh generation or numerical quadrature, and has also avoided the drawbacks of the global MAPS or MPS. The processes of localization increase the efficiency for solving linear systems, meanwhile maintaining appreciable accuracy of solutions. There is a very similar localized meshless numerical method called the local radial basis function differential quadrature method (LRBFDQ) [28,31]. In numerical algorithm, these two localized meshless numerical methods share many common features; however the basis function chosen has made them mathematically different. For the basis function approximating the solutions, LRBFDQ applies the differential quadrature (DQ), while the LMAPS uses the integrated radial basis functions (integrated RBF). Mathematically, integrating a function will obtain a solution function with an unknown constant, causing differences in solutions obtained from these two methods. Despite the differences, both schemes show very good performance in numerical simulation due to the same localization technique they shared. The technique itself will be further described in this paper and Ref. [2], and the comparison of related methods can be found in Ref. [32]. Finally the time domain and space domain are considered with different approaches of discretization. The time domain is discretized by the Crank–Nicolson scheme of FDM, where the Crank– Nicolson scheme is an implicit unconditional stable scheme with second-order of accuracy. For the space domain, the LMAPS is
79
chosen to discretize the space domain and approximate the solutions. We summarize here the techniques applied for this study: the Cole–Hopf transformation transforms the Burgers equations into a PDE satisfying diffusion equation, the Crank– Nicolson scheme approximates transient term and the LMAPS applied for space discretization. Combining these methods and integrating into the proposed numerical model of FD-LMAPS with the Cole–Hopf transformation, we will carry out numerical experiments to solve the multi-dimensional Burgers equations. There are four numerical experiments in this paper which are named a one-dimensional Burgers problem with propagating wave, a one-dimensional Burgers problem with diffusive N-wave, a two-dimensional Burgers problem in L-shaped domain, and a three-dimensional Burgers problem in an irregular domain. Further details of the proposed numerical model are provided in the following sections.
2. Governing equations and essential conditions 2.1. Governing equations The Burgers equations are quasi-linear PDEs which include nonlinear convection term and diffusion term. For the global domain Ω including the boundary ∂Ω, the Burgers equations in vector form listed as ,
, , , ∂U 1 þ ðU U∇ÞU ¼ ∇2 U ; Re ∂t , ,
ð1Þ ,
where U ðx ; tÞ ¼ ðu1 ; ⋯; ud Þ is the vector of velocity, x ¼ ðx1 ; ⋯; xd Þ A Rd , d ¼ 1; 2; 3, and Re is the Reynolds number. 2.2. Initial conditions and boundary conditions The initial conditions and boundary conditions are the essential conditions for obtaining unique solutions for the governing equations. Only Dirichlet-type boundary conditions are considered in this study, the reason will be explained in Section 3.2. The initial conditions in general form is depicted as , , ,, U ðx ; tÞ ¼ s ðx Þ; ð2Þ t¼0
and the boundary conditions is ,, , , ¼ β ðx ; tÞ: U ðx ; tÞ
ð3Þ
,
x A ∂Ω ,,
,,
Here s ðx Þ ¼ ðs1 ; ⋯; sd Þ is the given initial ,velocity, β ðx ; tÞ ¼ ðβ 1 ; ⋯; β d Þ is the given boundary conditions of U . 3. Theoretical model The Cole–Hopf transformation for the governing equations and essential conditions are introduced with some detailed derivations. The demonstrations, the applications, and restrictions of the Cole–Hopf transformation applied in this research will be followed. 3.1. Cole–Hopf transformation for the Burgers equations In order to transform the system of Burgers equations into diffusion equations, the Cole–Hopf transformation is introduced as follows [6,7]: ,
Uqþ
2 ∇q ¼ 0 and q a 0; Re ,
ð4Þ
where qðx ; tÞ is the undetermined variable after the transformation.
80
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92
Substituting Eq. (4) into the system of Burgers equations in Eq. (1), we have ∇q ∂q 1 2 ∂q 1 2 ∇ q ∇ ∇ q ¼ 0; ð5Þ q ∂t Re ∂t Re
straightforwardly, integrated basis functions corresponding to the Laplace operator are chosen to describe the particular solutions of Eq. (9).
the solution of following diffusion equation automatically satisfies Eq. (5):
4.2. Localized method of particular solutions
∂q 1 2 ∇ q ¼ 0: ∂t Re
In this paper, the LMAPS is applied to discretize the space , domain of the diffusion equation. Let x i A Ω; i ¼ 1; ⋯; N be the global points, where N is the total number of global points. From the MPS, q is approximated by a series of integrated basis functions of F,
ð6Þ
If the problem is well-posed, then the unique solution of Eq. (6) exists. Detailed derivation for the Cole–Hopf transformation in multidimensional case is given in the Appendix A.
N
,
,
j¼1
The essential conditions include both initial conditions and boundary conditions in this article. By the definition of the Cole–Hopf transformation as Eq. (4), the essential conditions in Eqs. (2) and (3) can be derived as follows: 2 , , ,, Initial condition : s ðx Þqðx ; tÞ þ ∇qðx ; tÞ ¼ 0; ð7Þ Re t¼0 , 2 , , , Boundary condition : β ðx ; tÞqðx ; tÞ þ ∇qðx ; tÞ ¼ 0: , Re x A ∂Ω
,
qðx ; tÞ ffi ∑ αj Fð‖x x j ‖Þ;
3.2. Cole–Hopf transformation for essential conditions
,
While providing Dirichlet type conditions of U , the essential conditions of Eqs. (7) and (8) become Robin type conditions of q. Therefore the proposed numerical scheme can discretize the transformed essential conditions, of q directly. However if the Neumann type conditions of U were provided for essential conditions, the transformed essential conditions will become nonlinear. The nonlinearity should be resolved for further analysis or other approximations. In this paper, only the Dirichlet type conditions will be considered. Therefore the governing equation and essential conditions to be solved are described by Eqs. (6)–(8).
4. Numerical methods
x A Ω:
ð10Þ
where ∇2 F ¼ f , f is the basis function and αj are the weighting coefficients to be determined. For the LMAPS, we can create a local influence domain ωi for , , each global point x i . The local points are defined as x i;k A ωi ; k ¼ 1; ⋯; NL (see Fig. 1). By the MPS, on each local influence domain ωi , we have NL
,
ð8Þ
,
,
,
qðx ; tÞ ffi ∑ αi;m Fð‖x x i;m ‖Þ; m¼1
,
x A ωi ;
ð11Þ
where αi;m is the weighting coefficients of the mth local point in the local influence domain ωi . By collocation method on Eq. (11), qi ¼ Fi αi ;
ð12Þ
T , , , , ; Fi ¼ Fð‖x i;k x i;m ‖Þ ; where qi ¼ q x i;1 ; t ; ⋯; q x i;NL ; t NLNL αi ¼ αi;1 ; ⋯; αi;NL T . If Fi is invertible, Eq. (12) can be rewritten as
αi ¼ Fi 1 qi :
ð13Þ
Now the Laplace operator ∇2 can be applied to Eq. (11) to obtain the following expression, ,
,
NL
,
,
∇2 qðx ; tÞ ffi ∇2 qðx i ; tÞ ¼ ∑ αi;m f ð‖x i x i;m ‖Þ ¼ f i αi :
ð14Þ
m¼1
This section introduces the implementation of numerical methods, the Crank–Nicolson scheme is used for discretization of the time domain and the LMAPS is applied for discretization of the space domain.
, , , , where f i ¼ f ð‖x i x i;1 ‖Þ; ⋯; f ð‖x i x i;NL ‖Þ . From Eqs. (13) and (14), we have
4.1. Crank–Nicolson scheme
where φi ¼ f i Fi 1 . In order to transfer the above localsystem into global system, it requires to map the local vector φi ¼ φi;1 ; ⋯; φi;NL into the global
The time-marching Crank–Nicolson scheme [10] is applied to discretize the time domain for numerical experiments in this study. The second-order accuracy and unconditional stable Crank–Nicolson scheme is selected for the time domain of diffusion equation, which can be written as follows: qn þ 1
Δt 2 n þ 1 Δt 2 n ∇ q ∇ q þ OðΔt 2 Þ; ¼ qn þ 2Re 2Re
,
∇2 qðx i ; tÞ ¼ f i Fi 1 qi ¼ φi qi ;
ð15Þ
ð9Þ
where Δt is the time interval, OðΔt 2 Þ is the truncation error, superscript of q denotes the number of time step and n marks the time step t n ¼ nΔt. Eq. (9) can be considered as a modified Helmholtz equation with time-dependent source term, the particular solutions of Eq. (9) can be described by integrated basis functions corresponding to either Laplace operator or modified Helmholtz operator. For particular solutions of modified Helmholtz equation, the modified Bessel functions make it difficult to obtain suitable RBF for applications. Since the integrated RBFs corresponding to the Laplace operator do not share the same problem and can be implemented to various PDEs other than Helmholtz equation more
,
Fig. 1. The diagram of global domain Ω, local influence domain ωi of point x i , global , , points x and local point x i;k .
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92
vector Φi where
Φi ¼ Φi;1 ; ⋯; Φi;N ; Φi;j ¼
8 > < 0;
, xj2 =
> : φi;k ;
, xj
ωi ,
¼ x i;k A ωi
:
ð16Þ
Therefore, the system extends into the global domain Ω, ∇2 q ¼ Φq;
where Φ ¼ Φ1 ; ⋯; ΦN
T
T , , , and q ¼ qðx 1 ; tÞ; ⋯; qðx N ; tÞ .
ð17Þ
Applying the similar procedures to the gradient operator ∇, by defining g ¼ ∇F, we have ,
,
NL
,
,
∇qðx ; tÞ ffi ∇qðx i ; tÞ ¼ ∑ αi;k gð‖x i x i;m ‖Þ ¼ gi αi ; m¼1
ð18Þ
,
ð19Þ ∇qðx i ; tÞ ¼ gi Fi 1 qi ¼ γi qi ; , , , , and γi ¼ gi Fi 1 . The where gi ¼ g ‖x i x i;1 ‖ ; ⋯; g ‖x i x i;NL ‖ local system can be transferred into global system by mapping the local vector γi ¼ γ i;1 ; ⋯; γ i;NL into the global vector 8 , > < 0; xj 2 = ωi ; Γi ¼ Γ i;1 ; ⋯; Γ i;N ; Γ i;j ¼ ð20Þ , , > : γ i;k ; x j ¼ x i;k A ωi ; ∇q ffi Γq:
ð21Þ
81
T Note that Γ ¼ Γ1 ; ⋯; ΓN and sub-elements Γi ; i ¼ 1; ⋯; N are defined as Eq. (20). Further details of the global extension of φi and γi can be found in Ref. [2]. Each linear term of the governing equations or the essential conditions can be discretized by the LMPS through similar processes. Further detailed numerical procedures for solving Burgers equations will be given in Section 5.
4.3. Integrated multiquadrics radial basis functions The integrated multiquadrics radial basis functions (integrated MQ-RBFs) are used for the weighting function in this work. The MQ-RBF and the integrated MQ-RBFs are defined as follows [32–34]. The MQ-RBF: f ðr Þ ¼
pffiffiffiffiffiffiffiffiffiffiffiffiffiffi r 2 þ c2 :
ð22Þ
The 1-dimensional integrated MQ-RBF: FðrÞ ¼
pffiffiffiffiffiffiffiffiffiffiffiffiffiffi i pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1h
2c2 þ r 2 r 2 þ c2 þ r ; r 2 þ c2 þ 3c2 r ln 6
The 2-dimensional integrated MQ-RBF: pffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 1 r 2 þ c2 þ c ; FðrÞ ¼ ðr 2 þ 4c2 Þ r 2 þ c2 c3 ln 9 3
ð23Þ
ð24Þ
Fig. 2. The evolution of u1 for Example 6.1 at different time step with (a) Re ¼ 1, Δt ¼ 10 2 , 21 global points and 3 local points, (b) Re ¼ 10, Δt ¼ 10 2 , 51 global points and 3 local points, (c) Re ¼ 100, Δt ¼ 10 2 , 501 global points and 3 local points.
82
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92
The 3-dimensional integrated MQ-RBF: pffiffiffiffiffiffiffiffiffiffi 8 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4 2 2 > 1 c < 24 ð2r 2 þ5c2 Þ r 2 þ c2 þ 8r ln r þ cr þ c ; FðrÞ ¼ > : c3 ; 3
5. Numerical procedures r 4 0;
ð25Þ
r ¼ 0;
where r is considered as the distance between the collocation points and constant c is the shape parameter of MQ and taken in the order of 1.0 in this study. There are various approaches proposed for obtaining optimal shape parameters [35–37], each claiming to be optimal and consistent, suitable for applications of MQ-RBF.
In this section the detailed procedures for solving the Burgers equations are provided by using the FD-LMAPS and the Cole–Hopf transformation. 5.1. Initial conditions ,
From Section 3.2, the initial conditions of velocity field U have been derived to Eq. (7) as the Robin type condition. The LMAPS and integrated MQ-RBF were applied in each term for discretiza, tion. For each global point x i , the local influence domain ωi can be
Fig. 3. The error evolution of Example 6.1 of u1 with different number of global points, (a) EL2 of Re ¼ 1, (b) ERel of Re ¼ 1, (c) EL2 ofRe ¼ 10, (d) ERel ofRe ¼ 10, (e) EL2 of Re ¼ 100, and (f) ERel of Re¼ 100.
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92
83
ð26Þ
dN equations for describing the initial conditions, where dN ¼ d U N. Hereby the global linear system is described as follows 2 Is þ Γs q ¼ 0; ð27Þ Re dNN
where qi ¼ qðx i ; tÞ. The above local system can be extended to global system. For d-dimensional problems, there are d initial conditions for each , global point of velocity field U . Therefore, N global points implies
where Is and Γs are defined as follows. For d-dimensional cases: 2 3 2 3 I1 Γ 6⋮7 6 7 ; Γs ¼ 4 ⋮ 5 : ð28Þ Is ¼ 4 5 Id Γ dNN
established as mentioned in Section 4.2 to obtain the following local linear system: ,,
s ðx Þq þ
2 2 , ∇q ¼ 0 ) sd ðx i Þqi þ γi qi ¼ 0; Re Re ,
dNN
Fig. 4. The evolution of u1 for Example 6.2 at different time step with (a) Re ¼ 1, Δt ¼ 10 3 , 501 global points and 3 local points, t A 10 2 ; 1 , and (b) Re ¼ 100, Δt ¼ 10 4 , 2 501 global points and 3 local points, tA (10 ,1).
Fig. 5. The error evolution of Example 6.2 of u1 with different number of global points, (a) EL2 of Re ¼ 1, t A 10 2 ; 10 , (b) ERel of Re ¼ 1, t A 10 2 ; 10 , (c) EL2 of Re ¼ 100,
2 2 t A 10 ; 1 , and (d) ERel of Re ¼ 100, t A (10 ,1).
84
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92 , xi.
Similar to the numerical process of the initial conditions, the , local influence domain is established for each boundary point x ib to obtain the following local linear system from Eq. (8),
,,
β ðx ; tÞq þ
2 2 , ∇q ¼ 0 ) β d ðx ib ; tÞqi þ γi qi ¼ 0: Re Re
ð32Þ
The above local system can also be extended to global system [2]. For d-dimensional problems, each boundary point provides d , boundary conditions of velocity field U , implies NB boundary points will result in d U NB equations. The global linear system of boundary conditions can be described as follows: 2 Ib þ Γb q ¼ 0; ð33Þ Re dNBN
Fig. 6. The L-shaped domain Ω for Example 6.3, a total number of 1027 points are shown in this figure. The black border is the boundary ∂Ω, the black squares are the boundary points and the red circles depict the global points. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Note that 2 3 Id;1 6 ⋮ 7 Id ¼ 4 5; Id;N
2
, 6 7 Id;i ¼ 40; ⋯; sd ðx i ; tÞ ; ⋯; 05 |fflfflfflffl{zfflfflfflffl}
:
ð39Þ
1N
After applying the Cole–Hopf transformation, there is a single variable q for each global point. It implies N unknowns of q needs , to be approximated for satisfying the velocity field U . The homogeneous solution of over-determined sparse linear system can be solved by the singular value decomposition of MA49 [8]. 5.2. Governing equations
It is necessary to introduce the boundary conditions into the system combined with governing equation to satisfy the velocity field. Therefore, the two linear systems are combined as follows # " " # Δt Φ I Δt Φ I þ 2Re nþ1 2Re q ¼ qn : ð36Þ 2 Γb Ib þ Re ½0 ðN þ dNBÞN
ðN þ dNBÞN
By using least-squares QR method [8], the numerical solver named MA49 is applied to obtain the solution of the sparse overdetermined linear system. Algorithm.
For space discretization of Eq. (9), the LMAPS scheme is , implemented. Let x i be the interior points and by collocation method, the following local system can be obtained.
Δt 2 n þ 1 Δt 2 n ∇ q ∇ q ¼ qn þ 2Re 2Re Δ t Δt ) qni þ 1 φ qn þ 1 ¼ qni þ φ qn : 2Re i i 2Re i i
where for each boundary point with a corresponding global point, and 2 3 2 3 Jd;1 , 6 7 6 7 Jd ¼ 4 ⋮ 5; Jd;ib ¼ 40; ⋯; β d ðx ib ; tÞ ; ⋯; 05 : ð35Þ |fflfflfflfflffl{zfflfflfflfflffl} Jd;NB i 1N
3
i
where Ib and Γb are defined as follows. For d-dimensional cases: 2 3 2 3T J1 6 7 6 7 ð34Þ Ib ¼ 4 ⋮ 5; Γb ¼ 4Γ1 ; ⋯; ΓNB ; ⋯; Γ1 ; ⋯; ΓNB 5 ; |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Jd Cycledtimes
qn þ 1
ð30Þ
The above local system can be assembled into a global system. It follows that Δt Δt Φ qn þ 1 ¼ Iþ Φ qn ; ð31Þ I 2Re 2Re NN NN where matrix I is the N N identity matrix. Notice that the terms on the right hand side are all known from the previous time step. Since the governing equation after the Cole–Hopf transformation consist only one equation for each point in the space domain, this sparse linear system is considered to be fully-determined. However it is still necessary to introduce the boundary conditions to satisfy the velocity field. 5.3. Boundary conditions ,
The Dirichlet type boundary conditions of velocity field U in Eq. (3) can be derived to the Robin type boundary condition as shown in Eq. (8). Each term is derived with the technique of the LMAPS and use integrated MQ-RBF for basis function. Let NB denotes the number of boundary points. For each boundary point , x ib A ∂Ω; ib ¼ 1; ⋯; NB, there exists a corresponding global point
Step 1: Construct and solve the linear system of the initial conditions as Eq. (27). Step 2 Construct the linear system of the governing equations and boundary conditions as in Eq. (36) and obtain solution of qn þ 1 . ^ nþ1 , according to the Step 3: Resolve numerical solution U definition of Cole–Hopf transformation. Step 4: Use solution qn þ 1 to replace qn , and repeat Step 2 and 3 for each time steps.
6. Numerical experiments In order to display the capacity of the proposed numerical model for solving the Burgers equation, four numerical experiments are considered. All numerical experiments are compared with analytical solutions for error analysis, and the essential conditions can be obtained by substituting the time and position of boundary points into the analytical solution. The error distribution applied in this research includes the maximum relative error ERel and the L2 norm error EL2 , which are defined as , , , ^ , Max U x i ; t U x i ; t 1 r i r N ; ð37Þ ERel ¼ , , Max U x ; t i 1rirN
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u ^ 2 u N , , , , u∑ U x i ; t U x i ; t u ui ¼ 1 u : EL2 ¼ u 2 N , t , U x ; t ∑ i
ð38Þ
i¼1
^ , Note that U is the analytical solution, U is the numerical solution and N is the number of test points. ,
85
number defines the dominancy between convection and diffusion. Higher Reynolds number leads to a stronger effect of the convection, and vice versa. In this example, when higher Reynolds number is adopted, the solution will be more similar to a step function or a Heaviside function propagating through the domain, making it more difficult to obtain a good numerical solution. The analytical solution of this experiment is given as [38]: 1
; x1 A ½0; 1; t A ð0; 2: 1 þ exp Reð2x1 t Þ=4
6.1. One-dimensional problem with propagating wave
u1 ðx1 ; tÞ ¼
The first example is set to verify the capability of proposed numerical model depicting a propagating wave. Here the Reynolds
Fig. 2(a–c) is the evolution of velocity u1 at different time steps. From Fig. 2(a) we can see when Re ¼ 1, the velocity u1 distributes
1
t=0.0000
u1 950.0 900.0 850.0 800.0 750.0 700.0 650.0 600.0 550.0 500.0 450.0 400.0 350.0 300.0 250.0 200.0 150.0 100.0 50.0
0.8
0.6
x2 0.4
0.2
0
1
0
0.2
0.4
x1
0.6
0.8
t=0.0000
0.6
x2 0.4
0.2
1
0
0.2
0.4
x1
0.6
0.8
0.6
x2 0.4
0.2
0
0.2
0.2
0
1
0.4
x1
0.6
0.8
1
0
0.2
0.4
x1
0.6
0.8
1
t=0.0500
u1 9.5 9.0 8.5 8.0 7.5 7.0 6.5 6.0 5.5 5.0 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0
0.8
0.6
x2 0.4
0.2
u1
0.8
0
0.4
0
0.095 0.090 0.085 0.080 0.075 0.070 0.065 0.060 0.055 0.050 0.045 0.040 0.035 0.030 0.025 0.020 0.015 0.010 0.005
950.0 900.0 850.0 800.0 750.0 700.0 650.0 600.0 550.0 500.0 450.0 400.0 350.0 300.0 250.0 200.0 150.0 100.0 50.0
x2
1
t=0.0000
u1
0.6
u1 9.5 9.0 8.5 8.0 7.5 7.0 6.5 6.0 5.5 5.0 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0
t=0.0005
0.8
1
0.8
0
1
1
0
0.2
0.4
x1
0.6
0.8
1
u1
t=5.0000
0.095 0.090 0.085 0.080 0.075 0.070 0.065 0.060 0.055 0.050 0.045 0.040 0.035 0.030 0.025 0.020 0.015 0.010 0.005
0.8
0.6
x2
0.4
0.2
0
ð39Þ
0
0.2
0.4
x1
0.6
0.8
1
Fig. 7. The evolution of u1 for Example 6.3, (a) Re ¼ 0:01, t ¼ 0, (b) Re ¼ 0:01, t ¼ 0:0005, (c) Re ¼ 1, t ¼ 0, (d) Re ¼ 1, t ¼ 0:05, (e) Re ¼ 100, t ¼ 0, and (f) Re ¼ 100, t ¼5.
86
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92
as a nearly straight line. When Reynolds number increases to 10 or 100 as shown in Fig. 2(b and c), an obvious wave front can be observed propagating through the computational domain. Adopting higher Reynolds number results in steeper wave fronts. Fig. 3(a–f) demonstrates the error evolution of different Reynolds number and error definition. Both ERel and EL2 yield acceptable accuracy, and better accuracy can be achieved by adding more global points. The comparison between ERel and EL2 has also shown that the distributions of the two errors have similar trend. The similarity in ERel and EL2 means the relative errors are evenly distributed throughout the domain, including the steep wave front of physical value u1 . This is an indication that our proposed numerical model is capable for solving steep wave front problems.
1
t=0.0000
-50.0 -100.0 -150.0 -200.0 -250.0 -300.0 -350.0 -400.0 -450.0 -500.0 -550.0 -600.0 -650.0 -700.0 -750.0 -800.0 -850.0 -900.0 -950.0
0.6
x2 0.4
0.2
1
0
0.2
0.4
x1
0.6
0.8
t=0.0000
0.6
x2 0.4
0.2
1
0.2
0.4
x1
0.6
0.8
-0.005 -0.010 -0.015 -0.020 -0.025 -0.030 -0.035 -0.040 -0.045 -0.050 -0.055 -0.060 -0.065 -0.070 -0.075 -0.080 -0.085 -0.090 -0.095
0.6
x2 0.4
0.2
0
0.2
0.4
x1
0.6
0.8
1
u2 -50.0 -100.0 -150.0 -200.0 -250.0 -300.0 -350.0 -400.0 -450.0 -500.0 -550.0 -600.0 -650.0 -700.0 -750.0 -800.0 -850.0 -900.0 -950.0
0.6
x2 0.4
0.2
0
1
0
0.2
0.4
x1
0.6
0.8
1
t=0.0500
u2 -0.5 -1.0 -1.5 -2.0 -2.5 -3.0 -3.5 -4.0 -4.5 -5.0 -5.5 -6.0 -6.5 -7.0 -7.5 -8.0 -8.5 -9.0 -9.5
0.8
0.6
x2 0.4
0.2
0
u2
t=0.0000
t=0.0005
0.8
1
0.8
0
1
u2 -0.5 -1.0 -1.5 -2.0 -2.5 -3.0 -3.5 -4.0 -4.5 -5.0 -5.5 -6.0 -6.5 -7.0 -7.5 -8.0 -8.5 -9.0 -9.5
0
The initial condition is particularly set at t ¼ 0:01, the ending time and space domain are defined differently for different cases.
1
0.8
0
The simulation of diffusive N-wave is another classical Burgers problem, which can be used to test the performance of proposed numerical model which describes strong diffusion phenomena. In this problem the Reynolds number affects the rate of diffusion and the effective range of the diffusive N-wave. The analytical solution is given by [38]: pffiffiffiffiffiffiffiffi 1=t expð Rex21 =4tÞ x1 pffiffiffiffiffiffiffiffi u1 ðx1 ; tÞ ¼ ð40Þ ; x1 A Ω; t 4 0: t 1 þ 1=t expð Rex21 =4tÞ
u2
0.8
0
6.2. One-dimensional problem with diffusive N-wave
1
0
0.2
0.4
x1
0.6
0.8
1
u2
t=5.0000
-0.005 -0.010 -0.015 -0.020 -0.025 -0.030 -0.035 -0.040 -0.045 -0.050 -0.055 -0.060 -0.065 -0.070 -0.075 -0.080 -0.085 -0.090 -0.095
0.8
0.6
x2 0.4
0.2
0
0
0.2
0.4
x1
0.6
0.8
1
Fig. 8. The evolution of u2 for Example 6.3, (a) Re ¼ 0:01, t ¼ 0, (b) Re ¼ 0:01, t ¼ 0:0005, (c) Re ¼ 1, t ¼ 0, (d) Re ¼ 1, t ¼ 0:05, (e) Re ¼ 100, t ¼ 0, and (f) Re ¼ 100, t¼ 5.
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92
Starting with Re ¼ 1, we adopt time range of t A ð0:01; 10 and space domain of x1 A ½ 10; 10. Fig. 4(a) shows clear evolution with diffusive N-wave pattern, and the waveform will flatten and spread throughout the space domain. For Re ¼ 100, the time range of t A ð0:01; 1 and space domain of x1 A ½ 1; 1 is adopted to simulate the diffusive N-wave evolution until the waveform is nearly horizontal. As shown in Fig. 4(b), the tendency of u1 for Re ¼ 100 is very similar to the case of Re ¼ 1, except the scale is now smaller. Therefore, it is no surprise to see the profiles of error
87
evolutions are also similar to those shown in Fig. 5(a–d), where both error distributions resulted with good accuracy. 6.3. Two-dimensional problem in L-shaped irregular domain From former examples the proposed numerical model is available to obtain good solution for 1D Burgers problem, now consider a more complicated problem with two-dimensional system equations in irregular domain. The diagram of L-shaped computational domain and point distribution is considered as Fig. 6, the analytical solutions for
Fig. 9. Error distribution for Example 6.3 of u1 with different number of points in space (a) EL2 of Re ¼ 0:01, Δt ¼ 10 5 , (b) ERel of Re ¼ 0:01, Δt ¼ 10 5 , (c) EL2 of Re ¼ 1, Δt ¼ 10 2 , (d) ERel of Re ¼ 1, Δt ¼ 10 2 , (e) EL2 of Re ¼ 100, Δt ¼ 10 1 , and (f) ERel of Re ¼ 100, Δt ¼10 1.
88
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92
this problem are: 8 10expð 5x1 þ 5x2 þ 50t > Re Þ > < u1 ðx1 ; x2 ; tÞ ¼ Re½1 þ expð 5x1 þ 5x2 þ 50tÞ Re
10expð 5x1 þ 5x2 þ 50t > Re Þ > : u2 ðx1 ; x2 ; tÞ ¼ Re½1 þ expð 5x1 þ 5x2 þ 50t Re Þ
1
;
ðx1 ; x2 Þ A Ω; t 4 0;
ð41Þ
where the ending time is selected independently for different Reynolds numbers in this problem. Fig. 7(a–f) and Fig. 8(a–f) demonstrate the evolution of velocity u1 and u2 with different Reynolds numbers. It can be observed that different Reynolds numbers lead to similar pattern of evolution, only different rates of evolution and magnitude of velocity fields.
t=0.00
1
ERel of u1 1.0E-03 7.0E-04 5.0E-04 2.0E-04 1.0E-04 5.0E-05
0.8
0.6
t=0.00 ERel of u2 1.0E-03 7.0E-04 5.0E-04 2.0E-04 1.0E-04 5.0E-05
0.8
0.6
x2
x2 0.4
0.4
0.2
0.2
0
1
0
0.2
0.4
x1
0.6
0.8
0
1
t=0.05
1
ERel of u1 1.0E-03 7.0E-04 5.0E-04 2.0E-04 1.0E-04 5.0E-05
0.8
0.6
0
0.2
0.4
x1
0.6
1
t=0.05 ERel of u2 1.0E-03 7.0E-04 5.0E-04 2.0E-04 1.0E-04 5.0E-05
0.8
0.6
x2
0.8
x2 0.4
0.4
0.2
0.2
0
1
0
0.2
0.4
x1
0.6
0.8
0
1
t=0.10
1
ERel of u1 1.0E-03 7.0E-04 5.0E-04 2.0E-04 1.0E-04 5.0E-05
0.8
0.6
0
0.2
0.4
x1
0.6
1
t=0.10 ERel of u2 1.0E-03 7.0E-04 5.0E-04 2.0E-04 1.0E-04 5.0E-05
0.8
0.6
x2
0.8
x2 0.4
0.4
0.2
0.2
0
0
0.2
0.4
x1
0.6
0.8
1
0
0
0.2
0.4
x1
0.6
0.8
1
Fig. 10. Space distribution of ERel for Example 6.3 with Re ¼ 1, Δt ¼ 10 3 , 1027 global points and 50 local points (a) ERel of u1 at t ¼ 0, (b) ERel of u2 at t ¼ 0, (c) ERel of u1 at t ¼ 0:05, (d) ERel of u2 at t ¼ 0:05, (e) ERel of u1 at t ¼ 0:1, and (f) ERel of u2 at t ¼0.1.
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92
89
Fig. 11. Error evolution for Example 6.3 of u1 with different time interval (a) EL2 of Re ¼ 1, and (b) ERel of Re¼ 1.
Detailed analysis of the proposed numerical model, including different time intervals, the number of global points, and the number of local points are examined in this example. Fig. 9(a–f) shows the error evolutions for the velocity u1 . From observation, adding global points or expanding local influence domain can achieve better accuracy. By comparing the numerical solutions of adding global points and expanding local influence domain, it is obvious that expanding local influence domain can obtain the same accuracy as adding global points. Next, the relative error distributions in space shown as Fig. 10(a–f), demonstrate both u1 and u2 have the same order of accuracy during different time. Furthermore, Fig. 10(a–f) exhibits that the relative error is higher while the location is closer to the boundary. This behavior is caused by the shape of the local influence area, which is different for each point near the boundaries. This will affect the accuracy due to certain amount of nearest-points which was chosen to depict the local influence area. Moreover the selection of local influence area may also be different for various variables, Fig. 10(c–f) shows the error distributions which are more vertically for u1 and horizontally for u2 , despite the same magnitude and similar space distributions are used for the solutions. Fig. 10 indicates further improvements may be implemented to the proposed model by applying a better method for selecting the local influence area. Finally, Fig. 11(a and b) displays the error distribution of Re ¼ 1 with 1027 global points, showing evidence that the accuracy can be also improved by using smaller time interval. 6.4. Three-dimensional problem in irregular domain For the final example, the three-dimensional Burgers problem in irregular domain is considered. The computational domain of this problem is given as Fig. 12: 8 9 x1 ¼ P sin ψ cos θ; > > < = ∂Ω ¼ ðx1 ; x2 ; x3 Þ x2 ¼ P sin ψ sin θ; 0 o ψ r π ; 0 o θ r 2π ; > > : ; x ¼ P cos ψ ; 3
ð42Þ where ellipsoid radius P is defined as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 P ¼ 3ð cos 3θ þ 8 sin θÞ1=3 :
ð43Þ
The analytical solution for this problem is given by 8 u ðx ; x ; x ; tÞ ¼ > > > 1 1 2 3 < u2 ðx1 ; x2 ; x3 ; tÞ ¼ > > > : u3 ðx1 ; x2 ; x3 ; tÞ ¼
2 expð t=ReÞ cos x1 sin x2 sin x3 Re 1 þ expð t=ReÞ sin x1 sin x2 sin x3 2 expð t=ReÞ sin x1 cos x2 sin x3 Re 1 þ expð t=ReÞ sin x1 sin x2 sin x3 2 expð t=ReÞ sin x1 sin x2 cos x3 Re 1 þ expð t=ReÞ sin x1 sin x2 sin x3
Fig. 12. The irregular domain for Example 6.4, a total number of 2547 points are shown in this figure. The gray shade is the boundary ∂Ω, the black spheres are the boundary points and the red spheres depict the global points. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
In Fig. 13(a–f), the plots show the velocity distribution of u1 ; u2 and u3 on the planes of x1 ¼ 1:5; x2 ¼ 1:5 and x3 ¼ 1:5. Fig. 13(a–c) shows the distinctive N-wave patterns for all three velocity components u1 ; u2 and u3 at t ¼ 0. While time elapses, the diffusive N-wave propagates through the domain, and the diffused solution at t ¼ 1 are shown in Fig. 13(d–f). The error evolutions in Fig. 14(a–d) indicate that for either Re ¼ 1 or Re ¼ 100, the numerical solutions obtained by the proposed model consistently showing acceptable accuracy. The error evolution rises during the late period of experiment due to the fact that the scale of the solution becomes very small. Therefore, even the smallest difference will be magnified in relative error, and the proposed numerical scheme is still able to achieve good accuracy under such circumstances. By further testing the problem with different numbers of global and local points, the coherent behavior of increasing accuracy while adding global points or expanding local influence domain are also demonstrated in error evolutions. 7. Conclusions
; ðx1 ; x2 ; x3 Þ A Ω; t 40:
ð44Þ
The proposed numerical model of FD-LMAPS using the Cole– Hopf transformation has proven the fine capacity for simulation of
90
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92
Fig. 13. Evolution of the velocity for Example 6.4 distributed on planes of x1 ¼ 1:5, x2 ¼ 1:5 and x3 ¼ 1:5. Conditions given as Re ¼ 1, 10,825 global points and 50 local points, (a) u1 at t ¼ 0, (b) u2 at t ¼ 0, (c) u3 at t ¼ 0, (d) u1 at t ¼ 1, (e) u2 at t ¼ 1, and (f) u3 at t ¼ 1.
multi-dimensional Burgers problems in regular and irregular domains with sufficient accuracy. Although with limitation to only solve problems with Dirichlet type boundary conditions, the characteristics of wave propagation and diffusive N-wave have been simulated with result achieving maximum relative error within the order of 10 3 . Four numerical experiments, namely, the one-dimensional problem with wave propagation, the onedimensional problem with diffusive N-wave, the two-dimensional problem in L-shaped irregular domain and the three-dimensional
problem in irregular domain, have been analyzed in this paper. The comparisons between numerical solutions and analytical solutions have given series of results with high accuracy. The error distributions also have shown consistent behavior of achieving better accuracy while adding global points and expanding local influence domain, which can also obtain better efficiency. Hereby we conclude that the proposed numerical model of FD-LMAPS with the Cole–Hopf transformation is a promising tool for future researches and applications.
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92
91
Fig. 14. Error evolution for Example 6.4 of velocity u1 (a) EL2 of Re ¼ 1, Δt ¼ 10 2 , (b) ERel of Re ¼ 1, Δt ¼ 10 2 , (c) EL2 of Re ¼ 100, Δt ¼ 10 2 , and (d) ERel of Re¼ 100, Δt ¼10 2.
Symbol
Appendix A
All variables in this paper are formatted in italic or regular style, only the vectors and matrices are formatted in arrow or bold style. Superscripts and subscripts
Detail derivations of Cole–Hopf transformation for multidimensional Burgers equations. The Burgers equation
b d i ib j k m n s
indicates the boundary. number of dimensions. ith point of global domain. ibth point of boundary. jth point of global domain. kth point of local domain. mth point of local domain. indicates nth time step. indicates the initial time step.
Acknowledgment Financial supports of this paper are provided by the National Science Counsel of Taiwan under the grant No. NSC100-3332-E002–020. It is greatly appreciated. The fourth author acknowledges the support of Distinguished Overseas Visiting Scholar Fellowship provided by the Ministry of Education of China.
, , , , ∂U 1 þ U U∇ U ¼ ∇2 U : Re ∂t
The Cole–Hopf transformation ,
Uqþ
2 ∇q ¼ 0andq a 0: Re
Cole-Hopf transformation for various differential terms 8 2 ∂q > udi ¼ Req > ∂xdi > > > 2 > ∂udi ∂q ∂q 2 ∂ q 2 > ¼ > > ∂t Req ∂xdi ∂t þ Req2 ∂xdi ∂t > > > 2 > 2 > ∂udi ∂q 2 ∂ q 2 > ¼ Req þ Req > 2 ∂x > ∂x2di di < ∂xdi : ∂udi ∂2 q ∂q ∂q 2 2 ¼ þ > Req ∂xdi ∂xdj Req2 ∂xdi ∂xdj > > ∂xdj > 3 > > 2 3 ∂q ∂2 q ∂q > 2 ∂ q 6 4 > ∂ u2di ¼ Req þ Req > 2 ∂x 2 Req3 ∂x ∂x3di ∂xdi > di ∂xdi di > > > 2 > > ∂2 udi ∂3 q ∂2 q ∂q ∂q ∂2 q ∂q ∂q 2 4 2 4 > > : ∂x2 ¼ Req ∂x x2 þ Req2 ∂xdi ∂xdj ∂xdj þ Req2 ∂xdi ∂x2 Req3 ∂xdi ∂xdj dj
di dj
dj
subscript di and dj depicts different dimension number d:
92
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92
Substituting the differential terms after Cole–Hopf transformation into the Burgers equation and obtain ∇q ∂q 1 2 ∂q 1 2 ∇ q ∇ ∇ q ¼ 0: q ∂t Re ∂t Re Using 1-D Burgers equation as an example ∂u1 ∂u1 ∂t þ u1 ∂x1 -
2
1 ∂ u1 ¼ Re ∂x2 1
" 2 # 2 ∂ q 2 ∂q ∂q 2 ∂q 2 ∂2 q 2 ∂q þ þ þ Req ∂x1 ∂t Req2 ∂x1 ∂t Req ∂x1 Req ∂x21 Req2 ∂x1 2
" 3 # 1 2 ∂3 q 6 ∂q ∂2 q 4 ∂q þ Re Req ∂x31 Req2 ∂x1 ∂x21 Req3 ∂x1 " 2 # ∂2 q 1 ∂q ∂q ∂q 2 ∂2 q 2 ∂q þ þ þ - ∂x1 ∂t q ∂x1 ∂t ∂x1 Req ∂x21 Req2 ∂x1 3 1 ∂3 q 3 ∂q ∂2 q 2 ∂q þ ¼ 3 2 2 Re ∂x1 Req ∂x1 ∂x1 Req ∂x1 3 ∂2 q 1 ∂q ∂q 2 ∂q ∂2 q 2 ∂q þ þ - ∂x1 ∂t q ∂x1 ∂t Req ∂x1 ∂x21 Req2 ∂x1 3 1 ∂3 q 3 ∂q ∂2 q 2 ∂q ¼ þ Re ∂x31 Req ∂x1 ∂x21 Req2 ∂x1 ¼
∂2 q 1 ∂q ∂q 1 ∂3 q 1 ∂q ∂2 q ¼ þ þ 3 ∂x1 ∂t q ∂x1 ∂t Re ∂x1 Req ∂x1 ∂x21 ! ! 1 ∂q ∂q 1 ∂2 q ∂ ∂q 1 ∂2 q ¼ 0: q ∂x1 ∂t Re ∂x21 ∂x1 ∂t Re ∂x21 -
It is obvious that diffusion equation satisfies the Burgers equations transformed via Cole–Hopf transformation ∂q 1 2 ∇ q ¼ 0: ∂t Re Therefore diffusion equation is applied as the governing equation after Cole–Hopf transformation in this study. References [1] Burgers JM. A mathematical model illustrating the theory of turbulence. Adv Appl Mech 1948;1:171–99. [2] Yao G, Chen CS, Kolibal J. A localized approach for the method of approximate particular solutions. Comput Math Appl 2011;61:2376–87. [3] Bahadir AR. A fully implicit finite-difference scheme for two-dimensional Burgers' equations. Appl Math Comput 2003;137:131–7. [4] Fu S, Hodges B. Time-centered split method for implicit discretization of unsteady advection problems. J Eng Mech 2009;135:256–64. [5] Young DL, Fan CM, Hu SP, Atluri SN. The Eulerian–Lagrangian method of fundamental solutions for two-dimensional unsteady Burgers' equations. Eng Anal Bound Elem 2008;32:395–412. [6] Hopf E. The partial differential equation ut þ uux ¼ μxx. Commun Pure Appl Math 1950;3:201–30. [7] Cole JD. On a quasilinear parabolic equation occurring in aerodynamics. Q Appl Math 1951;9:225–36. [8] HSL. A collection of FORTRAN codes for large scale scientific computation, 〈http://www.hsl.rl.ac.uk〉. [9] Kutluay S, Bahadir AR, Özdeş A. Numerical solution of one-dimensional Burgers equation: explicit and exact-explicit finite difference methods. J Comput Appl Math 1999;103:251–61.
[10] Kadalbajoo MK, Awasthi A. A numerical method based on Crank–Nicolson scheme for Burgers' equation. Appl Math Comput 2006;182:1430–42. [11] Öziş T, Aksan EN, Özdeş A. A finite element approach for solution of Burgers' equation. Appl Math Comput 2003;139:417–28. [12] Varoǧ lu E, Finn WD. Space-time finite element incorporating characteristics for the Burgers' equation. Int J Numer Methods Eng 1980;16:171–84. [13] Chino E, Tosaka N. Dual reciprocity boundary element analysis of timeindependent Burgers' equation. Eng Anal Bound Elem 1998;21:261–70. [14] Kakuda K, Tosaka N. The generalized boundary element approach to Burgers' equation. Int J Numer Methods Eng 1990;29:245–61. [15] Chen JK, Beraun JE. A generalized smoothed particle hydrodynamics method for nonlinear dynamic problems. Comput Method Appl M 2000;190:225–39. [16] Golberg MA, Chen CS. Improved multiquadric approximation for partial differential equations. Eng Anal Bound Elem 1996;18:9–17. [17] Hon YC, Mao XZ. An efficient numerical scheme for Burgers' equation. Appl Math Comput 1998;95:37–50. [18] Cheng AHD, Golberg MA, Kansas EJ, Zammito G. Exponential convergence and H-c multiquadric collocation method for partial differential equations. Numer Methods Partial Differ Equ 2003;19:571–94. [19] Li JC, Hon YC, Chen CS. Numerical comparisons of two meshless methods using radial basis functions. Eng Anal Bound Elem 2002;26:205–25. [20] Chen R, Wu Z. Applying multiquadric quasi-interpolation to solve Burgers' equation. Appl Math Comput 2006;172:472–84. [21] Young DL, Hu SP, Fan CM, Chen CW. Method of fundamental solutions and Cole–Hopf transformation for Burgers' equations. Comput Fluid Dyn J 2006;14:454–67. [22] Lin CY, Gu MH, Young DL. The time-marching method of fundamental solutions for multi-dimensional telegraph equations. CMC-Comput Mater Con 2010;18:43–68. [23] Tsai CC, Chen CS, Hsu TW. The method of particular solutions for solving axisymmetric polyharmonic and poly-Helmholtz equations. Eng Anal Bound Elem 2009;33:1396–402. [24] Wen PH, Chen CS. The method of particular solutions for solving scalar wave equations. Int J Numer Method Biomed Eng 2010;26:1878–89. [25] Atluri SN, Zhu T. A new meshless local Petrov–Galerkin (MLPG) approach in computational mechanics. Comput Mech 1998;22:117–27. [26] Lin H, Atluri SN. The meshless local Petrov–Galerkin (MLPG) method for solving incompressible Navier–Stokes equations. CMES-Comput Model Eng 2001;2:117–42. [27] Lee CK, Liu X, Fan SC. Local multiquadric approximation for solving boundary value problems. Comput Mech 2003;30:396–409. [28] Shu C, Ding H, Yeo KS. Computation of incompressible Navier–Stokes equations by local RBF-based differential quadrature method. CMES-Comput Model Eng 2005;7:195–205. [29] Zhou CH, Shu C. A local domain-free discretization method for simulation of incompressible flows over moving bodies. Int J Numer Meth Fl 2011;66:162–82. [30] Kosec G, Saler B. H-adaptive local radial basis function collocation meshless method. CMC-Comput Mater Con 2011;26:227–53. [31] Shu C, Ding H, Yeo KS. Local radial basis function-based differential quadrature method and its application to solve two dimensional incompressible Navier– Stokes equations. Comput Methods Appl Mech Eng 2003;192:941–54. [32] Yao G, Šarler B, Chen CS. A comparison of three explicit local meshless methods using radial basis functions. Eng Anal Bound Elem 2011;35:600–9. [33] Karageorghis A, Chen CS, Smyrlis YS. Matrix decomposition RBF algorithm for solving 3D elliptic problems. Eng Anal Bound Elem 2009;33:1368–73. [34] Gu MH, Fan CM, Young DL. The method of fundamental solutions for the multi-dimensional wave equations. J Mar Sci Tech-Taiw 2011;19:586–95. [35] Rippa S. An algorithm for selecting a good value for the parameter c in radial basis function interpolation. Adv Comput Math 1999;11:193–210. [36] Wertz J, Kansas EJ, Ling L. The role of the multiquadric shape parameters in solving elliptic partial differential equations. Comput Math Appl 2006;51:1335–48. [37] Chen CS, Fan CM, Wen PH. The method of particular solutions for solving certain partial differential equations. Numer Methods Partial Differ Equ 2012;28:506–22. [38] Polyanin AD, Zaitsev VF. Handbook of nonlinear partial differential equations. 1st ed.. Florida: Chapman and Hall/CRC; 2004.
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Localized method of approximate particular solutions with Cole–Hopf transformation for multi-dimensional Burgers equations C.Y. Lin a,b, M.H. Gu c, D.L. Young a,b,n, C.S. Chen d,e a
Department of Civil Engineering, National Taiwan University, Taipei 10617, Taiwan Hydrotech Research Institute, National Taiwan University, Taipei 10617, Taiwan c National Center for Research on Earthquake Engineering, Taipei 10668, Taiwan d Department of Mathematics, The University of Southern Mississippi, Hattiesburg, MS, USA e Department of Engineering Mechanics, Hohai University, Nanjing, Jiangsu 210098, China b
art ic l e i nf o
a b s t r a c t
Article history: Received 23 August 2013 Accepted 28 November 2013 Available online 25 December 2013
The Burgers equations depict propagating wave with quadratic nonlinearity, it can be used to describe nonlinear wave propagation and shock wave, where the nonlinear characteristics cause difficulties for numerical analysis. Although the solution approximation can be executed through iterative methods, direct methods with finite sequence of operation in time can solve the nonlinearity more efficiently. The resolution for nonlinearity of Burgers equations can be resolved by the Cole–Hopf transformation. This article applies the Cole–Hopf transformation to transform the system of Burgers equations into a partial differential equation satisfying the diffusion equation, and uses a combination of finite difference and the localized method of approximate particular solution (FD-LMAPS) for temporal and spatial discretization, respectively. The Burgers equations with behaviors of propagating wave, diffusive N-wave or within multi-dimensional irregular domain have been verified in this paper. Effectiveness of the FD-LMAPS has also been further examined in some experiments, and all the numerical solutions prove that the FD-LMAPS is a promising numerical tool for solving the multi-dimensional Burgers equations. & 2013 Elsevier Ltd. All rights reserved.
Keywords: Burgers equations Cole–Hopf transformation Diffusion equation Meshless methods Irregular domain Localized method of approximate particular solutions
1. Introduction The Burgers equations, are quasi-linear partial differential equations (PDEs) consisting an unsteady term, a nonlinear convection term and a diffusion term. The Burgers equations depict propagating wave with quadratic nonlinearity, it can be used to describe nonlinear wave propagation and shock wave [1]. The Burgers equations have been widely applied for mathematic analysis and problems of fluid mechanics. In this paper we approximate the solutions of the Burgers equations by the localized method of approximate particular solutions (LMAPS) [2] which is a novel developed numerical scheme using radial basis functions (RBFs). In order to overcome the nonlinearity of the Burgers equations, the implementation of the LMAPS to solve the resulting diffusion equation of applying Cole–Hopf transformation is considered in this paper. Basic solution approximation for nonlinear equations usually requires iterative approaches such as the Newton–Ralphson or Picard method [3,4]. However, the convergence of iterative approaches is not
n Corresponding author at: Department of Civil Engineering, National Taiwan University, Taipei 10617, Taiwan. Tel./fax: þ 886 2 2362 6114. E-mail address: [email protected] (D.L. Young).
0955-7997/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.enganabound.2013.11.019
always guaranteed. In contrast to iterative approaches, direct approaches for approximating the solution of the Burgers equations can overcome the nonlinearity through finite time sequence which avoids the issue of convergence. Such direct approaches can be easily found in the literature for solving Burgers equations, for example the well-known Eulerian–Lagrangian method [5], or the direct approach that will be applied in this paper — the Cole–Hopf transformation [6,7]. The Cole–Hopf transformation was introduced by Hopf and Cole [6,7] to deal with the nonlinearity of the viscous Burgers equations. After applying the Cole–Hopf transformation, the system of Burgers equations can be transformed into a single variable PDE satisfying the diffusion equation. The analytical solution of the diffusion equation can be simply found by solving the Cauchy problem in regular domain. Therefore, it provides a method to obtain analytical solutions for the Burgers equations. Although the Cole–Hopf transformation was originally developed for deriving analytical solutions, in this research the Cole–Hopf transformation is only applied to deal with the nonlinearity of Burgers equations for further numerical processes. For the numerical implementation, it is necessary for the essential conditions, i.e. initial and boundary conditions, to be transformed according to the definition of the Cole-Hopf transformation. This implies both the governing equation and the essential
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92
conditions are transformed into the single variable form. Although the Burgers equations have been transformed into single variable form, the number of equations that depict the essential conditions remains the same, which indicates the linear system generated for dealing with multi-dimensional problems will be overdetermined. Hereby a numerical linear solver named MA49 using least-squares QR method [8] with capability for dealing with the sparse over-determined system is applied in this article. After overcoming the nonlinearity by the Cole–Hopf transformation, further solution approximation is required. For numerical analysis, there are well-known and widely applied mesh-dependent numerical methods, such as the finite difference method (FDM) [3,9,10], the finite element method (FEM) [11,12] or the boundary element method (BEM) [13,14]. All these classical numerical methods are popular and effective but require difficult procedures of mesh generation or numerical quadrature. In order to avoid such difficulties, various meshless numerical methods such as the smoothed particle hydrodynamics (SPH) [15], the multiquadrics collocation method (MQ) [16–20], the method of fundamental solutions (MFS) [5,21,22], the method of particular solutions (MPS) [23,24] and the meshless local Petrov–Galerkin method (MLPG) [25,26], to name just a few, have been developed. For global domain-type meshless methods with RBFs, the resultant interpolation matrix of MQ or the MPS is often dense and ill-conditioned which makes the computation inefficient and unstable. These disadvantages cause the global type numerical methods to be restrained for large-scale computational problems. Hereby extension researches on localized formulation of global domain-type meshless numerical methods have been intensively studied among recent literature [2,27–31]. In this paper we focus on the localized formulation to modify the method of approximate particular solutions (MAPS), also called the localized method of approximate particular solutions (LMAPS). The main concept of localization was introduced by Yao, Chen, and Kolibal [2]. The LMAPS is based on the global MAPS with modification of choosing limited number of weighting points within the entire local domain. The breakthrough of the localization technique allows dense matrix of the original global linear system to reduce into a sparse linear system. This will decrease the risks of ill-conditioned linear system and lessen computational memory loading. Therefore the LMAPS as a meshless method has the advantages without requiring any mesh generation or numerical quadrature, and has also avoided the drawbacks of the global MAPS or MPS. The processes of localization increase the efficiency for solving linear systems, meanwhile maintaining appreciable accuracy of solutions. There is a very similar localized meshless numerical method called the local radial basis function differential quadrature method (LRBFDQ) [28,31]. In numerical algorithm, these two localized meshless numerical methods share many common features; however the basis function chosen has made them mathematically different. For the basis function approximating the solutions, LRBFDQ applies the differential quadrature (DQ), while the LMAPS uses the integrated radial basis functions (integrated RBF). Mathematically, integrating a function will obtain a solution function with an unknown constant, causing differences in solutions obtained from these two methods. Despite the differences, both schemes show very good performance in numerical simulation due to the same localization technique they shared. The technique itself will be further described in this paper and Ref. [2], and the comparison of related methods can be found in Ref. [32]. Finally the time domain and space domain are considered with different approaches of discretization. The time domain is discretized by the Crank–Nicolson scheme of FDM, where the Crank– Nicolson scheme is an implicit unconditional stable scheme with second-order of accuracy. For the space domain, the LMAPS is
79
chosen to discretize the space domain and approximate the solutions. We summarize here the techniques applied for this study: the Cole–Hopf transformation transforms the Burgers equations into a PDE satisfying diffusion equation, the Crank– Nicolson scheme approximates transient term and the LMAPS applied for space discretization. Combining these methods and integrating into the proposed numerical model of FD-LMAPS with the Cole–Hopf transformation, we will carry out numerical experiments to solve the multi-dimensional Burgers equations. There are four numerical experiments in this paper which are named a one-dimensional Burgers problem with propagating wave, a one-dimensional Burgers problem with diffusive N-wave, a two-dimensional Burgers problem in L-shaped domain, and a three-dimensional Burgers problem in an irregular domain. Further details of the proposed numerical model are provided in the following sections.
2. Governing equations and essential conditions 2.1. Governing equations The Burgers equations are quasi-linear PDEs which include nonlinear convection term and diffusion term. For the global domain Ω including the boundary ∂Ω, the Burgers equations in vector form listed as ,
, , , ∂U 1 þ ðU U∇ÞU ¼ ∇2 U ; Re ∂t , ,
ð1Þ ,
where U ðx ; tÞ ¼ ðu1 ; ⋯; ud Þ is the vector of velocity, x ¼ ðx1 ; ⋯; xd Þ A Rd , d ¼ 1; 2; 3, and Re is the Reynolds number. 2.2. Initial conditions and boundary conditions The initial conditions and boundary conditions are the essential conditions for obtaining unique solutions for the governing equations. Only Dirichlet-type boundary conditions are considered in this study, the reason will be explained in Section 3.2. The initial conditions in general form is depicted as , , ,, U ðx ; tÞ ¼ s ðx Þ; ð2Þ t¼0
and the boundary conditions is ,, , , ¼ β ðx ; tÞ: U ðx ; tÞ
ð3Þ
,
x A ∂Ω ,,
,,
Here s ðx Þ ¼ ðs1 ; ⋯; sd Þ is the given initial ,velocity, β ðx ; tÞ ¼ ðβ 1 ; ⋯; β d Þ is the given boundary conditions of U . 3. Theoretical model The Cole–Hopf transformation for the governing equations and essential conditions are introduced with some detailed derivations. The demonstrations, the applications, and restrictions of the Cole–Hopf transformation applied in this research will be followed. 3.1. Cole–Hopf transformation for the Burgers equations In order to transform the system of Burgers equations into diffusion equations, the Cole–Hopf transformation is introduced as follows [6,7]: ,
Uqþ
2 ∇q ¼ 0 and q a 0; Re ,
ð4Þ
where qðx ; tÞ is the undetermined variable after the transformation.
80
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92
Substituting Eq. (4) into the system of Burgers equations in Eq. (1), we have ∇q ∂q 1 2 ∂q 1 2 ∇ q ∇ ∇ q ¼ 0; ð5Þ q ∂t Re ∂t Re
straightforwardly, integrated basis functions corresponding to the Laplace operator are chosen to describe the particular solutions of Eq. (9).
the solution of following diffusion equation automatically satisfies Eq. (5):
4.2. Localized method of particular solutions
∂q 1 2 ∇ q ¼ 0: ∂t Re
In this paper, the LMAPS is applied to discretize the space , domain of the diffusion equation. Let x i A Ω; i ¼ 1; ⋯; N be the global points, where N is the total number of global points. From the MPS, q is approximated by a series of integrated basis functions of F,
ð6Þ
If the problem is well-posed, then the unique solution of Eq. (6) exists. Detailed derivation for the Cole–Hopf transformation in multidimensional case is given in the Appendix A.
N
,
,
j¼1
The essential conditions include both initial conditions and boundary conditions in this article. By the definition of the Cole–Hopf transformation as Eq. (4), the essential conditions in Eqs. (2) and (3) can be derived as follows: 2 , , ,, Initial condition : s ðx Þqðx ; tÞ þ ∇qðx ; tÞ ¼ 0; ð7Þ Re t¼0 , 2 , , , Boundary condition : β ðx ; tÞqðx ; tÞ þ ∇qðx ; tÞ ¼ 0: , Re x A ∂Ω
,
qðx ; tÞ ffi ∑ αj Fð‖x x j ‖Þ;
3.2. Cole–Hopf transformation for essential conditions
,
While providing Dirichlet type conditions of U , the essential conditions of Eqs. (7) and (8) become Robin type conditions of q. Therefore the proposed numerical scheme can discretize the transformed essential conditions, of q directly. However if the Neumann type conditions of U were provided for essential conditions, the transformed essential conditions will become nonlinear. The nonlinearity should be resolved for further analysis or other approximations. In this paper, only the Dirichlet type conditions will be considered. Therefore the governing equation and essential conditions to be solved are described by Eqs. (6)–(8).
4. Numerical methods
x A Ω:
ð10Þ
where ∇2 F ¼ f , f is the basis function and αj are the weighting coefficients to be determined. For the LMAPS, we can create a local influence domain ωi for , , each global point x i . The local points are defined as x i;k A ωi ; k ¼ 1; ⋯; NL (see Fig. 1). By the MPS, on each local influence domain ωi , we have NL
,
ð8Þ
,
,
,
qðx ; tÞ ffi ∑ αi;m Fð‖x x i;m ‖Þ; m¼1
,
x A ωi ;
ð11Þ
where αi;m is the weighting coefficients of the mth local point in the local influence domain ωi . By collocation method on Eq. (11), qi ¼ Fi αi ;
ð12Þ
T , , , , ; Fi ¼ Fð‖x i;k x i;m ‖Þ ; where qi ¼ q x i;1 ; t ; ⋯; q x i;NL ; t NLNL αi ¼ αi;1 ; ⋯; αi;NL T . If Fi is invertible, Eq. (12) can be rewritten as
αi ¼ Fi 1 qi :
ð13Þ
Now the Laplace operator ∇2 can be applied to Eq. (11) to obtain the following expression, ,
,
NL
,
,
∇2 qðx ; tÞ ffi ∇2 qðx i ; tÞ ¼ ∑ αi;m f ð‖x i x i;m ‖Þ ¼ f i αi :
ð14Þ
m¼1
This section introduces the implementation of numerical methods, the Crank–Nicolson scheme is used for discretization of the time domain and the LMAPS is applied for discretization of the space domain.
, , , , where f i ¼ f ð‖x i x i;1 ‖Þ; ⋯; f ð‖x i x i;NL ‖Þ . From Eqs. (13) and (14), we have
4.1. Crank–Nicolson scheme
where φi ¼ f i Fi 1 . In order to transfer the above localsystem into global system, it requires to map the local vector φi ¼ φi;1 ; ⋯; φi;NL into the global
The time-marching Crank–Nicolson scheme [10] is applied to discretize the time domain for numerical experiments in this study. The second-order accuracy and unconditional stable Crank–Nicolson scheme is selected for the time domain of diffusion equation, which can be written as follows: qn þ 1
Δt 2 n þ 1 Δt 2 n ∇ q ∇ q þ OðΔt 2 Þ; ¼ qn þ 2Re 2Re
,
∇2 qðx i ; tÞ ¼ f i Fi 1 qi ¼ φi qi ;
ð15Þ
ð9Þ
where Δt is the time interval, OðΔt 2 Þ is the truncation error, superscript of q denotes the number of time step and n marks the time step t n ¼ nΔt. Eq. (9) can be considered as a modified Helmholtz equation with time-dependent source term, the particular solutions of Eq. (9) can be described by integrated basis functions corresponding to either Laplace operator or modified Helmholtz operator. For particular solutions of modified Helmholtz equation, the modified Bessel functions make it difficult to obtain suitable RBF for applications. Since the integrated RBFs corresponding to the Laplace operator do not share the same problem and can be implemented to various PDEs other than Helmholtz equation more
,
Fig. 1. The diagram of global domain Ω, local influence domain ωi of point x i , global , , points x and local point x i;k .
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92
vector Φi where
Φi ¼ Φi;1 ; ⋯; Φi;N ; Φi;j ¼
8 > < 0;
, xj2 =
> : φi;k ;
, xj
ωi ,
¼ x i;k A ωi
:
ð16Þ
Therefore, the system extends into the global domain Ω, ∇2 q ¼ Φq;
where Φ ¼ Φ1 ; ⋯; ΦN
T
T , , , and q ¼ qðx 1 ; tÞ; ⋯; qðx N ; tÞ .
ð17Þ
Applying the similar procedures to the gradient operator ∇, by defining g ¼ ∇F, we have ,
,
NL
,
,
∇qðx ; tÞ ffi ∇qðx i ; tÞ ¼ ∑ αi;k gð‖x i x i;m ‖Þ ¼ gi αi ; m¼1
ð18Þ
,
ð19Þ ∇qðx i ; tÞ ¼ gi Fi 1 qi ¼ γi qi ; , , , , and γi ¼ gi Fi 1 . The where gi ¼ g ‖x i x i;1 ‖ ; ⋯; g ‖x i x i;NL ‖ local system can be transferred into global system by mapping the local vector γi ¼ γ i;1 ; ⋯; γ i;NL into the global vector 8 , > < 0; xj 2 = ωi ; Γi ¼ Γ i;1 ; ⋯; Γ i;N ; Γ i;j ¼ ð20Þ , , > : γ i;k ; x j ¼ x i;k A ωi ; ∇q ffi Γq:
ð21Þ
81
T Note that Γ ¼ Γ1 ; ⋯; ΓN and sub-elements Γi ; i ¼ 1; ⋯; N are defined as Eq. (20). Further details of the global extension of φi and γi can be found in Ref. [2]. Each linear term of the governing equations or the essential conditions can be discretized by the LMPS through similar processes. Further detailed numerical procedures for solving Burgers equations will be given in Section 5.
4.3. Integrated multiquadrics radial basis functions The integrated multiquadrics radial basis functions (integrated MQ-RBFs) are used for the weighting function in this work. The MQ-RBF and the integrated MQ-RBFs are defined as follows [32–34]. The MQ-RBF: f ðr Þ ¼
pffiffiffiffiffiffiffiffiffiffiffiffiffiffi r 2 þ c2 :
ð22Þ
The 1-dimensional integrated MQ-RBF: FðrÞ ¼
pffiffiffiffiffiffiffiffiffiffiffiffiffiffi i pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1h
2c2 þ r 2 r 2 þ c2 þ r ; r 2 þ c2 þ 3c2 r ln 6
The 2-dimensional integrated MQ-RBF: pffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 1 r 2 þ c2 þ c ; FðrÞ ¼ ðr 2 þ 4c2 Þ r 2 þ c2 c3 ln 9 3
ð23Þ
ð24Þ
Fig. 2. The evolution of u1 for Example 6.1 at different time step with (a) Re ¼ 1, Δt ¼ 10 2 , 21 global points and 3 local points, (b) Re ¼ 10, Δt ¼ 10 2 , 51 global points and 3 local points, (c) Re ¼ 100, Δt ¼ 10 2 , 501 global points and 3 local points.
82
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92
The 3-dimensional integrated MQ-RBF: pffiffiffiffiffiffiffiffiffiffi 8 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4 2 2 > 1 c < 24 ð2r 2 þ5c2 Þ r 2 þ c2 þ 8r ln r þ cr þ c ; FðrÞ ¼ > : c3 ; 3
5. Numerical procedures r 4 0;
ð25Þ
r ¼ 0;
where r is considered as the distance between the collocation points and constant c is the shape parameter of MQ and taken in the order of 1.0 in this study. There are various approaches proposed for obtaining optimal shape parameters [35–37], each claiming to be optimal and consistent, suitable for applications of MQ-RBF.
In this section the detailed procedures for solving the Burgers equations are provided by using the FD-LMAPS and the Cole–Hopf transformation. 5.1. Initial conditions ,
From Section 3.2, the initial conditions of velocity field U have been derived to Eq. (7) as the Robin type condition. The LMAPS and integrated MQ-RBF were applied in each term for discretiza, tion. For each global point x i , the local influence domain ωi can be
Fig. 3. The error evolution of Example 6.1 of u1 with different number of global points, (a) EL2 of Re ¼ 1, (b) ERel of Re ¼ 1, (c) EL2 ofRe ¼ 10, (d) ERel ofRe ¼ 10, (e) EL2 of Re ¼ 100, and (f) ERel of Re¼ 100.
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92
83
ð26Þ
dN equations for describing the initial conditions, where dN ¼ d U N. Hereby the global linear system is described as follows 2 Is þ Γs q ¼ 0; ð27Þ Re dNN
where qi ¼ qðx i ; tÞ. The above local system can be extended to global system. For d-dimensional problems, there are d initial conditions for each , global point of velocity field U . Therefore, N global points implies
where Is and Γs are defined as follows. For d-dimensional cases: 2 3 2 3 I1 Γ 6⋮7 6 7 ; Γs ¼ 4 ⋮ 5 : ð28Þ Is ¼ 4 5 Id Γ dNN
established as mentioned in Section 4.2 to obtain the following local linear system: ,,
s ðx Þq þ
2 2 , ∇q ¼ 0 ) sd ðx i Þqi þ γi qi ¼ 0; Re Re ,
dNN
Fig. 4. The evolution of u1 for Example 6.2 at different time step with (a) Re ¼ 1, Δt ¼ 10 3 , 501 global points and 3 local points, t A 10 2 ; 1 , and (b) Re ¼ 100, Δt ¼ 10 4 , 2 501 global points and 3 local points, tA (10 ,1).
Fig. 5. The error evolution of Example 6.2 of u1 with different number of global points, (a) EL2 of Re ¼ 1, t A 10 2 ; 10 , (b) ERel of Re ¼ 1, t A 10 2 ; 10 , (c) EL2 of Re ¼ 100,
2 2 t A 10 ; 1 , and (d) ERel of Re ¼ 100, t A (10 ,1).
84
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92 , xi.
Similar to the numerical process of the initial conditions, the , local influence domain is established for each boundary point x ib to obtain the following local linear system from Eq. (8),
,,
β ðx ; tÞq þ
2 2 , ∇q ¼ 0 ) β d ðx ib ; tÞqi þ γi qi ¼ 0: Re Re
ð32Þ
The above local system can also be extended to global system [2]. For d-dimensional problems, each boundary point provides d , boundary conditions of velocity field U , implies NB boundary points will result in d U NB equations. The global linear system of boundary conditions can be described as follows: 2 Ib þ Γb q ¼ 0; ð33Þ Re dNBN
Fig. 6. The L-shaped domain Ω for Example 6.3, a total number of 1027 points are shown in this figure. The black border is the boundary ∂Ω, the black squares are the boundary points and the red circles depict the global points. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Note that 2 3 Id;1 6 ⋮ 7 Id ¼ 4 5; Id;N
2
, 6 7 Id;i ¼ 40; ⋯; sd ðx i ; tÞ ; ⋯; 05 |fflfflfflffl{zfflfflfflffl}
:
ð39Þ
1N
After applying the Cole–Hopf transformation, there is a single variable q for each global point. It implies N unknowns of q needs , to be approximated for satisfying the velocity field U . The homogeneous solution of over-determined sparse linear system can be solved by the singular value decomposition of MA49 [8]. 5.2. Governing equations
It is necessary to introduce the boundary conditions into the system combined with governing equation to satisfy the velocity field. Therefore, the two linear systems are combined as follows # " " # Δt Φ I Δt Φ I þ 2Re nþ1 2Re q ¼ qn : ð36Þ 2 Γb Ib þ Re ½0 ðN þ dNBÞN
ðN þ dNBÞN
By using least-squares QR method [8], the numerical solver named MA49 is applied to obtain the solution of the sparse overdetermined linear system. Algorithm.
For space discretization of Eq. (9), the LMAPS scheme is , implemented. Let x i be the interior points and by collocation method, the following local system can be obtained.
Δt 2 n þ 1 Δt 2 n ∇ q ∇ q ¼ qn þ 2Re 2Re Δ t Δt ) qni þ 1 φ qn þ 1 ¼ qni þ φ qn : 2Re i i 2Re i i
where for each boundary point with a corresponding global point, and 2 3 2 3 Jd;1 , 6 7 6 7 Jd ¼ 4 ⋮ 5; Jd;ib ¼ 40; ⋯; β d ðx ib ; tÞ ; ⋯; 05 : ð35Þ |fflfflfflfflffl{zfflfflfflfflffl} Jd;NB i 1N
3
i
where Ib and Γb are defined as follows. For d-dimensional cases: 2 3 2 3T J1 6 7 6 7 ð34Þ Ib ¼ 4 ⋮ 5; Γb ¼ 4Γ1 ; ⋯; ΓNB ; ⋯; Γ1 ; ⋯; ΓNB 5 ; |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Jd Cycledtimes
qn þ 1
ð30Þ
The above local system can be assembled into a global system. It follows that Δt Δt Φ qn þ 1 ¼ Iþ Φ qn ; ð31Þ I 2Re 2Re NN NN where matrix I is the N N identity matrix. Notice that the terms on the right hand side are all known from the previous time step. Since the governing equation after the Cole–Hopf transformation consist only one equation for each point in the space domain, this sparse linear system is considered to be fully-determined. However it is still necessary to introduce the boundary conditions to satisfy the velocity field. 5.3. Boundary conditions ,
The Dirichlet type boundary conditions of velocity field U in Eq. (3) can be derived to the Robin type boundary condition as shown in Eq. (8). Each term is derived with the technique of the LMAPS and use integrated MQ-RBF for basis function. Let NB denotes the number of boundary points. For each boundary point , x ib A ∂Ω; ib ¼ 1; ⋯; NB, there exists a corresponding global point
Step 1: Construct and solve the linear system of the initial conditions as Eq. (27). Step 2 Construct the linear system of the governing equations and boundary conditions as in Eq. (36) and obtain solution of qn þ 1 . ^ nþ1 , according to the Step 3: Resolve numerical solution U definition of Cole–Hopf transformation. Step 4: Use solution qn þ 1 to replace qn , and repeat Step 2 and 3 for each time steps.
6. Numerical experiments In order to display the capacity of the proposed numerical model for solving the Burgers equation, four numerical experiments are considered. All numerical experiments are compared with analytical solutions for error analysis, and the essential conditions can be obtained by substituting the time and position of boundary points into the analytical solution. The error distribution applied in this research includes the maximum relative error ERel and the L2 norm error EL2 , which are defined as , , , ^ , Max U x i ; t U x i ; t 1 r i r N ; ð37Þ ERel ¼ , , Max U x ; t i 1rirN
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u ^ 2 u N , , , , u∑ U x i ; t U x i ; t u ui ¼ 1 u : EL2 ¼ u 2 N , t , U x ; t ∑ i
ð38Þ
i¼1
^ , Note that U is the analytical solution, U is the numerical solution and N is the number of test points. ,
85
number defines the dominancy between convection and diffusion. Higher Reynolds number leads to a stronger effect of the convection, and vice versa. In this example, when higher Reynolds number is adopted, the solution will be more similar to a step function or a Heaviside function propagating through the domain, making it more difficult to obtain a good numerical solution. The analytical solution of this experiment is given as [38]: 1
; x1 A ½0; 1; t A ð0; 2: 1 þ exp Reð2x1 t Þ=4
6.1. One-dimensional problem with propagating wave
u1 ðx1 ; tÞ ¼
The first example is set to verify the capability of proposed numerical model depicting a propagating wave. Here the Reynolds
Fig. 2(a–c) is the evolution of velocity u1 at different time steps. From Fig. 2(a) we can see when Re ¼ 1, the velocity u1 distributes
1
t=0.0000
u1 950.0 900.0 850.0 800.0 750.0 700.0 650.0 600.0 550.0 500.0 450.0 400.0 350.0 300.0 250.0 200.0 150.0 100.0 50.0
0.8
0.6
x2 0.4
0.2
0
1
0
0.2
0.4
x1
0.6
0.8
t=0.0000
0.6
x2 0.4
0.2
1
0
0.2
0.4
x1
0.6
0.8
0.6
x2 0.4
0.2
0
0.2
0.2
0
1
0.4
x1
0.6
0.8
1
0
0.2
0.4
x1
0.6
0.8
1
t=0.0500
u1 9.5 9.0 8.5 8.0 7.5 7.0 6.5 6.0 5.5 5.0 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0
0.8
0.6
x2 0.4
0.2
u1
0.8
0
0.4
0
0.095 0.090 0.085 0.080 0.075 0.070 0.065 0.060 0.055 0.050 0.045 0.040 0.035 0.030 0.025 0.020 0.015 0.010 0.005
950.0 900.0 850.0 800.0 750.0 700.0 650.0 600.0 550.0 500.0 450.0 400.0 350.0 300.0 250.0 200.0 150.0 100.0 50.0
x2
1
t=0.0000
u1
0.6
u1 9.5 9.0 8.5 8.0 7.5 7.0 6.5 6.0 5.5 5.0 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0
t=0.0005
0.8
1
0.8
0
1
1
0
0.2
0.4
x1
0.6
0.8
1
u1
t=5.0000
0.095 0.090 0.085 0.080 0.075 0.070 0.065 0.060 0.055 0.050 0.045 0.040 0.035 0.030 0.025 0.020 0.015 0.010 0.005
0.8
0.6
x2
0.4
0.2
0
ð39Þ
0
0.2
0.4
x1
0.6
0.8
1
Fig. 7. The evolution of u1 for Example 6.3, (a) Re ¼ 0:01, t ¼ 0, (b) Re ¼ 0:01, t ¼ 0:0005, (c) Re ¼ 1, t ¼ 0, (d) Re ¼ 1, t ¼ 0:05, (e) Re ¼ 100, t ¼ 0, and (f) Re ¼ 100, t ¼5.
86
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92
as a nearly straight line. When Reynolds number increases to 10 or 100 as shown in Fig. 2(b and c), an obvious wave front can be observed propagating through the computational domain. Adopting higher Reynolds number results in steeper wave fronts. Fig. 3(a–f) demonstrates the error evolution of different Reynolds number and error definition. Both ERel and EL2 yield acceptable accuracy, and better accuracy can be achieved by adding more global points. The comparison between ERel and EL2 has also shown that the distributions of the two errors have similar trend. The similarity in ERel and EL2 means the relative errors are evenly distributed throughout the domain, including the steep wave front of physical value u1 . This is an indication that our proposed numerical model is capable for solving steep wave front problems.
1
t=0.0000
-50.0 -100.0 -150.0 -200.0 -250.0 -300.0 -350.0 -400.0 -450.0 -500.0 -550.0 -600.0 -650.0 -700.0 -750.0 -800.0 -850.0 -900.0 -950.0
0.6
x2 0.4
0.2
1
0
0.2
0.4
x1
0.6
0.8
t=0.0000
0.6
x2 0.4
0.2
1
0.2
0.4
x1
0.6
0.8
-0.005 -0.010 -0.015 -0.020 -0.025 -0.030 -0.035 -0.040 -0.045 -0.050 -0.055 -0.060 -0.065 -0.070 -0.075 -0.080 -0.085 -0.090 -0.095
0.6
x2 0.4
0.2
0
0.2
0.4
x1
0.6
0.8
1
u2 -50.0 -100.0 -150.0 -200.0 -250.0 -300.0 -350.0 -400.0 -450.0 -500.0 -550.0 -600.0 -650.0 -700.0 -750.0 -800.0 -850.0 -900.0 -950.0
0.6
x2 0.4
0.2
0
1
0
0.2
0.4
x1
0.6
0.8
1
t=0.0500
u2 -0.5 -1.0 -1.5 -2.0 -2.5 -3.0 -3.5 -4.0 -4.5 -5.0 -5.5 -6.0 -6.5 -7.0 -7.5 -8.0 -8.5 -9.0 -9.5
0.8
0.6
x2 0.4
0.2
0
u2
t=0.0000
t=0.0005
0.8
1
0.8
0
1
u2 -0.5 -1.0 -1.5 -2.0 -2.5 -3.0 -3.5 -4.0 -4.5 -5.0 -5.5 -6.0 -6.5 -7.0 -7.5 -8.0 -8.5 -9.0 -9.5
0
The initial condition is particularly set at t ¼ 0:01, the ending time and space domain are defined differently for different cases.
1
0.8
0
The simulation of diffusive N-wave is another classical Burgers problem, which can be used to test the performance of proposed numerical model which describes strong diffusion phenomena. In this problem the Reynolds number affects the rate of diffusion and the effective range of the diffusive N-wave. The analytical solution is given by [38]: pffiffiffiffiffiffiffiffi 1=t expð Rex21 =4tÞ x1 pffiffiffiffiffiffiffiffi u1 ðx1 ; tÞ ¼ ð40Þ ; x1 A Ω; t 4 0: t 1 þ 1=t expð Rex21 =4tÞ
u2
0.8
0
6.2. One-dimensional problem with diffusive N-wave
1
0
0.2
0.4
x1
0.6
0.8
1
u2
t=5.0000
-0.005 -0.010 -0.015 -0.020 -0.025 -0.030 -0.035 -0.040 -0.045 -0.050 -0.055 -0.060 -0.065 -0.070 -0.075 -0.080 -0.085 -0.090 -0.095
0.8
0.6
x2 0.4
0.2
0
0
0.2
0.4
x1
0.6
0.8
1
Fig. 8. The evolution of u2 for Example 6.3, (a) Re ¼ 0:01, t ¼ 0, (b) Re ¼ 0:01, t ¼ 0:0005, (c) Re ¼ 1, t ¼ 0, (d) Re ¼ 1, t ¼ 0:05, (e) Re ¼ 100, t ¼ 0, and (f) Re ¼ 100, t¼ 5.
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92
Starting with Re ¼ 1, we adopt time range of t A ð0:01; 10 and space domain of x1 A ½ 10; 10. Fig. 4(a) shows clear evolution with diffusive N-wave pattern, and the waveform will flatten and spread throughout the space domain. For Re ¼ 100, the time range of t A ð0:01; 1 and space domain of x1 A ½ 1; 1 is adopted to simulate the diffusive N-wave evolution until the waveform is nearly horizontal. As shown in Fig. 4(b), the tendency of u1 for Re ¼ 100 is very similar to the case of Re ¼ 1, except the scale is now smaller. Therefore, it is no surprise to see the profiles of error
87
evolutions are also similar to those shown in Fig. 5(a–d), where both error distributions resulted with good accuracy. 6.3. Two-dimensional problem in L-shaped irregular domain From former examples the proposed numerical model is available to obtain good solution for 1D Burgers problem, now consider a more complicated problem with two-dimensional system equations in irregular domain. The diagram of L-shaped computational domain and point distribution is considered as Fig. 6, the analytical solutions for
Fig. 9. Error distribution for Example 6.3 of u1 with different number of points in space (a) EL2 of Re ¼ 0:01, Δt ¼ 10 5 , (b) ERel of Re ¼ 0:01, Δt ¼ 10 5 , (c) EL2 of Re ¼ 1, Δt ¼ 10 2 , (d) ERel of Re ¼ 1, Δt ¼ 10 2 , (e) EL2 of Re ¼ 100, Δt ¼ 10 1 , and (f) ERel of Re ¼ 100, Δt ¼10 1.
88
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92
this problem are: 8 10expð 5x1 þ 5x2 þ 50t > Re Þ > < u1 ðx1 ; x2 ; tÞ ¼ Re½1 þ expð 5x1 þ 5x2 þ 50tÞ Re
10expð 5x1 þ 5x2 þ 50t > Re Þ > : u2 ðx1 ; x2 ; tÞ ¼ Re½1 þ expð 5x1 þ 5x2 þ 50t Re Þ
1
;
ðx1 ; x2 Þ A Ω; t 4 0;
ð41Þ
where the ending time is selected independently for different Reynolds numbers in this problem. Fig. 7(a–f) and Fig. 8(a–f) demonstrate the evolution of velocity u1 and u2 with different Reynolds numbers. It can be observed that different Reynolds numbers lead to similar pattern of evolution, only different rates of evolution and magnitude of velocity fields.
t=0.00
1
ERel of u1 1.0E-03 7.0E-04 5.0E-04 2.0E-04 1.0E-04 5.0E-05
0.8
0.6
t=0.00 ERel of u2 1.0E-03 7.0E-04 5.0E-04 2.0E-04 1.0E-04 5.0E-05
0.8
0.6
x2
x2 0.4
0.4
0.2
0.2
0
1
0
0.2
0.4
x1
0.6
0.8
0
1
t=0.05
1
ERel of u1 1.0E-03 7.0E-04 5.0E-04 2.0E-04 1.0E-04 5.0E-05
0.8
0.6
0
0.2
0.4
x1
0.6
1
t=0.05 ERel of u2 1.0E-03 7.0E-04 5.0E-04 2.0E-04 1.0E-04 5.0E-05
0.8
0.6
x2
0.8
x2 0.4
0.4
0.2
0.2
0
1
0
0.2
0.4
x1
0.6
0.8
0
1
t=0.10
1
ERel of u1 1.0E-03 7.0E-04 5.0E-04 2.0E-04 1.0E-04 5.0E-05
0.8
0.6
0
0.2
0.4
x1
0.6
1
t=0.10 ERel of u2 1.0E-03 7.0E-04 5.0E-04 2.0E-04 1.0E-04 5.0E-05
0.8
0.6
x2
0.8
x2 0.4
0.4
0.2
0.2
0
0
0.2
0.4
x1
0.6
0.8
1
0
0
0.2
0.4
x1
0.6
0.8
1
Fig. 10. Space distribution of ERel for Example 6.3 with Re ¼ 1, Δt ¼ 10 3 , 1027 global points and 50 local points (a) ERel of u1 at t ¼ 0, (b) ERel of u2 at t ¼ 0, (c) ERel of u1 at t ¼ 0:05, (d) ERel of u2 at t ¼ 0:05, (e) ERel of u1 at t ¼ 0:1, and (f) ERel of u2 at t ¼0.1.
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92
89
Fig. 11. Error evolution for Example 6.3 of u1 with different time interval (a) EL2 of Re ¼ 1, and (b) ERel of Re¼ 1.
Detailed analysis of the proposed numerical model, including different time intervals, the number of global points, and the number of local points are examined in this example. Fig. 9(a–f) shows the error evolutions for the velocity u1 . From observation, adding global points or expanding local influence domain can achieve better accuracy. By comparing the numerical solutions of adding global points and expanding local influence domain, it is obvious that expanding local influence domain can obtain the same accuracy as adding global points. Next, the relative error distributions in space shown as Fig. 10(a–f), demonstrate both u1 and u2 have the same order of accuracy during different time. Furthermore, Fig. 10(a–f) exhibits that the relative error is higher while the location is closer to the boundary. This behavior is caused by the shape of the local influence area, which is different for each point near the boundaries. This will affect the accuracy due to certain amount of nearest-points which was chosen to depict the local influence area. Moreover the selection of local influence area may also be different for various variables, Fig. 10(c–f) shows the error distributions which are more vertically for u1 and horizontally for u2 , despite the same magnitude and similar space distributions are used for the solutions. Fig. 10 indicates further improvements may be implemented to the proposed model by applying a better method for selecting the local influence area. Finally, Fig. 11(a and b) displays the error distribution of Re ¼ 1 with 1027 global points, showing evidence that the accuracy can be also improved by using smaller time interval. 6.4. Three-dimensional problem in irregular domain For the final example, the three-dimensional Burgers problem in irregular domain is considered. The computational domain of this problem is given as Fig. 12: 8 9 x1 ¼ P sin ψ cos θ; > > < = ∂Ω ¼ ðx1 ; x2 ; x3 Þ x2 ¼ P sin ψ sin θ; 0 o ψ r π ; 0 o θ r 2π ; > > : ; x ¼ P cos ψ ; 3
ð42Þ where ellipsoid radius P is defined as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 P ¼ 3ð cos 3θ þ 8 sin θÞ1=3 :
ð43Þ
The analytical solution for this problem is given by 8 u ðx ; x ; x ; tÞ ¼ > > > 1 1 2 3 < u2 ðx1 ; x2 ; x3 ; tÞ ¼ > > > : u3 ðx1 ; x2 ; x3 ; tÞ ¼
2 expð t=ReÞ cos x1 sin x2 sin x3 Re 1 þ expð t=ReÞ sin x1 sin x2 sin x3 2 expð t=ReÞ sin x1 cos x2 sin x3 Re 1 þ expð t=ReÞ sin x1 sin x2 sin x3 2 expð t=ReÞ sin x1 sin x2 cos x3 Re 1 þ expð t=ReÞ sin x1 sin x2 sin x3
Fig. 12. The irregular domain for Example 6.4, a total number of 2547 points are shown in this figure. The gray shade is the boundary ∂Ω, the black spheres are the boundary points and the red spheres depict the global points. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
In Fig. 13(a–f), the plots show the velocity distribution of u1 ; u2 and u3 on the planes of x1 ¼ 1:5; x2 ¼ 1:5 and x3 ¼ 1:5. Fig. 13(a–c) shows the distinctive N-wave patterns for all three velocity components u1 ; u2 and u3 at t ¼ 0. While time elapses, the diffusive N-wave propagates through the domain, and the diffused solution at t ¼ 1 are shown in Fig. 13(d–f). The error evolutions in Fig. 14(a–d) indicate that for either Re ¼ 1 or Re ¼ 100, the numerical solutions obtained by the proposed model consistently showing acceptable accuracy. The error evolution rises during the late period of experiment due to the fact that the scale of the solution becomes very small. Therefore, even the smallest difference will be magnified in relative error, and the proposed numerical scheme is still able to achieve good accuracy under such circumstances. By further testing the problem with different numbers of global and local points, the coherent behavior of increasing accuracy while adding global points or expanding local influence domain are also demonstrated in error evolutions. 7. Conclusions
; ðx1 ; x2 ; x3 Þ A Ω; t 40:
ð44Þ
The proposed numerical model of FD-LMAPS using the Cole– Hopf transformation has proven the fine capacity for simulation of
90
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92
Fig. 13. Evolution of the velocity for Example 6.4 distributed on planes of x1 ¼ 1:5, x2 ¼ 1:5 and x3 ¼ 1:5. Conditions given as Re ¼ 1, 10,825 global points and 50 local points, (a) u1 at t ¼ 0, (b) u2 at t ¼ 0, (c) u3 at t ¼ 0, (d) u1 at t ¼ 1, (e) u2 at t ¼ 1, and (f) u3 at t ¼ 1.
multi-dimensional Burgers problems in regular and irregular domains with sufficient accuracy. Although with limitation to only solve problems with Dirichlet type boundary conditions, the characteristics of wave propagation and diffusive N-wave have been simulated with result achieving maximum relative error within the order of 10 3 . Four numerical experiments, namely, the one-dimensional problem with wave propagation, the onedimensional problem with diffusive N-wave, the two-dimensional problem in L-shaped irregular domain and the three-dimensional
problem in irregular domain, have been analyzed in this paper. The comparisons between numerical solutions and analytical solutions have given series of results with high accuracy. The error distributions also have shown consistent behavior of achieving better accuracy while adding global points and expanding local influence domain, which can also obtain better efficiency. Hereby we conclude that the proposed numerical model of FD-LMAPS with the Cole–Hopf transformation is a promising tool for future researches and applications.
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92
91
Fig. 14. Error evolution for Example 6.4 of velocity u1 (a) EL2 of Re ¼ 1, Δt ¼ 10 2 , (b) ERel of Re ¼ 1, Δt ¼ 10 2 , (c) EL2 of Re ¼ 100, Δt ¼ 10 2 , and (d) ERel of Re¼ 100, Δt ¼10 2.
Symbol
Appendix A
All variables in this paper are formatted in italic or regular style, only the vectors and matrices are formatted in arrow or bold style. Superscripts and subscripts
Detail derivations of Cole–Hopf transformation for multidimensional Burgers equations. The Burgers equation
b d i ib j k m n s
indicates the boundary. number of dimensions. ith point of global domain. ibth point of boundary. jth point of global domain. kth point of local domain. mth point of local domain. indicates nth time step. indicates the initial time step.
Acknowledgment Financial supports of this paper are provided by the National Science Counsel of Taiwan under the grant No. NSC100-3332-E002–020. It is greatly appreciated. The fourth author acknowledges the support of Distinguished Overseas Visiting Scholar Fellowship provided by the Ministry of Education of China.
, , , , ∂U 1 þ U U∇ U ¼ ∇2 U : Re ∂t
The Cole–Hopf transformation ,
Uqþ
2 ∇q ¼ 0andq a 0: Re
Cole-Hopf transformation for various differential terms 8 2 ∂q > udi ¼ Req > ∂xdi > > > 2 > ∂udi ∂q ∂q 2 ∂ q 2 > ¼ > > ∂t Req ∂xdi ∂t þ Req2 ∂xdi ∂t > > > 2 > 2 > ∂udi ∂q 2 ∂ q 2 > ¼ Req þ Req > 2 ∂x > ∂x2di di < ∂xdi : ∂udi ∂2 q ∂q ∂q 2 2 ¼ þ > Req ∂xdi ∂xdj Req2 ∂xdi ∂xdj > > ∂xdj > 3 > > 2 3 ∂q ∂2 q ∂q > 2 ∂ q 6 4 > ∂ u2di ¼ Req þ Req > 2 ∂x 2 Req3 ∂x ∂x3di ∂xdi > di ∂xdi di > > > 2 > > ∂2 udi ∂3 q ∂2 q ∂q ∂q ∂2 q ∂q ∂q 2 4 2 4 > > : ∂x2 ¼ Req ∂x x2 þ Req2 ∂xdi ∂xdj ∂xdj þ Req2 ∂xdi ∂x2 Req3 ∂xdi ∂xdj dj
di dj
dj
subscript di and dj depicts different dimension number d:
92
C.Y. Lin et al. / Engineering Analysis with Boundary Elements 40 (2014) 78–92
Substituting the differential terms after Cole–Hopf transformation into the Burgers equation and obtain ∇q ∂q 1 2 ∂q 1 2 ∇ q ∇ ∇ q ¼ 0: q ∂t Re ∂t Re Using 1-D Burgers equation as an example ∂u1 ∂u1 ∂t þ u1 ∂x1 -
2
1 ∂ u1 ¼ Re ∂x2 1
" 2 # 2 ∂ q 2 ∂q ∂q 2 ∂q 2 ∂2 q 2 ∂q þ þ þ Req ∂x1 ∂t Req2 ∂x1 ∂t Req ∂x1 Req ∂x21 Req2 ∂x1 2
" 3 # 1 2 ∂3 q 6 ∂q ∂2 q 4 ∂q þ Re Req ∂x31 Req2 ∂x1 ∂x21 Req3 ∂x1 " 2 # ∂2 q 1 ∂q ∂q ∂q 2 ∂2 q 2 ∂q þ þ þ - ∂x1 ∂t q ∂x1 ∂t ∂x1 Req ∂x21 Req2 ∂x1 3 1 ∂3 q 3 ∂q ∂2 q 2 ∂q þ ¼ 3 2 2 Re ∂x1 Req ∂x1 ∂x1 Req ∂x1 3 ∂2 q 1 ∂q ∂q 2 ∂q ∂2 q 2 ∂q þ þ - ∂x1 ∂t q ∂x1 ∂t Req ∂x1 ∂x21 Req2 ∂x1 3 1 ∂3 q 3 ∂q ∂2 q 2 ∂q ¼ þ Re ∂x31 Req ∂x1 ∂x21 Req2 ∂x1 ¼
∂2 q 1 ∂q ∂q 1 ∂3 q 1 ∂q ∂2 q ¼ þ þ 3 ∂x1 ∂t q ∂x1 ∂t Re ∂x1 Req ∂x1 ∂x21 ! ! 1 ∂q ∂q 1 ∂2 q ∂ ∂q 1 ∂2 q ¼ 0: q ∂x1 ∂t Re ∂x21 ∂x1 ∂t Re ∂x21 -
It is obvious that diffusion equation satisfies the Burgers equations transformed via Cole–Hopf transformation ∂q 1 2 ∇ q ¼ 0: ∂t Re Therefore diffusion equation is applied as the governing equation after Cole–Hopf transformation in this study. References [1] Burgers JM. A mathematical model illustrating the theory of turbulence. Adv Appl Mech 1948;1:171–99. [2] Yao G, Chen CS, Kolibal J. A localized approach for the method of approximate particular solutions. Comput Math Appl 2011;61:2376–87. [3] Bahadir AR. A fully implicit finite-difference scheme for two-dimensional Burgers' equations. Appl Math Comput 2003;137:131–7. [4] Fu S, Hodges B. Time-centered split method for implicit discretization of unsteady advection problems. J Eng Mech 2009;135:256–64. [5] Young DL, Fan CM, Hu SP, Atluri SN. The Eulerian–Lagrangian method of fundamental solutions for two-dimensional unsteady Burgers' equations. Eng Anal Bound Elem 2008;32:395–412. [6] Hopf E. The partial differential equation ut þ uux ¼ μxx. Commun Pure Appl Math 1950;3:201–30. [7] Cole JD. On a quasilinear parabolic equation occurring in aerodynamics. Q Appl Math 1951;9:225–36. [8] HSL. A collection of FORTRAN codes for large scale scientific computation, 〈http://www.hsl.rl.ac.uk〉. [9] Kutluay S, Bahadir AR, Özdeş A. Numerical solution of one-dimensional Burgers equation: explicit and exact-explicit finite difference methods. J Comput Appl Math 1999;103:251–61.
[10] Kadalbajoo MK, Awasthi A. A numerical method based on Crank–Nicolson scheme for Burgers' equation. Appl Math Comput 2006;182:1430–42. [11] Öziş T, Aksan EN, Özdeş A. A finite element approach for solution of Burgers' equation. Appl Math Comput 2003;139:417–28. [12] Varoǧ lu E, Finn WD. Space-time finite element incorporating characteristics for the Burgers' equation. Int J Numer Methods Eng 1980;16:171–84. [13] Chino E, Tosaka N. Dual reciprocity boundary element analysis of timeindependent Burgers' equation. Eng Anal Bound Elem 1998;21:261–70. [14] Kakuda K, Tosaka N. The generalized boundary element approach to Burgers' equation. Int J Numer Methods Eng 1990;29:245–61. [15] Chen JK, Beraun JE. A generalized smoothed particle hydrodynamics method for nonlinear dynamic problems. Comput Method Appl M 2000;190:225–39. [16] Golberg MA, Chen CS. Improved multiquadric approximation for partial differential equations. Eng Anal Bound Elem 1996;18:9–17. [17] Hon YC, Mao XZ. An efficient numerical scheme for Burgers' equation. Appl Math Comput 1998;95:37–50. [18] Cheng AHD, Golberg MA, Kansas EJ, Zammito G. Exponential convergence and H-c multiquadric collocation method for partial differential equations. Numer Methods Partial Differ Equ 2003;19:571–94. [19] Li JC, Hon YC, Chen CS. Numerical comparisons of two meshless methods using radial basis functions. Eng Anal Bound Elem 2002;26:205–25. [20] Chen R, Wu Z. Applying multiquadric quasi-interpolation to solve Burgers' equation. Appl Math Comput 2006;172:472–84. [21] Young DL, Hu SP, Fan CM, Chen CW. Method of fundamental solutions and Cole–Hopf transformation for Burgers' equations. Comput Fluid Dyn J 2006;14:454–67. [22] Lin CY, Gu MH, Young DL. The time-marching method of fundamental solutions for multi-dimensional telegraph equations. CMC-Comput Mater Con 2010;18:43–68. [23] Tsai CC, Chen CS, Hsu TW. The method of particular solutions for solving axisymmetric polyharmonic and poly-Helmholtz equations. Eng Anal Bound Elem 2009;33:1396–402. [24] Wen PH, Chen CS. The method of particular solutions for solving scalar wave equations. Int J Numer Method Biomed Eng 2010;26:1878–89. [25] Atluri SN, Zhu T. A new meshless local Petrov–Galerkin (MLPG) approach in computational mechanics. Comput Mech 1998;22:117–27. [26] Lin H, Atluri SN. The meshless local Petrov–Galerkin (MLPG) method for solving incompressible Navier–Stokes equations. CMES-Comput Model Eng 2001;2:117–42. [27] Lee CK, Liu X, Fan SC. Local multiquadric approximation for solving boundary value problems. Comput Mech 2003;30:396–409. [28] Shu C, Ding H, Yeo KS. Computation of incompressible Navier–Stokes equations by local RBF-based differential quadrature method. CMES-Comput Model Eng 2005;7:195–205. [29] Zhou CH, Shu C. A local domain-free discretization method for simulation of incompressible flows over moving bodies. Int J Numer Meth Fl 2011;66:162–82. [30] Kosec G, Saler B. H-adaptive local radial basis function collocation meshless method. CMC-Comput Mater Con 2011;26:227–53. [31] Shu C, Ding H, Yeo KS. Local radial basis function-based differential quadrature method and its application to solve two dimensional incompressible Navier– Stokes equations. Comput Methods Appl Mech Eng 2003;192:941–54. [32] Yao G, Šarler B, Chen CS. A comparison of three explicit local meshless methods using radial basis functions. Eng Anal Bound Elem 2011;35:600–9. [33] Karageorghis A, Chen CS, Smyrlis YS. Matrix decomposition RBF algorithm for solving 3D elliptic problems. Eng Anal Bound Elem 2009;33:1368–73. [34] Gu MH, Fan CM, Young DL. The method of fundamental solutions for the multi-dimensional wave equations. J Mar Sci Tech-Taiw 2011;19:586–95. [35] Rippa S. An algorithm for selecting a good value for the parameter c in radial basis function interpolation. Adv Comput Math 1999;11:193–210. [36] Wertz J, Kansas EJ, Ling L. The role of the multiquadric shape parameters in solving elliptic partial differential equations. Comput Math Appl 2006;51:1335–48. [37] Chen CS, Fan CM, Wen PH. The method of particular solutions for solving certain partial differential equations. Numer Methods Partial Differ Equ 2012;28:506–22. [38] Polyanin AD, Zaitsev VF. Handbook of nonlinear partial differential equations. 1st ed.. Florida: Chapman and Hall/CRC; 2004.