ICSV17, Cairo, Egypt, 1822 July 2010 1
AEROACOUSTICS AT LOW MACH NUMBER MERGING
FLOWS AT DUCT JUNCTION
Garret C.Y. Lam, S.K. Tang
Department of Building Service Engineering, The Hong Kong Polytechnic University, Hung
Hom, Kowloon, Hong Kong, P.R. China
email: garret.lam@polyu.edu.hk
Randolph C.K. Leung
Department of Mechanical Engineering, The Hong Kong Polytechnic University, Hung Hom,
Kowloon, Hong Kong, P.R. China
Complicated ductworks are common in modern airconditioned buildings, in which duct
junctions are often found. Air flows merge at these junctions and initiate shear layers at junc
tion edge due to the different velocities and directions of flows upstream of a junction. These
shear layers eventually roll up to form vortical structures, whose interactions produce addi
tional noise in the ductwork. Previous investigations mainly focused on the propagation of
noise at duct junctions. However, the aeroacoustical generation of noise, and its propagation,
of merging flow at duct junction are seldom studied. Thus, this paper aims to report a nu
merical investigation of the aeroacoustics of merging flow at a Tjunction and this is done by
solving NavierStokes equations with the conservation element and solution element (CE/SE)
scheme. This method has already been demonstrated as a successful tool in investigating
various aeroacoustics problems such as jet screeching noise, noise generated at cavities. The
scheme is built on a comparatively simple, yet noble, framework which saves computational
resources while it yields accurate results. The first part of this paper presents a study of wave
propagation at Tjunctions, which illustrates the capability of CE/SE method in capturing the
acoustics accurately. The second part shows an investigation of 90° impinging merging flows
at Tjunctions. Since the flows in ducts are usually turbulent in practical situation, MILES
approach with wall modelling is adopted in the present study to account for the effects of tur
bulence. From the instantaneous vorticity at Tjunction, the vortices are found shedding at
the edge downstream of Tjunction. Furthermore, identifying the regions where the acoustic
or aerodynamic effect dominates is always a task in onestep aeroacoustics simulation. A
simple method using two neighbouring time signals is also introduced to identity these re
gions and its effectiveness is reported in this paper.
17
th
International Congress on Sound and Vibration (ICSV17), Cairo, Egypt, 1822 July 2010
2
1. Introduction
In all modern airconditioned buildings, ventilation systems are installed to convey the treated
air to different zones in the buildings and extract the used air through complicated ductworks, thus
provide comfort for occupants. However, this system also serves as a noise source and the noise
generated propagate along these ductworks. Thus, this deteriorates the indoor acoustical environ
ment.
Duct junctions, at which a smaller duct joins the main air ducts, are commonly found in these
ductworks and therefore act as an important role in the sound transmission. In the past several dec
ades, numerous researches were carried out to investigate the noise transmission across these junc
tions. Miles [1] used conformal mapping to evaluate the noise transmission across duct junction at
low frequency. Bruggeman [2] also studied the propagation of low frequency sound at Tjunctions
with experiments. Dubos et al. [3] applied modal decomposition to investigate the problems. Tang
[4, 5, 6] further extended the researches to the higher frequency range.
On the other hand, the noise generation by flows at duct junction are seldom studied. The
generation is mainly due to the flow merging at duct junctions. Studies for the ductwork aeroacous
tics mainly focused on the cavity flow, backward facing step rather than the merging flow at duct
junction. Some other studies such as Krijger et al. [7] examined the steady flow behaviour of merg
ing flow, but not its acoustics. Furthermore, these studies usually only covered the lowReynolds
number regime.
The present study, thus, aims at investigating the mechanism of the noise generation at Tduct
in a more realistic manner by onestep aeroacoustics simulations. This was done by solving the
NavierStoke Equations with the conservation element and solution element (CE/SE) method [8] in
twodimensional manner. This method has already been demonstrated as a successful tool in inves
tigating various aeroacoustics problems such as jet screeching noise, noise generated at cavities.
Since the practical flow in ductworks are usually turbulent. Effects of turbulence are accounted for
by applying the MILES approach [9] while logarithmic wall law is applied in the near wall region.
In the following section, the NavierStokes equations will first be mentioned. The numerical
scheme adopted will then be introduced and followed by demonstrating the acoustical behaviour of
excitation in Tduct. Finally, the aeroacoustics of Tduct will be analyzed and discussed.
2. Governing Equations
The unsteady full NavierStokes equations for compressible flow in two dimensions can be
written in strong conservation form as
(
)
(
)
0=
∂
∂
+
∂
∂
+
∂
∂
y
gg
x
ff
t
u
vmmvmmm

(1)
with
[
]
( )
[ ]
( )[ ]
[ ]
[ ]
xyyxyyyxyvm
xxyxxxyxxvm
T
m
T
m
T
m
qvug
qvuf
vpepvuvvg
upeuvpuuf
evuu
++=
++=
++=
++=
=
ττττ
ττττ
ρρρ
ρρρ
ρρρ
,,,0
,,,0
,,,
,,,
,,,
2
2
, (2)
where
ρ
, p, u and v are the density, thermodynamic pressure, velocities in x and ydirections re
spectively and m is the index ranging from 1 to 4 in this paper. The quantities e,
µ
, q
x
, q
y
τ
xx,
τ
xy
and
τ
yy
are the total energy, the fluid viscosity, heat flux and stress components respectively as defined
by,
17
th
International Congress on Sound and Vibration (ICSV17), Cairo, Egypt, 1822 July 2010
3
(
)
21
22
vup
e
+
+
−
=
ρ
γ
, (3)
ref
ref
TsT
Ts
T
/
/1
)2/3(
+
+
=
µ
, (4)
( )
( )
y
T
MP
q
x
T
MP
q
r
y
r
x
∂
∂
−
=
∂
∂
−
=
22
1
,
1
γ
µ
γ
µ
, (5)
∂
∂
−
∂
∂
=
∂
∂
−
∂
∂
=
∂
∂
−
∂
∂
=
x
u
y
v
x
v
y
u
y
v
x
u
xyxyxx
2
Re
1
3
2
,
Re
1
,2
Re
1
3
2
τττ
, (6)
where γ is the specific heat ratio, M is Mach number, P
γ
is the Prandtl number, Re is the Reynolds
number and T is the temperature. These equations are nondimensionalised with the characteristic
flow velocity, density, temperature, viscosity and length, i.e. U
ref
,
ρ
ref
, T
ref
,
µ
ref
and L of the prob
lem. Except otherwise stated, all quantities in the forthcoming discussions are dimensionless.
3. CE/SE Method
CE/SE method [8] was originally developed to be a nondissipative, yet comparatively simple
framework. It minimizes the applications of artificial conditions to numerical schemes or adhoc
parameters, which may not be consistent with the flow physics. Its unique features includes (i) the
identical treatment of space and time dimension, (ii) enforcing the spacetime flux conservation
locally and globally, (iii) solving the shock problem without using Riemann solver, (iv) extending
to multidimensional problem easily, (v) simple implementation of boundary conditions. Brief de
scription of this method is provided below and detailed formulation can be referred to references.
3.1 Conservation of spacetime flux in CE/SE method
CE/SE scheme utilizes the conservation of spacetime flux as shown in Figure 1. The spatial
dimension x, y together with the time t forms a three dimensional Euclidean space, E
3
. The conser
vation of spacetime flux governs an arbitrary volume V and can be written as
0
)(
=⋅
∫
VS
dsH, (7)
where
(
)
mvmmvmm
uggff,,
−
−
=
H is the spacetime flux vector for m = 1, 2, 3, 4 and S(V) denotes
the surface around an arbitrary spacetime region V and ds = ndσ, with n and dσ being the unit nor
mal vector pointing outward and the area of a surface element S(V) respectively.
Figure 1. Arbitrary spacetime volume V Figure 2. Decomposition of two
in two dimensional case. dimensional domain
3.2 Formulation of CE/SE method
CE/SE scheme first divides the domain into nonoverlapping conservation element (CE) and
solution element (SE) according to the basic geometry of the mesh. Triangle is taken as the basic
Spacetime volume V
x
y
t
ds
n
A
C E
B
D
F
G
17
th
International Congress on Sound and Vibration (ICSV17), Cairo, Egypt, 1822 July 2010
4
shape here. Figure 2 shows the decomposition of the two dimensional computational domain. Tri
angle ACE is the base of the mesh with its geometrical centre G. Point B, D, and F, are the geomet
rical centres of the adjacent triangles. In Figure 3, they show the definition of CE and SE. Here we
assume that the point A, B, C, D, E, F and G are at nth time level while the superscript ‘′’ and ‘″’
denote the points at n+1/2 th and n+1 th time level respectively.
Figure 3a. Conservation element CE. Figure 3b. Solution element SE.
The CE is then formed by extruding the hexagon ABCDEF in the time dimension with height
Δ
t/2, which is the polyhedral ABCDEFGA
′
B
′
C
′
D
′
E
′
F
′
G
′
as shown in Figure 3a. This CCE, where
flux conservation is enforced, can be further divided into 3 basic CEs (BCE) ABGFA
′
B
′
G
′
F
′
,
BCDGB
′
C
′
D
′
G
′
and DEFGD
′
E
′
F
′
G
′
.
On the other hand, the SE of this element at n+1/2 th time level, denoted by SE(G
*
, n+1/2), is
formed by the 4 planes A
′
B
′
C
′
D
′
E
′
F
′
G
′
, BGB
″
G
″
, DGD
″
G
″
and FGF
″
G
″
. The solution point G
*
′,
where the flow variables of this SE is stored, is the geometrical centre of the hexagon ABCDEF at
n+1/2 th time level. Within this SE, the flow variables can be expressed in terms of values at G
*
′
by first order Taylors expansion,
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
2/1
*
*
*
*
**
2/1
**
*:,,
+
′′
′
′′′
+
′′
−+−+−+=
′
n
G
mtG
G
myG
G
mx
G
m
n
GGm
ttyyxxGtyx
φφφφφ
, (8)
where
φ
m
denotes the flow variables u
m
, f
m
and g
m
,
φ
mx
,
φ
my
,
φ
mt
are the derivatives of
φ
m
in x, y and t
respectively and
(
)
2/1
**
,,
+
′′
n
GG
tyx is the coordinates of G
*
′ in E
3
space. The viscous terms f
vm
and g
vm
are assumed to be constant within SE. Moreover, the flow variables f
m
, g
m
, f
vm
and g
vm
can be ex
pressed in terms of u
m
, u
mx
and u
my
. Therefore, the spacetime flux leaving a surface can be calcu
lated as the flow variables at the centroid of the surface multiplied by the area of the surface, which
is clearly a function of u
m
, u
mx
and u
my
at the solution point, e.g. flux leaving BGB
′
G
′
is
(
)
(
)
(
)
(
)
[
]
*
**
,,
′
′′
′′
=
G
my
G
mx
G
m
GBBG
uuuflux
ϕ
, (9)
In the time marching scheme of CE/SE method, the flow variables at nth time level are as
sumed to be known. Consider the CE in Figure 3a, the flux leaving the 3 BCEs ABGFA
′
B
′
G
′
F
′
,
BCDGB
′
C
′
D
′
G
′
and DEFGD
′
E
′
F
′
G
′
are the functions of the solution u
m
at SE(A
*
, n), SE(C
*
, n) and
SE(E
*
, n) respectively, which are known according to the assumptions. Moreover, the flux leaving
top surface A
′
B
′
C
′
D
′
E
′
F
′
G
′
is just the solution value
(
)
*′G
m
u multiplied by its area. Thus, the solu
tion u
m
at G*
′
at the n+1/2 th time level can be determined by summing the flux leaving the 3 BCEs
and then dividing it by the area of A
′
B
′
C
′
D
′
E
′
F
′
G
′
due to the flux conservation of CE.
After all the solution values at the n+1/2 th time level are calculated, all the spatial derivatives
are determined by simplified Courant Number Insensitive Scheme [11], which is a central differenc
x
y
t
D
B
F
G
A
’
B
’
C
’
D
’
E
’
F
’
G
’
B
”
D
”
F
”
2
t
Δ
2
t
Δ
G
D
E
B
C
F
A
’
B
’
C
’
D
’
E
’
F
’
G
’
A
2
t
Δ
17
th
International Congress on Sound and Vibration (ICSV17), Cairo, Egypt, 1822 July 2010
5
ing scheme with a controlling parameter governed by the local Courant number for each element.
This scheme minimizes the numerical dissipations caused in a mesh of nonuniform grid size.
4. Acoustics of Tduct
The acoustics of Tduct is important in building service applications as the noise propagates
in the ductworks and also serves a good testing case for the CE/SE scheme on its capability of cap
turing the acoustics. Tang et al. [6] demonstrated that the plane wave theory only predicts correctly
for very low frequency propagation in a Tduct as some higher modes might be excited at the junc
tion, whose intensities may be even stronger than that of plane wave mode.
The configuration of the simulation is shown as in Figure 4. The duct width W and speed of
sound are taken as the reference length and speed respectively. The duct length L/W is 20. The
origin is set as the left corner of the junction as shown. Buffer zone of length 10W is added to the
three ends of the Tduct. The mesh consists of 307214 triangular elements with the mesh size of
0.025. Simulation was carried with time step size of 0.001. Plane waves are excited in the vertical
duct beneath buffer zone in the absence of flow. The initial condition of the simulation is set as
ρ
=
ρ
ref
, u = v = 0, p = p
ref
= 1/
γ
, where
γ
is the specific heat ratio of air. The boundary condition in
CE/SE method adopts ghost cell approach and sliding wall boundary condition is applied [12, 13].
The outflow condition is set as Type 1 NRBC condition incorporated with back pressure p
ref
as
stated by Loh[12], i.e.
φ
g
=
φ
i
except p
g
= p
ref
,
φ
xg
=
φ
xi
,
φ
yg
=
φ
yi
, where
φ
denotes the flow variables,
superscript ‘g’ and ‘i’ denote the ghost cell and the neighbourhood interior cell respectively.
Figure 4. Configuration of Tduct.
Table 1. Plane wave amplitudes calculated by CE/SE, Tang and Lam [6].
Wave number
kW/
π
Vertical duct
– CE/SE
Vertical duct
– Tang
Main duct –
CE/SE
Main duct –
Tang
0.25 0.35 0.36 0.65 0.65
0.5 0.45 0.46 0.62 0.63
0.667 0.57 0.58 0.55 0.57
0.75 0.64 0.63 0.51 0.52
1.25 0.31 0.31 0.20 0.20
2 0.91 0.98 0.02 0.00
2.5 0.62 0.63 0.13 0.12
2.7 0.80 0.83 0.11 0.09
3.2 0.56 0.58 0.09 0.10
Table 1 compares the plane wave amplitude in different section of the Tduct calculated by
CE/SE scheme and the work of Tang and Lam [6]. The left most column is the normalized wave
number k times the duct width W divided by
π
and the shown amplitudes are normalized with that
of excitation. In this table, the results of CE/SE scheme matched previous results quite well. This
indicates that CE/SE correctly predicted the amount of energy that plane wave mode remained
L=20
Buffer Zone
Length = 10
O
x
y
L=20
L=20
W=1
17
th
International Congress on Sound and Vibration (ICSV17), Cairo, Egypt, 1822 July 2010
6
when it pass through a duct junction. This proves the capability of CE/SE schemes in calculating
the acoustic propagation in Tducts.
5. Aeroacoustics of Tduct
In this case, the configuration is the same as that adopted in Section 4 except different inlet
boundary conditions are implemented and no acoustical excitation is introduced in the vertical duct.
Essentially, it series to calculate the noise generation by the Tduct. Two air streams with maxi
mum Mach number Ma = 0.1 flow from the left inlet and the top inlet of the Tduct and merge at
the duct junction. The Reynolds number of the flow takes a practical value of , Re ≈ 2.3×10
5
. Fur
thermore, a fully developed turbulent profile is adopted for both inlet velocity profiles. The outlet
boundary condition is still the same as that implemented in the previous section, but the wall condi
tion is changed to noslip one.
In order to capture correctly the unsteady boundary layer development and the vortex shed
ding, a new mesh with 461828 triangular elements is applied and the mesh points are clustered near
the wall region to capture the steep velocity profile. The smallest mesh size is 0.001 near the wall,
which corresponds to ∼16.8 wall units, and the time step size is 2.5×10
5
for this case.
At the Reynolds number chosen, the simulated flow is turbulent flow, thus MILES introduced
by Boris et al. [9] is adopted to accounting for the turbulent effect. It relies on the internal dissipa
tion of a numerical scheme to mimic the dissipation of the unresolved turbulent scale, i.e. no sub
grid terms are used to which CE/SE scheme is ready to adopt [13]. Despite the implementation of
MILES, the turbulent dissipation near noslip wall region is still not fully accounted. Therefore, a
simple wall modelling is applied to calculate the dissipation by wall at the near wall region. It is
based on the classical logarithmic law profile for the velocity with the van Driest compressibility
correction. The velocity at the boundary cell is inverted by iterative Newton method to calculate
the friction velocity, which leads to the turbulent dissipation by wall.
5.1 Unsteady merging flows in Tduct
Figure 5 shows the instantaneous vorticity downstream of the junction at t = 146. In this fig
ure, a thin shear layer is extended from the upstream corner of the junction to the middle of the
horizontal duct section. The vortex, however, is not shed by this shear layer and this rather happens
in the wake region at the corner downstream to the junction at a frequency of 0.44. This implies
that the vortex shedding mechanism follows the stability of wake flow and is in agreement with
Juniper [14] in his analysis of absolute instability of a confined shear flow. The primary and secon
dary vortexes then collide the lower wall at 4 < x <5 and bounce back to the upper wall. The sec
ondary vortex is conveyed to the lower wall again in the downstream while the primary vortex
sticks to the upper wall travelling downstream.
Figure 5. Instantaneous vorticity at time = 146; Zoom: Downstream of duct junction.
5.2 Identification of acoustical and aerodynamic regions
In the onestep simulation, different types of flow disturbances at same frequency, due to
acoustics or vortex shedding are usually encountered at the same location because the acoustics is
generated by the underlying flow disturbance. Different types of disturbances cannot be separated
by only comparing the spectrum of FFT analysis. A simple method of using the phase difference
17
th
International Congress on Sound and Vibration (ICSV17), Cairo, Egypt, 1822 July 2010
7
between two adjacent time signals can serves this purpose effectively and it relies on the different
propagation speeds of different types of disturbance like the two microphone method [15]. Thus, it
is suitable to apply in the low Mach number flows.
Consider two points P and Q in the flow with
Δ
x denoting their separation and the direction of
a travelling disturbance is assumed to be parallel to PQ. Their phases of that disturbance after per
forming FFT is denoted by
θ
P
and
θ
Q
respectively and
θ
PQ
, the phase difference is defined as
θ
PQ
=
θ
P

θ
Q
. If
θ
PQ
is equal to the phase difference calculated by the speed of sound, (
θ
PQ
)
a
, that distur
bance can be acoustical disturbance. On the other hand, when
θ
PQ
is equal to that calculated by the
nominal flow convection speed (
θ
PQ
)
u
, that disturbance can be counted as aerodynamic disturbance.
Moreover, in the frequency range (0 < f < 2), (
θ
PQ
)
a
is very small compared to (
θ
PQ
)
u
due to Ma <<
1 as shown in Figure 6. This implies that in a region with mixed disturbances of acoustics and
aerodynamics, an increasing effect of acoustics will shift
θ
PQ
towards (
θ
PQ
)
a
.
In this merging flow with Ma = 0.1, all the time signals are extracted along the middle line of
the horizontal duct with
Δ
x = x
P
– x
Q
= 0.2, where x
P
and x
Q
are the xcoordinates of locations P and
Q respectively. Figure 6 and 7 show the calculated
θ
PQ
at different locations along the middle of
horizontal duct and only those peak frequencies whose magnitudes are no less than 6dB the maxi
mum amplitude are shown. It shows that in the upstream region, i.e. x < 0, all the travelling distur
bances are acoustics in nature based on the above criterion. In the region 0 ≤ x ≤ 1, some
θ
PQ
devi
ate from (
θ
PQ
)
a
but some still stick to it as shown in Figure 7. Types of disturbances change rapidly
in this small region of junction. In the further downstream region, i.e. x > 1, the disturbances are
found to be aerodynamics in nature and this is the same situation at x ≈ 21. This indicates that the
duct length of 20 is still not far enough for the farfield assumption to be valid.
From the above results, the introduced method is simple, yet effective to identify the type of
dominant effect in the flow region. It certainly suits the applications in low Mach number flows as
the large disparity of the characteristic speeds.
Figure 6.
θ
PQ
at x < 0. Figure 7.
θ
PQ
at x > 0.
Conclusions
In the present study, the aeroacoustics of merging flow at Tduct junction has been numeri
cally studied with the CE/SE scheme, which has also been demonstrated in the calculation of acous
tical excitation in Tduct as an appropriate tool for investigating internal acoustics. The simulation
of merging flows in Tduct was also presented and the vortex shedding mechanism is found to fol
low the stability of wake flow. Finally, a method of identifying different regions where the acoustic
or aerodynamic effect dominates was introduced and was shown its effectiveness.
17
th
International Congress on Sound and Vibration (ICSV17), Cairo, Egypt, 1822 July 2010
8
6. Acknowledgements
Support received from the Research Grants Council of the HKSAR Government through grant
PolyU 5278/06E and from the Hong Kong Polytechnic University through grant J–BB2C is grate
fully acknowledged.
REFERENCES
[1] J.W. Miles, “The diffraction of sound due to right angled joints in rectangular ducts,” Journal
of the Acoustical Society of America 19, 572579 (1947).
[2] J. C. Bruggeman, “The propagation of lowfrequency sound in a twodimensional duct system
with T joint and right angle bends: Theory and experiments,” Journal of the Acoustical Soci
ety of America 82, 10451051 (1987).
[3] V. Dubos, J. Kergomard, A. Khettabi, J.P. Balmont, D.H. Keefe, and C.J. Nederveen, “Theory
of sound propagation in a duct with a branched tube using modal decomposition,” Acta Acus
tica united with Acustica 85, 153169 (1999).
[4] S.K. Tang and F.Y.C. Li, “On low frequency sound transmission loss of double sidebranches:
A comparison between theory and experiment,” Journal of the Acoustical Society of America
113, 32153225 (2003).
[5] S.K. Tang, “Sound transmission characteristics of Teejunctions and the associated length
corrections,” Journal of the Acoustical Society of America 115, 218227 (2004).
[6] S.K. Tang and G.C.Y. Lam, “On sound propagation from a slanted side branch into an infi
nitely long rectangular duct,” Journal of the Acoustical Society of America 124, 19211929
(2008).
[7] J.K.B. Krijger, B. Hillen, H.W. Hoogstraten, and M.P.M.G. Van Der Raadt, “Steady two
dimensional merging flow from two channels into a single channel,” Applied Scientific Re
search 47, 233246 (1990).
[8] S. C. Chang, “The Method of SpaceTime Conservation Element and Solution Element – a
New Approach for Solving the NavierStokes and Euler Equations,” Journal of Computa
tional Physics 119, 295–324 (1995).
[9] J.P. Boris, F.F. Grinstein, E.S. Oran, and R.L. Kolbe, “New insights into large eddy simula
tion,” Fluid Dynamics Research 10, 199227 (1992).
[10] B.S. Venkatachari, G.C. Cheng, B.K. Soni, S.C. Chang, “ Validation and verification of Cou
rant number insensitive CE/SE method for transient viscous flow simulations,” Mathematics
and Computers in Simulation 78, 653670 (2008).
[11] J.C. Yen and D.A. Wagner, “Computational aeroacoustics using a simplified courant number
insensitive CE/SE method,” AIAA paper 20052820, 11
th
AIAA/CEAS Aeroacoustics Confer
ence, Monterey, California, USA, 2325 May 2005.
[12] C.Y. Loh, “ On a nonreflecting boundary condition for hyperbolic conservation laws,” AIAA
Paper 20033975, 16
th
AIAA Computational Fluid Dynamics Conference and the 33
rd
Fluid
Dynamics Conference and Exhibit, Orlando, Florida, USA, 2326 Jun 2003.
[13] C.Y. Loh, and L.S. Hultgren, “Jet screech noise computation,” AIAA Journal 44, 992997
(2006).
[14] M.P. Juniper, “The effect of confinement on the stability of twodimensional shear flows,”
Journal of Fluid Mechanics 565, 171195 (2006)
[15] J.Y. Chung and D.A. Blaser, “Transfer function method of measuring induct acoustic proper
ties. I. Theory,” Journal of the Acoustical Society of America 68, 907913 (1980).
Enter the password to open this PDF file:
File name:

File size:

Title:

Author:

Subject:

Keywords:

Creation Date:

Modification Date:

Creator:

PDF Producer:

PDF Version:

Page Count:

Preparing document for printing…
0%
Σχόλια 0
Συνδεθείτε για να κοινοποιήσετε σχόλιο