Amano R.S., Sunden B. (Eds.) Computational Fluid Dynamics and Heat Transfer: Emerging Topics
Подождите немного. Документ загружается.
Sunden CH010.tex 10/9/2010 15: 22 Page 386
386 Computational Fluid Dynamics and Heat Transfer
(say, 7.0) to guarantee the incompressibility of the fluid. The speed of sound, c
s
,
can be formulated as:
c
2
s
=
γp
0
ρ
0
= αU
2
(7)
where, U is the characteristic or maximum fluid velocity. The choice of α is a
compromise: it should be an adequate value to avoid making p
0
(or c
s
) too large,
which will force the time advancement of the simulation to become prohibitively
slow,andalsotoavoidyieldingMachnumbersthatcanviolatetheincompressibility
conditionofthefluid;α takes100inthe casestudiestobe presented, whichensures
the density variation be less than 1.0%.
Velocity rate of change
The SPH formulation of viscous effects in fluid flow requires careful modeling to
satisfy its adaptivity and robustness. In this respect, Morris et al. [8] presented the
following expression, equation (8), to simulate viscous diffusion.
$
1
ρ
∇·µ∇
v
%
i
=
j=i
m
j
(µ
i
+ µ
j
)v
ij
ρ
i
ρ
j
1
|r
ij
|
∂w
ij
∂r
i
(8)
Thisexpressionhas proventobe highly accurateinparticular forthesimulations of
lowReynoldsnumberplanarshearflows[7, 8].Thevariableµ denotesthedynamic
viscosity of the fluid.
In what concerns equation (2), the pressure gradient terms in the Navier–
Stokes (N–S) equations can be formulated using the following symmetrized [6]
SPH version:
−
1
ρ
∇p
i
=−
j=i
m
j
p
i
ρ
2
i
+
p
j
ρ
2
j
∇
i
w
ij
(9)
Therefore, the SPH momentum equation yields
dv
i
dt
=−
j=i
m
j
⎡
⎢
⎢
⎢
⎣
p
j
ρ
2
j
+
p
i
ρ
2
i
∇
i
w
ij
−
1.0
ρ
i
ρ
j
µ
i
+ µ
j
∇
i
w
ij
· r
ij
r
2
ij
+ 0.01h
2
v
ij
⎤
⎥
⎥
⎥
⎦
+ F
i
(10)
where, F is the body force per unit mass (m/s
2
). The quantity 0.01h
2
is used to
prevent singularities for r
ij
≈0.
Position update
ToassuretheSPHparticlesmovewithavelocityconsistentwiththeaveragevelocity
of its neighboring particles it used the variant proposed by Morris [10]
dr
i
dt
= v
i
+ ε
j
m
j
v
ij
¯ρ
ij
w
ij
(11)
Here, ¯ρ
ij
=(ρ
i
+ρ
j
)/2.0 and the factor ε is chosen as ε =0.5.
Sunden CH010.tex 10/9/2010 15: 22 Page 387
SPH – a versatile multiphysics modeling tool 387
Rate of change of temperature
Clearyet al.[16]derivedaSPHheattransferequationforsolidifyingmetals,which
can be rewritten for a fluid of constant heat capacity as:
4.0m
j
ρ
i
ρ
j
λ
i
λ
j
λ
i
+ λ
j
r
ij
·∇
i
w
ij
r
2
ij
+ 0.01h
2
T
ij
−
2.0m
j
ρ
i
ρ
j
µ
i
µ
j
v
ij
· r
ij
2
µ
i
+ µ
j
·
r
ij
·∇
i
w
ij
r
2
ij
+ 0.01h
2
2
c
p
dT
i
dt
=
j=i
(12)
where,λ and c
p
are theheat conductivityand the specific heatat constant pressure,
respectively. Due to theincompressibility of thefluid, pressure workwas neglected
and the second term on the right-hand side denotes the contribution from viscous
heat dissipation.
Time advancement
The4thorderRunge–Kuttaintegrationalgorithmisusedtoperformtheintegrations
of equations (4), (10), (11), and (12).The time step t is adaptive and determined
by the Courant–Friedrichs–Lewy (CFL) condition [8] given by:
t = min
⎛
⎝
c
1
h
max
i
(
15.0U,|v
i
|
)
,c
2
?
h
max acceleration
,
c
3
h
max sound speed
⎞
⎠
(13)
The Courant “safety” factors c
1
, c
2
, and c
3
are given as c
1
=0.125, c
2
=0.25, and
c
3
=0.25, respectively. In the right-hand side of equation (13), a quantity 15.0U
is introduced in the first constraint so that the three constraints may have similar
influence upon the calculation.
Smoothing kernel
High accuracy three-splines quintic kernel is used
w
ij
= n
d
1.0
h
β
⎧
⎪
⎪
⎨
⎪
⎪
⎩
(
3.0 −s
)
5
− 6.0
(
2.0 −s
)
5
+ 15.0
(
1.0 −s
)
5
if 0 ≤ s < 1.0
(
3.0 −s
)
5
− 6.0
(
2.0 −s
)
5
if 1.0 ≤ s < 2.0
(3.0 −s)
5
if 2.0 ≤ s < 3.0
0ifs ≥ 3.0
(14)
where, s =|r
ij
|/h, β is the dimensionality of the problem, and n
d
a normalization
constant, whichis determined by the normalization property [15] ofthe smoothing
kernel
&
w(|r −r
|,h)dr
= 1.0 (15)
For the present implementation, β =2 and n
d
=7.0/478.0π.
Sunden CH010.tex 10/9/2010 15: 22 Page 388
388 Computational Fluid Dynamics and Heat Transfer
Boundary
d
A
d
B
v
A
v
B
B
A
f
x
y
Fluid channelWall region
∆p
∆p
∆p/2
∆p/2
Figure 10.1. Illustration of initial particle arrangement and treatment of solid wall
boundary.
The smoothing length h is given as:
h = 1.001 ×
(
x
)
2
+
(
y
)
2
(16)
where, x and y are the initial spacing in the x and y direction, respectively,
between two neighboring particles, and x =y. Thus, the interpolation of the
field quantity with respect to one fluid particle extends to a region where about 61
neighboring particles are included.
Solid wall boundary
The treatment of solid wall boundaries is illustrated in Figure 10.1. Taking into
account the influence region of the employedquintic smoothing kernel, four layers
of‘dumb’wallparticles arearranged outsideofthe fluidchannelandparallelto the
realwallboundarywiththefirstlayerp/2awayfromthe wallplane.These‘dumb’
particles have the same spatial separation p as fluid particles. The density of a
‘dumb’particleis initiallyset equaltothereference fluiddensity, ρ
0
, anditremains
unchanged. Boundary particles interact with the fluid particles by contributing
to their density variations, and by prescribing viscous forces on the nearby fluid
particles.
Toensure thenon-penetratingconditionofsolidwallboundary, thecalculations
aremade withboundaryparticlesexertingcentral forcesonfluid particles [17].For
a pair of boundary and fluid particle separated by a distance r, the force per unit
mass f (r) is assumed to have the 12-6 Lennard–Jones form [18]:
f (r) = D
r
0
r
12
−
r
0
r
6
r
r
2
(17)
f (r) is set to zero if r > r
0
so that the force is purely repulsive. r
0
is taken to be
0.45p. During the calculations, only the component of f(r) normal to the solid
wallisconsidered; D has theunitsof(m/s)
2
andis taken as10.0F ×Characteristic-
Geometric-Length. Varying the magnitude of D has no particular influence on the
calculation results.
Nonslip wall boundary condition is guaranteed by using the method proposed
by Morris et al. [8]. For each fluid particleA, the normal distance d
A
to the wall is
Sunden CH010.tex 10/9/2010 15: 22 Page 389
SPH – a versatile multiphysics modeling tool 389
y
d
u
x
(y, ∞)
U
0
x
B0 B0B 0
−d
0
Figure 10.2. Combined Poiseuille and Couette flow.
calculatedandforeachboundaryparticleB,d
B
isobtained.Thenanartificialveloc-
ity v
B
=−(d
B
/d
A
)v
A
is applied to the boundar y particle B with an assumption of
zero velocity condition on the boundary plane itself. This artificial velocity v
B
is
used to calculate the viscous force applied on the nearby fluid particles, but not
used to locate the boundary particle position. Boundary particles are motionless
or moving at specifiedvelocities.The relative velocitybetween fluidand boundary
particles is:
v
AB
= v
A
− v
B
= χv
A
(18)
where
χ = min
χ
max
,1.0 +
d
B
d
A
(19)
χ
max
is used to exclude the extremely large values when fluid particles get too
close to the wall. In Morris et al. [8] good results were obtained for χ
max
=1.5
when simulating low Reynolds number planar shear flows. In the present study,
χ
max
=7.0is specified, whichis thelargestvalueof χ for theinitial regularparticle
arrangement.
10.2.2 Benchmarking
The validity of the above-outlined SPH formulation is examined for two bench-
marks dealing with heat and fluid flow, namely: combined transient Poiseuille and
Couette flow, and steady heatconvection for combinedPoiseuilleand Couette flow
conditions; these benchmarks were selected because analytical solutions for them
are available.
Transient combined Poiseuille and Couette flow
Thistestcaseis thefluidflowbetweentwoparallelinfiniteplates, locatedaty =−d
and y =d, respectively, as shown in Figure 10.2. The system is initially at rest. At
time t =0, the upper plate starts to move at constant velocity U
0
along the positive
x-axis. Concurrently, a body force F is aligned parallel to the x-axis. The resultant
velocity profile is dependent on a non-dimensional quantity, designated as B
B =
d
2
νU
0
F (20)
Sunden CH010.tex 10/9/2010 15: 22 Page 390
390 Computational Fluid Dynamics and Heat Transfer
and its time-dependent series solution is:
u
x
(y,t) =
F
2ν
(d
2
− y
2
) +
U
0
2d
(y + d)
−
∞
n=0
(−1)
n
16Fd
2
νπ
3
(2n +1)
3
cos
πy
2d
(2n +1)
exp
−
(2n +1)
2
π
2
ν
4d
2
t
+
∞
n=1
(−1)
n
2U
0
nπ
sin
nπ
2d
(y + d)
exp
−
n
2
π
2
ν
4d
2
t
(21)
The first and third terms in the right-hand side of equation (21) refer to the contri-
bution by the applied body force F and the other two terms result from the motion
of the upper plate.When time t approaches infinite, the steady-state solution takes
the form:
u
x
(y,∞) =
B
2
U
0
1 −
y
2
d
2
+
1
2
U
0
1 +
y
d
(22)
The physical proper ties of the fluid are specified as, ν =1.0×10
−6
m
2
/s and
ρ =1000.0kg/m
3
.Thevalueofd istakenas0.0005mandU
0
5.0×10
−5
m/s.Fifty
SPH fluidparticlesareset tospan thechannel.The bodyforce F is adjustedto give
different values of B. The results are presented in Figures 10.3a and b.
The overall agreement between simulations and exact analytical solutions is
excellent: the largest deviation between SPH simulation and series solution is
observed at the position close to the motionless wall, and it is less than 2.0%.
Steady heat convection for combined Poiseuille and Couette flow conditions
The second test case is conducted on the steady convective heat transfer between
two parallel infinite plates, which are kept at constant temperatures: T
e
and T
w
(T
e
> T
w
),respectively.TheupperplateismovingataconstantvelocityU
0
yielding
the Couette flow condition, as shown in Figure 10.4. Concurrently, a body force F
may be applied too.
The exact solution for the steady-state temperature profile is given as:
T
∗
(
y
)
=
1
2
(1 +y
∗
) +
Br
8
(1 −y
∗2
) −
Br
6
B(y
∗
− y
∗3
) +
Br
12
B
2
(1 −y
∗4
) (23)
where, Br stands for the Brinkman number which is the product of the Prandtl
number Pr by the Eckert number Ec (Br =Pr.Ec), y
∗
and T
∗
are the normalized
position and temperature, respectively, as defined below:
Br =
µc
p
λ
U
2
0
c
p
(T
e
− T
w
)
(24)
Sunden CH010.tex 10/9/2010 15: 22 Page 391
SPH – a versatile multiphysics modeling tool 391
4.8
(a)
4.0
3.2
2.4
1.6
0.8
0.0
−1.0
−1.0
−0.5
0.0
2.0
5.0
B = 10.0
0
y/d
u/U
0
0.5 1.0
4.8
5.6
(b)
4.0
3.2
2.4
1.6
0.8
0.0
−1.0
−1.0
−0.5
0.0
2.0
5.0
B = 10.0
0
y/d
u/U
0
0.5 1.0
Figure 10.3. Results for combined Poiseuille and Couette flow at: (a) time=
0.212 seconds, and (b)time approaching infinite. (Symbols represent
the SPH predictions and the full lines the analytical solutions).
y
∗
=
y
d
(25)
T
∗
(y) =
T(y) − T
w
T
e
− T
w
(26)
The parameter B is the ratio between the contributions from the body force to the
fluidflowandthemotionofupperplate, asexpressedinequation(20).Thefirstterm
intheright-handsideofequation(23)isduetotheheatconductionbetweenthetwo
Sunden CH010.tex 10/9/2010 15: 22 Page 392
392 Computational Fluid Dynamics and Heat Transfer
y
d
u U
0
, T T
e
u 0, T T
w
u u
x
(y)
U
0
x
−d
0
Body force F
T T(y)
Figure 10.4. Steady heat convection for combined Poiseuille and Couette flow
condition.
parallel plates, whereas the last three terms describe the viscous heat dissipation,
which for moderate values of B and small values of Br is negligible; however, for
high viscosityfluids, as expected, may besignificant. In thefollowing calculations,
to test the ability of the SPH methodology in predicting heat dissipation, the fluid
is assumed to have a relatively large viscosity of 1.0 Pa·s. Moreover, thegeometric
dimension is increased to d =0.05m and the moving velocity of the upper plate is
U
0
=0.05m/s. Other parameters specified are: heat conductivity, λ =0.1W/m·K;
heat capacity, c
p
=4,200 J/kg·K; density, ρ =1,000.0kg/m
3;
and the temperature
of the lower plate, T
w
=0
◦
C. By varying the magnitude of F and the upper plate
temperature T
e
, different B and Br values are obtained.
As before, 50 SPH fluid particles are set to span the channel. The results for
this particular case are depicted in Figures 10.5a and b, and, as observed, the SPH
predictions for the non-dimensional temperature compare well with the analytical
solution; the maximum deviation is less than 1%.
The heat convection process and the heat dissipation characteristics are also
reproduced well by the SPH simulations.
10.3 Control of the Onset ofTurbulence in MHD Fluid Flow
Magnetohydrodynamic (MHD) flow control finds g rowing importance in many
industrial settings, in particular in the metallurgical industry. The electromagnetic
force imposed by a stationary magnetic field may suppress the flow instabilities
of electrically conducting fluids and creates favorable flow conditions in fluid flow
devices, such as tundishes, as described by Szekley and Ilegbusi [19]. The work
reported by Hartmann andLazarus[20] was thefirst oneto experimentallyidentify
thechangeindragandthesuppressionofturbulencewhenamagneticfieldisapplied
toturbulentliquid metalflows. Fromthen on, numerousstudies haveaddressed the
interactions between the applied transverse (wall-normal) magnetic field and the
fluidflow,asreportedinLeeandChoi[21],Krasnovetal.[22], JiandGardner[23].
A magnetic field applied in the transverse flow direction is effective in reducing
the turbulence fluctuations and suppressing the near-wall streamwise vorticity; the
mean velocity profile is modified, therefore yielding a drag increase. In contrast to
Sunden CH010.tex 10/9/2010 15: 22 Page 393
SPH – a versatile multiphysics modeling tool 393
(a)
1.2
1.6
0.8
T
*
T
*
0.0
0.4
−1.0 −0.5
0.0
2.0
5.0
Br = 0.001
10.0
y/d
0.5 1.0
(b)
6.0
8.0
10.0
4.0
0.0
2.0
−1.0
−0.5
0.0
0.4
Br = 0.0001
1.0
y/d
0.5 1.0
Figure 10.5. Results for temperature profiles when (a) B =0 and (b) B =10.0.
(Symbols represent the SPH predictions and the full lines the
analytical solutions).
the extensive research effort dealing with MHD flows with an applied transverse
magnetic field, little attention was given to the flows interacting with a streamwise
magnetic field. The application of a streamwise magnetic field, to a flow, however,
is able to restrain the velocity fluctuations in the transverse plane of that flow and
thetransitiontoturbulencecanbedistinctlydelayed,orevensuppressedintheflow.
The ability of a streamwise magnetic field of controlling the onset of turbulence
in an electrically conducting fluid flow was recently corroborated by Jiang et al.
[24]. To achieve the most effective flow control with the least amount of energy
consumption isthe goal ofan MHD flow controlstrategy, whichcan be realizedby
combining the two above-mentioned magnetic fields.
Sunden CH010.tex 10/9/2010 15: 22 Page 394
394 Computational Fluid Dynamics and Heat Transfer
10.3.1 MHD flow control
Motionofelectricallyconductingfluidsacrossamagneticfieldinducestheso-called
Lorentz force, which, when properly handled, may suppress the flow instability
and restrain the flow turbulence. A magnetic field applied in the transverse (wall-
normal) direction of a fluid flow creates a ponderomotive body force and changes
themeanvelocityprofile duetotheformationoftheHartmannlayers atelectrically
insulating boundaries normal to the magnetic field.A magnetic field applied to the
flow in thestreamwise direction restrainsthevelocity fluctuationsin the transverse
plane of the fluid flow and prevents the onset of turbulence. In comparison to
a streamwise magnetic field, a magnetic field applied in the transverse direction
of the flow acts directly on the main flow and may prove to be more flexible in
controlling the flow; however, it increases the flow resistance and needs additional
energy consumption. Therefore, a reasonable control strategy may be realized by
bringing an angle between the magnetic field and the flow direction within the
rangeofthetwoabove-mentionedmagneticsupplyconditions, i.e.,0
◦
(streamwise)
and 90
◦
(transverse). In this way, both velocity fluctuations in the transverse and
streamwisedirection arerestrainedconcurrentlyandamoreeffective controleffect
may be achieved with some advantage in what concerns external energy usage.
The transition to turbulence of electrically conducting fluid flow in a two-
dimensional (2-D) channel constrained by two parallel infinite plates, which are
assumedelectricallyinsulated, asdepicted inFigure 10.6,is chosento illustratethe
underlying concept. Fluid enters the channel with a uniform velocity U
0
, initially.
ForaReynoldsnumberexceedingacriticalvalue(≈2,000),theturbulence, after
ashortentrylength,setsin. Forthesakeofrestrainingthetransitiontoturbulence,a
uniformmagneticfieldBisappliedatanangleθ tothemainfluidflowdirection.The
presentstudy, byadjustingtheangleθ,aimstoachievethebestrestrainingeffect.To
comparethe resultsobtained fordifferentvaluesofθ, a bodyforceF
B
isprescribed
in the x direction to balance the produced ponderomotive electromagnetic force
so as to maintain the maximum flow velocity (≈U
0
) approximately equal for the
casesofinterest.Therestraining effectis quantitatively evaluatedbythe timewhen
the turbulence sets in, or by measuring the magnitude of the “turbulence-related”
quantities – the mean velocity ¯u and the two velocity cor relations in the x and y
direction, i.e., u
¯
u
and v
¯
v
at a specified time.
B
θ
F
B
U
0
u (y, t )
y
x
Figure 10.6. MHD flow control for 2-D channel flow.
Sunden CH010.tex 10/9/2010 15: 22 Page 395
SPH – a versatile multiphysics modeling tool 395
10.3.2 MHD modeling
The MHD flow is governed by a set of coupled partial differential equations that
express the conservation of mass, momentum and the interaction between the flow
and the applied magnetic field. In the present work, it is assumed that the hydro-
dynamic Reynolds number Re is much larger than the magnetic Reynolds number
Re
m
, namely
Re =
UL
ν
>> Re
m
=
U
m
L
m
η
(27)
Here, η =(σµ
0
)
−1
, where ν, σ, and µ
0
are the kinematic viscosity, electrical
conductivity, and magnetic permeability of free space (µ
0
=4π ×10
−7
H/m),
respectively; U and L are the characteristic velocity and length of the fluid flow,
respectively; the subscript m denotes the corresponding magnetic quantities. This
assumptionguarantees thefluctuationsofthemagneticfieldduetothefluid motion
arevery smallascomparedtothe appliedfieldB; therefore,the totalmagneticfield
can be taken as uniform and time-independent, and the induced electrical current j
is given by:
j = σ(∇φ + V × B) (28)
where,V isthe fluid velocity vector and φ is the induced electrical potential, which
is produced by the interaction between the applied magnetic field, and the flow
vorticity w (=∇ ×V).The presentstudy does not aimto investigate theturbulence
patternorstructureofMHDturbulentflows[21–23],ofinteresthereistoinvestigate
the time duration prior to the onset of turbulence of the flow. The vorticity itself,
w , is considered to be negligible, and the induced electrical potential is assumed to
be zero, which yields [25]
j = σ(V × B) (29)
The Lorentz force defined per unit of volume, F takes the form
F = j ×B = σ(V × B) × B (30)
Accordingly, the MHD equations for the above-mentioned transient 2-D (x −y)
planar channel flow can be written as follows:
Dρ
Dt
=−ρ
∂u
∂x
+
∂v
∂y
(31)
ρ
Du
Dt
=−
∂p
∂x
+ µ
∂
2
u
∂x
2
+
∂
2
v
∂y
2
− σ( B sinθ)
2
u +F
B
(32)
ρ
Dv
Dt
=−
∂p
∂y
+ µ
∂
2
u
∂x
2
+
∂
2
v
∂y
2
− σ( B cosθ)
2
v (33)