Strojniški vestnik - Journal of Mechanical Engineering 61(2015)4, 254-264 © 2015 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2014.2185
Original Scientific Paper
Received for review: 2014-09-19 Received revised form: 2015-03-03 Accepted for publication: 2015-03-11
BEM-Based Algorithm for URANS Simulations of Flow over a Square Cylinder
Janez Lupse* - Leopold Skerget - Jure Ravnik
University of Maribor, Faculty of Mechanical Engineering, Slovenia
The boundary element method (BEM) is extended for the discretisation of the incompressible Navier-Stokes equations in the velocity-vorticity form for the URANS formulation. BEM uses fundamental solutions as weighting functions, thus the approach enables inclusions of some of the governing physics in the method itself. For testing purposes the flow around the square cylinder in the two-dimensional channel with a blockage ratio of 1/8 was chosen as the appropriate test case. The flow was calculated for the Reynolds numbers ranging from 100 to 1000. At Reynolds number 1000, the governing equations were solved directly in a quasi DNS manner and in the RANS form using the Spalart-Allmaras turbulence model within an URANS simulation.
Keywords: boundary element method, velocity-vorticity form, incompressible flow, URANS, square cylinder Highlights
•	Extension of the boundary element based algorithm is proposed for simulation of turbulent flows.
•	Velocity vorticty formulation of Navier - Stokes equations is employed.
•	Simulation of flow over a square cylinder inside a channel at blockage of 1/8.
•	Simulations were performed within the unsteady laminar regime of the flow and within the low turbulent regime of the flow.
•	Comparison between the unsteady RANS model and the quasi-DNS simulation.
0 INTRODUCTION
The boundary element method has been seen as a very successful method for solving various problems. In general, this method has an advantage over other numerical methods when the fundamental solution to the governing equations exists, thus effectively reducing the dimension of the problem [1]. Using the description of the governing physics regarding the fundamental solution, a stable and highly accurate numerical scheme is obtained, even on coarse numerical meshes. Unlike most of the other discretisation methods, using fundamental solutions this method relies on capturing the governing physics instead of the brute force. Furthermore, the derivatives of the transported variables, normal to the boundary of the domain or sub-domain, are a part of the solution of the system. This enables accurate representation of these values and the corresponding integral parameters such as the Nusselt number for example. Unfortunately, the governing equations for the fluid flow in general are of convective-diffusive type with important and non-trivial source terms. Even when the convective-diffusive fundamental solution is used, discretisation of the whole solution domain is needed. In addition, when a fundamental solution that includes more of the governing physics of the flow is used (with the exception of the Laplace fundamental solution), the integrals that arise from the discretisation of the governing equations need to
be re-evaluated every time the corresponding field of variables (e.g., velocity) is changed. However, the inclusion of the flow physics into the discretisation usually leads to coarser numerical meshes.
Whilst there are various approaches for the discretisation of domain terms (dual reciprocity method [2], the green element method [3], etc.) arising regarding the discretisations of the Navier-Stokes equations the so-called boundary-domain integral method (BDIM) [4] is used in our case. This approach gives discretisation of the solution domain similar to the finite element method, with each cell representing a separate sub-domain. The appropriate compatibility and equilibrium conditions have to be applied in order to interconnect the sub-domains. This results in an over-determined linear system of equations with a non-symmetric sparse system matrix that is solved by the LSQR solver [5] using diagonal-preconditioning.
In Lupse et al. [6] the method was extended for solving the Reynolds averaged Navier-Stokes equations (RANS). This article is aimed at testing the developed algorithm for the intrinsically unsteady flows that can also be described using the unsteady RANS (URANS) approach. In regard to this purpose the flow over a square cylinder inside a channel was chosen as a relevant test case.
The confined flow around the square cylinder is a challenging test for numerical methods as a multitude of interesting flow phenomena are present within relatively simple geometry. This type of flow
254
*Corr. Author's Address: University of Maribor, Faculty of Mechanical Engineering, Smetanova 17, 2000 Maribor, Slovenia, janez.lupse@yahoo.ca
Strojniški vestnik - Journal of Mechanical Engineering 61(2015)4, 254-264
is important during many engineering applications ranging from the cooling of electronics to building aerodynamics. Although free-flow around the circular cylinders has been the focus of extensive research [7], less effort has been put into researching the confined flow around bluff bodies, especially around the square cylinders. Additional parameters need to be defined for the confined flow, namely those inflow velocity profiles and blockage ratios of the channel that do not appear in free-flow cases. Davis et al. [8] performed experiments and numerical simulations of the flows for two different blockage ratios (B = 1/4 and B = 1/6) and a wide variety of Reynolds numbers ranging from 100 and up to 1800. Later, [9] performed numerical simulations for blockage ratio of 1/8 and Reynolds numbers of the flow up to 900. A common argument about simulations of [9] is that the meshes used were too coarse around the vicinity of the cylindrical obstacle, hence the discrepancies of the obtained results with other authors. Breuer et al. [10] tested two different numerical methods on this type of flow at the blockage ratio of 1/8 for Reynolds numbers up to 300. Their results are widely considered as the benchmark for the laminar part of the flow.
The results of [10] were later confirmed by the so-called multi-particle collision dynamics method tested at Reynolds numbers lower than 130 [11]. Galletti et al. [12] performed a test of a proper orthogonal decomposition-based model on the confined flow over a square cylinder and obtained comparable results to [10] having a blockage ratio of 1/8. Later the flow was also simulated for various laminar Reynolds numbers and extended those simulations to different blockage ratios and non-isothermal flows [13] and [14].
The flow chosen for the test case featured many interesting phenomena at different flow regimes. At the Reynolds number equal to 1 creeping flow is observed. The flow separates from the trailing edge of the cylinder with any increase in the Reynolds number. Although two recirculation bubbles are present and increase in length as the Reynolds number of the flow increases, the flow is steady up to a critical point. Authors in literature have reported different values for this critical number but agree that it lies somewhere between Re = 54 [15] and Re = 70 [16]. After this point the flow becomes unsteady but periodic. Separation occurs at the trailing edge at first and with increasing Reynolds number of the flow moves to the leading edge of the cylinder. Flow separation causes a distinct pattern of flow, known as the Von Karman vortex street. Spatial effects occur at about the Reynolds number 300 [10]. As the Reynolds number of the flow is defined by the cylinder diameter
and not the channel width it is necessary to have in mind that with a used blockage ratio of 1/8 at about this Reynolds number, transition to turbulence in the channel flow is slowly starting to occur. Thus in reality, the flow beyond the Reynolds number 300 is be-coming increasingly turbulent. Consequentially, for the inlet profile the turbulent velocity profile would be more appropriate passed these Reynolds numbers of the flow. The data available in literature for the Reynolds numbers up to 1000 [8] is obtained by imposing the standard parabolic velocity profile. As such, the parabolic velocity profile in this work, used at the inlet boundary condition, was used purely for the testing purposes of the developed numerical algorithm and code in order to enable easier and more accurate comparison with the available data of [9].
1 GOVERNING EQUATIONS
Whilst the more common approach to describing the fluid flow is the use of the Navier-Stokes equations in their primitive forms, the velocity-vorticity formulation is used in our work. In this formulation the flow equations are split into kinematic and kinetic parts. In general, the approach gives us the advantage of eliminating the pressure gradient from the governing equations. However, the boundary conditions for the vorticity field in Eqs. (1) and (2) are generally unknown intuitively, unlike the velocity boundary condition in the momentum transport equation.
The flow kinetics is described by the vorticity transport equation obtained from the momentum transport equation by applying the curl operator [17]. For the incompressible planar flow the following equation is obtained:
da
^+U
dt ' dx
= vn
d a
dxjdxj
1 dpgi
Po
dx,
(1)
where ra is used for the vorticity, Uj for the components of the velocity vector and Xj the spatial coordinates; v0 the kinematic viscosity, p0 the fluid density and g the components of the gravity acceleration vector. The vorticity vector is treated as a scalar within the planar flow since only the component of the vorticity normal to the plane of the flow has a non-zero value. In the RANS form of Eq. (1), an additional source term is obtained on the right hand side and the RANS form of the vorticity transport equation is then [6]:
da TT da
— + Uj-= v,
dt 1 dxj
d2 a
i dpgi i /
dxi dxi p0 dxi p0 dxi
(2)
BEM-Based Algorithm for URANS Simulations of Flow over a Square Cylinder
255
Strojniski vestnik - Journal of Mechanical Engineering 61(2015)4, 254-264
where the additional source termfm is defined by:
e dvtœ dvt r 2 dvt dUi Ji = -ßij-+ 2eij-m + 2--. (3 )
dxj	dxj dxj dxj
and constants of the used model see the original reference [20].
2 BOUNDARY-DOMAIN INTEGRAL EQUATIONS
The effective viscosity is now comprised of the molecular part v0 and the modelled part vt . The overbar indicating the averaging process is omitted for the sake of brevity.
The kinematics part is formed by the Poisson equation that links the velocity and vorticity fields. It is derived at from the continuity equation and the vorticity definition [1]. For the planar flow the following form of the kinematics equation is obtained:
d2Ui dm n
-— + eij-= 0.
dx;dx; dxh
(4)
j j
k
Prior to discretisation using the boundary element method, the governing equations need to be rewritten in integral form. If we use Q for the solution domain and r for its boundary, then the integral form of the kinematics Eq. (4) can be written as
c m u du*d T-\itu*d r
ea [mnu dr — e i G)—— dQ.
J J J	'J J
Q dxJ
(8)
There exist very efficient algorithms for solving the velocity-vorticity forms of the equations [18] but they are only applicable for special cases using simplified geometry. Care has to be taken however since Eq. (4) admits solutions in which neither solenoidality of the velocity field nor vorticity definition hold true [19]. In order to obtain a general solution, an additional equation has to be included [1]:
dUi	dUj
-+ ei1n ¡m = -e3-
dn 3 3 1 dt
(5)
where the normal and tangential derivatives in respect to domain boundary are defined as:
dn
DX :
(6)
where we use Ç for the location of the source point, c for the geometric coefficient and u* to denote the Laplace fundamental solution.
* 1 1
u = — ln —. 2n r
(9)
In Eq. (9) r is used for the distance between the integration and source points. The detailed derivation of the integral form of Eq. (8) can be found in [21]. Eq. (8) is then used for calculation of the velocity field inside the solution domain. In order to calculate the boundary conditions for the vorticity transport equation in a general case Eq. (4) has to be combined with Eq. (5). A new integral form of the kinematics Eq. (4) is obtained;
c OU (^w U f* 'r:
and n and t are the unit normal and unit tangent vectors.
Finally, in order to close the system of equations, the Spalart-Allmaras model [20] is used as a closure model. The main equation of this closure model without the trip term is:
D = cH [1 - Ft 2 JL
Dt	a dx ;
i ~\ dv (Vo .
j /


a
b2
dx;
Cw1F" '
C
b1
21 Ft2
(7)
Eq. (7) is the transport equation for the so-called modified turbulent viscosity from which modelled viscosity vt, as used in Eq. (2) is obtained by simple algebraic expression. For the details on the functions
-ejjUj ^dF-eJ^^ (10)
a dxj
with the tangential derivative defined as in Eq. (6).
The vorticity transport equation is a non-stationary convective-diffusive equation with nonlinear source terms. As such, several fundamental solutions can be used as weighting functions. If the parabolic diffusion fundamental solution is used the following integral form of the vorticity transport equation is obtained [6];
du	dco
c(E^miE,, tk)+ [ffl-dr- [— u dT =
V W V k! r dn r dn

'0 n
= n
2
n
256
Lupse, J. - Skerget, L. - Ravnlk, J.
Strojniški vestnik - Journal of Mechanical Engineering 61(2015)4, 254-264
J mUjnjU*dT - J mU
du
\
Vr
a dxj J
voPo
Jpginju*d r-Jpgi
du
d a
J
* \
Vr
(
3X;
J J
J fmnju"d r-J f
du
d a
* \ da
Vr
a dxJ J
(11)
where I is used for the model source terms. For the Spalart-Allmaras model it is defined as:
I = C
^ dV
b2
dx,
\ J /
+ Cbi [1 - Ft 2 ]v-
C F - Cb1 F ^w1rw 2 rt2 K
V
d
(14)
The index k is used for denoting the current time step, thus tk denotes the current time. Instead of the parabolic diffusion fundamental solution the convection-diffusion fundamental solution can also be used. The following form of the equation is obtained [6];
c tk )+\wQfd r-f|Vd r =
-1 -2 )u"d
"o n
-V-\œ{Uj -Uj0)rijU dT-vo r
±HUj -Ujo ))dn
voPo
Jpgirju*d r-Jpgi
du
\
Vr
dx,-j y
J /rijUd r-J /
du
d n
* \ d n
Vr
n dxj y
(12)
where we use Pj = -2/At, p2 = -1/(2A) and Q* = du*/dn + (Uj0nj0/v0)u*, for the time discretisation. As the transport equation of the turbulence model used is of the same type as the vorticity transport equation, a similar integral equation is derived at:
[c (Ms, tk + d r-jfVd r =
O f	*	O f *
h--I vk-luk-ldQh--I Iu dQ-
O
"o Q
(
"o Q
du
\
¡vUnUdT-¡vU.—dQ
J	O dx,
Q
Vr
O
dv
jvt-—nu
Vr dxj
dr- (vt dq
Q dxj dxj .
(13)
In the final discrete form of the equation the source term I is split into the production and destruction terms and linearised to improve the numerical stability.
3 THE NUMERICAL ALGORITHM 3.1 Discrete Form
In order to obtain the solutions to the governing equations, the solution domain is divided into C domain cells also called sub-domains comprised of 9 nodes (Fig. 1) and E boundary elements;

(15)
C=1
e=1
The classic single domain BEM is used for the calculation of boundary vorticity values by Eq. (10). However, all other equations are discretised using the domain decomposition approach, also called subdomain BEM. The basic idea behind this approach is to split the solution domain into multiple sub-domains, similar to the cells used in the finite element methods. The sub-domains are then interconnected by the appropriate compatibility and equilibrium conditions corresponding to each sub-domain. This approach leads to an overdetermined sparse system matrix and in comparison with the classical BEM enables large savings in CPU time and computer memory.
Fig. 1. Sub-domain:the dot represents the function node and x the double node (node with function and its normal derivative as an unknown)
1
v
0
2
0
1
—e
o
o
v
o
BEM-Based Algorithm for URANS Simulations of Flow over a Square Cylinder
257
Strojniski vestnik - Journal of Mechanical Engineering 61(2015)4, 254-264
As both the vorticity transport and the equations of the turbulence models are of convective-diffusive type, the general discrete form is presented in Eq. (16) where 9 is used for the corresponding field function (vorticity, turbulent kinetic energy, etc.). The number of nodes in each element is marked with n and the general source term S is introduced that includes the sources obtained from time discretisation and various body source terms depending on the equation solved (fr for the vorticity transport equation, I for the turbulence model etc.)
c(tm, t)+t a }T {t}"-t {ge }T & j" =
e=1	e=1	Idn
-t Rj } j + t {dcj } M" +
1
e=1
c=1
+C j f "Ctidcj f -b fl"
-tb r
c=1
(16)
The integrals introduced in Eq. (16) are defined
as:
he = f®n d^r, ge = \xnu'dr, r on	r
d *
aetj = f ®nnju*dr, dCtJ = f Tn —dn,
' re	' nc dXj
bc = f t nudn,	(17)
with the functions x, ® and ¥ in our case representing the constant, quadratic and biquadratic interpolations.
3.2 The Solution Procedure
A brief schematic of the algorithm is presented in Fig. 2. All equations are solved sequentially within a global nonlinear loop - depending whether the turbulence model is included the loop terminates after the solution of the kinetic equation for the domain vorticity values or after solutions of equations of the turbulence model. The main convergence criterion is the convergence of the vorticity field with additional monitoring of the convergence of the turbulence model's values. As the kinetics equations are highly nonlinear, an under relaxation procedure is introduced to enable convergence of the solution. When the
convergence of the solution meets the predefined accuracy the next time-step is initiated and the solution process repeated.
In regard to the velocity-vorticity formulation of the Navier-Stokes equations, boundary conditions for the vorticity transport Eq. (2) need to be supplied. They are obtained from the proper velocity boundary conditions and an initial vorticity field by solving Eq. (10). Due to the unknown initial value of the vorticity field, the solution of Eq. (10) has to be coupled with the vorticity transport equation (Eqs. (11) or (12)) in an iterative way. The Laplace fundamental solution is used in Eq. (10), as the equation is of the Poisson type. Unfortunately, in order to preserve the solenoidality of the velocity field, the equation needs to be solved in a global sense and its discretisation leading to a fully populated system matrix. It is possible to reduce the memory requirements of the algorithm by introducing fast wavelet, or some other similar method for the matrix compression [22] with negligible loss of accuracy.
Fig. 2. Flowchart of the numerical algorithm for the solution of governing equations in the velocity-vorticity form
With the boundary values of the vorticity field known, Eq. (8) is solved for the domain velocity field. As Eq. (8) is the Poisson equation for the velocity field, the simple Laplace fundamental solution is used. The sub-domain integral method is used ([4] and [6]) for the discretisation of Eq. (8) and all the following Eqs. (11) or (12) and (13). This approach yields a rectangular system matrix, solved using the LSQR type solver and a diagonal pre-conditioner [5]. In general the sub-domains can be comprised of various numbers of nodes, depending on the interpolation used for the field function values. In our case, the sub-domains used contained nine mesh nodes, thus
0
n
258
Lupse, J. - Skerget, L. - Ravnlk, J.
Strojniški vestnik - Journal of Mechanical Engineering 61(2015)4, 254-264
enabling quadratic interpolation of the function values. The normal derivative (flux) values that are also part of the solution are interpolated with constant elements, so the same mesh nodes can be used as for the function values (see Fig. 1). Any higher order interpolation of the flux values would require either the inclusion of additional nodes or evaluation of the corner values of fluxes within each sub-domain, which is not recommended as the procedure can introduce additional errors.
The next step in the algorithm is the solving of the vorticity transport equations (Eqs. (11) and (12)) that provide the domain vorticity field values and thus complete the information needed for describing the laminar or (in the case of direct numerical simulation) the turbulent isothermal incompressible flow also. The algorithm iterates the steps described until satisfactory convergence within the vorticity field is achieved. Whilst in the first two parts of the algorithm (boundary vorticity values and domain velocity field calculations), the Laplace fundamental solution is used because of the mathematical simplicities of the equations in question, thus more advanced fundamental solutions can be used for the derivation of the integral statements of the vorticity transport equation. In our case the parabolic diffusion and the
convection-diffusion fundamental solutions were implemented, as can be seen from Eqs (11) and (12).
In cases of non-isothermal flows or the RANS/ URANS simulations, additional equations for the conservation of the energy and/or the turbulence models need to be solved. Those equations are of the same type as the vorticity transport equation and are thus solved in a similar fashion. It also has to be noted that when the low-Reynolds turbulence models are used, a classical linearisation procedure of the dissipation terms in the model's equations is used.
4 DESCRIPTION OF THE TEST CASE
The geometry of the test problem consisted of a plane channel with a square obstacle placed inside. For the unit length, the obstacle diameter D was taken (see Fig. 3), thus we defined the Reynolds number of the test case as Re = DUc / v o, where Uc was the centreline velocity at the inlet section. The blockage ratio of the test problem was defined by: B = D / h = 1/8 = 0.125, where h was the channel's height.
At the solid walls and on the surface of the square cylinder, no-slip boundary condition was prescribed. In order to enable better comparison with other authors [8] and [10], the parabolic velocity profile was

a)
		t	
Inlet	h	□ ID X ' Xvx=vy=0	Outlet
			
b)
X
c)
Fig. 3. a) Geometry and boundary conditions of the test case, b) Mesh c: around the cylinder and c) Mesh c: zoom in; for finest numerical mesh used for the calculation, comprised of approximately 14700 sub-domains and 56342 mesh nodes (mesh lines connect the function
nodes of sub-domains as presented in Fig. 1)
BEM-Based Algorithm for URANS Simulations of Flow over a Square Cylinder
259
Strojniski vestnik - Journal of Mechanical Engineering 61(2015)4, 254-264
prescribed at the inlet. The convective type boundary condition (as described in [23]) was used at the outlet for the velocity field whilst for all other variables the zero normal derivative was prescribed.
The flow was firstly calculated at the Reynolds numbers 100 and 150 on three different meshes, with the focus on increasing mesh density within the vicinity of the obstacle (see Fig. 3).
The coarsest mesh (mesh a) was comprised of around 4700 sub-domains (19132 mesh nodes), the mid-sized mesh (mesh b) had about 10700 subdomains (43586 mesh nodes) and the finest one (mesh c) had about 14700 sub-domains (56342 mesh nodes). As the finest mesh had the greatest resolution within the vicinity of the cylinder, only this one was used in the calculation of flow at the Reynolds number 1000. The time-step used for simulations was At = 2.5 *10-3 and At = 2*10-2 for URANS simulations, normalised as t* = (t U) / D. t* is non-dimensional time, t real time and U average velocity at the inlet of the channel.
4.1 Laminar Flow
Firstly the flow was calculated at the Reynolds numbers 100 and 150. For the initial values of the higher Reynolds simulation, the values were used as obtained by the Re = 100 simulation. Vortex shedding was initiated without introducing any artificial perturbations. Streamlines for quarters of the time period are shown in Fig. 4 for both flow cases. At the lower laminar Reynolds number the separation of the flow occurred at the cylinder trailing edge (see Fig. 4a). However, at the Reynolds number 150, the flow was also separated from the leading edges, as can be seen in Fig. 4b. The separations in both cases then caused periodic vortex shedding, also known as the Von Karman vortex street.
In Figs 5a and b the time variations of the drag and lift coefficients, evaluated from the forces acting on the cylinder are presented and the periodical nature of the flow can be clearly seen. The vortex shedding frequency f was evaluated from the lift coefficient
a)
b)
Fig. 4. Streamlines of the flow at a) Re = 100 and b) Re = 150
a)
Time step	y	Time step
Fig. 5. Time variation of the lift and drag coefficients for the Reynolds number; a) 100 and b) 150
260
Lupse, J. - Skerget, L. - Ravnlk, J.
Strojniški vestnik - Journal of Mechanical Engineering 61(2015)4, 254-264
time series and used for the calculation of the Strouhal number, we defined by:
St =
f ■ D
(18)
Table 1 presents the results for both the Strouhal number and the drag coefficient, together with comparison with available data from literature. General increases for both integral parameters were observed with the increase in the Reynolds number. When taken into account the results from the literature are quite widespread, the calculated Strouhal number matched the given data quite well, whilst the drag coefficient value seems to be systematically lower than the values from the literature.
A detailed comparison between the velocity profiles at the Reynolds number 100 was made with the results of Breuer et al. [10]. In order to make this comparison possible, the velocity profiles were
extracted at approximately the same moment as described in the cited authors work; when the cross-stream velocity Uy at an axial position of 10 cylinder diameters downstream the cylinder, changed its sign from negative to positive. The comparisons between the velocity profiles are shown in Figs. 6a and b, where good agreement can be seen for meshes b and c used for the calculation of the flow, with minor deviations of the profiles towards the outlet section of the channel. Those deviations could most likely be attributed to the usage of relatively coarse mesh far away from the obstacle and a shorter computational domain as compared to [10].
Additional error was also produced by the finite length of time-step thus enabling only approximate determination of the time for the extraction of the velocity profiles.
The comparisons between the velocity profiles at the Reynolds number 150 are shown in Figs. 6c and
b) vy; Re = 100
c) vx; Re = 150
d) vy; Re = 150
Fig. 6. Instantaneous velocity profiles for Re = 100 and Re = 150 at the middle of the channel at a certain moment (see text for explanation), compared with the results of [9]
BEM-Based Algorithm for URANS Simulations of Flow over a Square Cylinder
261
Strojniski vestnik - Journal of Mechanical Engineering 61(2015)4, 254-264
d. As no data for the profiles could be found in the literature, the comparison was only made between meshes b and c.
At the higher Reynolds number the velocity profile was qualitatively similar to the profiles at the Reynolds number 100, with the higher maximum values within the velocity field.
Table 1. Comparison between the Strouhai number and drag coefficient for Reynolds 100 and 150
Re	Mesh	St	Cd
	Breuer[10]	0.139	1.35
	Lamura [11]	0.133	1.38
	Turki [13] and [14]	0.134	1.40
100	Galletti [12]	0.138	-
	Mesh a	0.129	1.33
	Mesh b	0.135	1.31
	Mesh c	0.137	1.27
	Breuer[10]	0.147	1.33
	Turki [13] and [14]	0.136	1.38
150	Galletti [12]	0.144	-
	Mesh a	0.133	1.43
	Mesh b	0.141	1.39
	Mesh c	0.147	1.29
200	Mesh c	0.151	1.37
250	Mesh c	0.155	1.46
4.2 Turbulent Flow
In order to obtain initial values for the flow field, the calculations were initiated at the Reynolds number 500 on the most refined grid for a few hundred timesteps in order to provide initial values for higher Reynolds calculations. The Reynolds number was then increased to the final value of Re = 1000. As we have defined the Reynolds number in the same manner as in the reference literature (with obstacle width and maximum inlet velocity), it has to be noted that the Reynolds number based on the channel width and bulk velocity is around 5300. Although physically less appropriate, the parabolic velocity profile was used at the inlet in order to enable easier and more accurate comparisons between the experimental and numerical data of [8]. The flow was calculated for direct solutions of the governing equations (dubbed the quasi-direct numerical simulation (qDNS) due to the mesh being insufficiently refined) and the modified governing equations in the URANS mode.
In Figs. 7a and b the time series for the drag and lift coefficient are presented for the direct solution of the governing equations for both the parabolic-diffusive and the convective-diffusive fundamental solutions. From these series, the loss of periodicity
12000 13000 14000
Time step c) URANS; Drag and lift
Fig. 7. Time series of lift and drag coefficient at Re = 1000 for qDNS and URANS simulations
262
Lupse, J. - Skerget, L. - Ravnlk, J.
Strojniški vestnik - Journal of Mechanical Engineering 61(2015)4, 254-264
of the flow can be observed in the qDNS simulation. This is also confirmed by the experimental results of [8] at blockage ratios of 1/4 and 1/6. As the equations were solved in a direct manner, the small structures in the flow were limited only by a time-step and a mesh resolution. In Fig. 8a the streamlines are presented for the qDNS simulation showing a variety of different size flow structures. Fig. 7c presents the time series of the drag and lift coefficients for the URANS simulation. After a transient time, the smaller structures in the flow are dampened out by the imposed RANS model and the flow again becomes periodic. This is further shown in Fig. 8b regarding streamlines for the URANS simulation. Only large structures of the flow are observed. Having only large structures captured by the simulation, the simulated flow is qualitatively similar in behaviour to the laminar flow. Furthermore the dampening of the smaller structures by the increased viscosity of the URANS model seemed to lower the average value of the drag coefficient in comparison with the direct simulation result.
The flow at the Reynolds number 1000 showed highly nonlinear behaviour thus the frequency for the calculation of Strouhal number was obtained by FFT analysis of the lift coefficient signal (Fig. 7b). The drag coefficient was obtained by averaging the signal (Fig. 7a).
Table 2. Comparison between the Strouhal number and the average drag coefficient for the flow at the Reynolds number of 1000
	B	St	Cd
Davis exp.	1/4	0.179	2.37
Davis num.	1/4	0.144	1.88
Davis exp.	1/6	0.153	2.11
Davis num.	1/6	0.131	1.76
qDNS	1/8	0.157	1.94
URANS-SA	1/8	0.114	1.53
A comparison of the data from [8], found in the literature, was made between both the Strouhal number and the drag coefficient in Table 2. Although the ratio of the cylinder diameter to the channel width (B) was similar but not equal, no other data at this Reynolds number could be found for the geometry used in our simulations. The Strouhal number and the drag coefficient values in Table 2 given by the URANS simulation show values noticeably lower than in the reference data, most likely due to the imposed model dampening out smaller flow structures with the artificially increased eddy viscosity from the model.
a)
bj	X
Fig. 8. Streamlines at Re = 1000; a) in qDNS simulation and b) in URANS simulation
5 CONCLUSIONS
Whilst many other discretisation methods received a lot more attention in the past, the usage of the boundary element method for the discretisation of the governing equations of the fluid flow is not widely known. The uses of fundamental solutions and with that the appropriate descriptions of some of the governing physics make this method appealing for use in fluids and could with more refinement come close or even surpass modern discretisation methods. Furthermore, as no artificial diffusivity is present within the discretisation scheme, as a direct consequence of the usage of Green's functions, the use of the DNS type simulations would be possible without additional changes of the algorithms.
This paper derived at URANS type turbulence models for the velocity-vorticity forms of the Navier-Stokes equations. The governing equations are discretized by the boundary element method. Boundary values of vorticity are determined as a part of the algorithm. Vorticity transport and
BEM-Based Algorithm for URANS Simulations of Flow over a Square Cylinder
263
Strojniski vestnik - Journal of Mechanical Engineering 61(2015)4, 254-264
turbulence model equations are solved using a domain decomposition approach.
The algorithm is tested on the flow over the square obstacle in a channel at laminar and turbulent flow regimes. Velocity profiles, the estimated Strouhal number and drag coefficient for the laminar flow regime were compared with the benchmark results of other authors and proved the validity of the algorithm.
The turbulent flow regime was simulated at Reynolds number 1000 using quasi DNS and URANS approaches using Spalart-Allmaras turbulence model. Time traces of the drag and lift coefficient have been presented as well as average values compared with experimental data. Good agreement was observed.
6 REFERENCES
[1]	Skerget, L., Hribersek, M., Zunic, Z. (2003). Natural convection flows In complex cavities by BEM. International Journal of Numerical Methods for Heat & Fluid Flow, vol. 13, no. 6, p. 720-735, DOI:10.1108/09615530310498394.
[2]	Popov, V., Power, H. (1999). The DRM-MD integral equation method: an efficient approach for the numerical solution of domain dominant problems. International Journal for Numerical Methods in Engineering, vol. 44, no. 3, p. 327-353, DOI:10.1002/(SICI)1097-0207(19990130)44:3<327::AID-NME506>3.0.M;2-6.
[3]	Taigbenu, A.E. (1995). The green element method. International Journal for Numerical Methods in Engineering, vol. 38, no. 13, p. 2241-2263, DOI:10.1002/nme.1620381307.
[4]	Skerget, L., Rek, Z. (1995). Boundary-domain integral method using a velocity-vorticity formulation. Engineering Analysis with Boundary Elements, vol. 15, no. 4, p. 359-370, DOI:10.1016/0955-7997(95)00036-N.
[5]	Paige, C.C., Saunders, M.A. (1982). LSQR: An algorithm for sparse linear equations and sparse least squares. ACM Transactions on Mathematical Software (TOMS), vol. 8, no. 1, p. 43-71, DOI:10.1145/355984.355989.
[6]	Lupse, J., Skerget, L., Ravnik, J. (2014). Velocity-vorticity RANS turbulence modeling by boundary element method. Engineering Analysis with Boundary Elements, vol. 39, p. 4452, DOI:10.1016/j.enganabound.2013.11.004.
[7]	Ternik, P., Rudolf, R. (2014). Laminar forced convection heat transfer characteristics from a heated cylinder in water based nanofluids. International Journal of Simulation Modelling, vol. 13, no. 3, p. 312-322, DOI:10.2507/IJSIMM13(3)5.271.
[8]	Davis, F.W., Moore, F., Purtell, P. (1984). A numerical-experimental study of confined flow around rectangular cylinders. Physics of Fluids, vol. 27, no. 1, DOI:10.1063/1.864486.
[9]	Mukhopadhyay, A., Biswas, G., Sundararajan, T. (1992). Numerical investigation of confined wakes behind a square cylinder in a channel. International Journal for Numerical Methods in Fluids, vol. 14, no. 12, p. 1473-1484, DOI:10.1002/fld.1650141208.
[10]	Breuer, M., Bernsdorf, J., Zeiser, T., Durst, F. (2000). Accurate computations of the laminar flow past a square cylinder based
on two different methods: lattice-boltzmann and finite volume. International Journal of Heat and Fluid Flow, vol. 21, no. 2, p. 186-196, DOI:10.1016/S0142-727X(99)00081-8.
[11]	Lamura, A., Gompper, G. (2002). Numerical study of the flow around a cylinder using multiparticle collision dynamics. The European Physical Journal E: SoftMatter and Biological Physics, vol. 9, no. 5, p. 477-485, DOI:10.1140/epje/i2002-10107-0.
[12]	Galletti, B., Bruneau, C.H., Zannetti, L., lollo, A. (2004). Low-order modelling of laminar flow regimes past a confined square cylinder. Journal of Fluid Mechanics, vol. 503, p. 161170, DOI:10.1017/S0022112004007906.
[13]	Turki, S., Abbassi, H., Ben Nasrallah, S. (2003). Two-dimensional laminar fluid flow and heat transfer in a channel with a built-in heated square cylinder. International Journal of Thermal Sciences, vol. 42, no. 12, p. 1105-1113, DOI:10.1016/S1290-0729(03)00091-7.
[14]	Turki, S., Abbassi, H., Nasrallah, S.B. (2003). Effect of the blockage ratio on the flow in a channel with a built-in square cylinder. Computational Mechanics, vol. 33, no. 1, p. 22-29, DOI:10.1007/s00466-003-0496-2.
[15]	Kelkar, K.M., Patankar, S.V. (1992). Numerical prediction of vortex shedding behind square cylinders. International Journal of Numerical Methods for Fluid, vol. 14, no. 3, p. 327-341, DOI:10.1002/fld.1650140306.
[16]	Okajima, A. (1982). Strouhal numbers of rectangular cylinders. Journal of Fluid Mechanics, vol. 123, p. 379-398, DOI:10.1017/S0022112082003115.
[17]	Skerget, L., Samec, N. (2005). BEM for the two-dimensional plane compressible fluid dynamics. Engineering Analysis with Boundary Elements, vol. 29, no. 1, p. 41-57, DOI:10.1016/j. enganabound.2004.08.007.
[18]	Kim, J., Moin, P., Moser, R. (1987). Turbulence statistics in fully developed channel flow at low Reynolds number. Journal of Fluid Mechanics, vol. 177, no. 1, p. 133-166, DOI:10.1017/ S0022112087000892.
[19]	Wu, J.C., Thompson, J.F. (1973). Numerical solutions of time-dependent incompressible Navier-Stokes equations using an integro-differential formulation. Computers & Fluids, vol. 1, no.
2,	p. 197-215, DOI:10.1016/0045-7930(73)90018-2.
[20]	Spalart, P.R., Allmaras, S.R. (1994). A one-equation turbulence model for aerodynamic flows. La Recherche Aérospatiale, vol. 1, no. 1, p. 5-21, DOI:10.2514/6.1992-439.
[21]	Skerget, L., Ravnik, J. (2009). BEM simulation of compressible fluid flow in an enclosure induced by thermoacoustic waves. Engineering Analysis with Boundary Elements, vol. 33, no. 4, p. 561-571, DOI:10.1016/j.enganabound.2008.08.003.
[22]	Ravnik, J., Skerget, L., Hribersek, M. (2004). The wavelet transform for BEM computational fluid dynamics. Engineering Analysis with Boundary Elements, vol. 28, no. 11, p. 13031314, DOI:10.1016/j.enganabound.2004.05.002.
[23]	Orlanski, I. (1976). A simple boundary condition for unbounded hyperbolic flows. Journal of Computational Physics, vol. 21, no.
3,	p. 251-269, DOI:10.1016/0021-9991(76)90023-1.
264
Lupse, J. - Skerget, L. - Ravnlk, J.