APP下载

Effects of finite-size heavy particles on the turbulent flows in a square duct*

2017-04-26ZhaowuLin林昭武XuemingShao邵雪明ZhaoshengYu余钊圣LianpingWang王连平

水动力学研究与进展 B辑 2017年2期
关键词:连平

Zhao-wu Lin (林昭武), Xue-ming Shao (邵雪明), Zhao-sheng Yu (余钊圣), Lian-ping Wang (王连平)

1.State Key Laboratory of Fluid Power and Mechatronic Systems, Department of Mechanics, Zhejiang University, Hangzhou 310027, China, E-mail: lin3030612023@163.com

2.Department of Mechanical Engineering, University of Delaware, Newark, Delaware 19716, USA

Effects of finite-size heavy particles on the turbulent flows in a square duct*

Zhao-wu Lin (林昭武)1, Xue-ming Shao (邵雪明)1, Zhao-sheng Yu (余钊圣)1, Lian-ping Wang (王连平)2

1.State Key Laboratory of Fluid Power and Mechatronic Systems, Department of Mechanics, Zhejiang University, Hangzhou 310027, China, E-mail: lin3030612023@163.com

2.Department of Mechanical Engineering, University of Delaware, Newark, Delaware 19716, USA

A parallel direct-forcing fictitious domain method is applied in fully-resolved numerical simulations of particle-laden turbulent flows in a square duct. The effects of finite-size heavy particles on the mean secondary flow, the mean streamwise velocity, the root-mean-square velocity fluctuation, and the particle concentration distribution are investigated at the friction Reynolds number of 150, the particle volume fraction of 2.36%, the particle diameter of 0.1 duct width, and the Shields number ranging from 1.0 to 0.2. Our results show that the particle sedimentation breaks the up-down symmetry of the mean secondary vortices, and results in a stronger secondary-flow circulation which transports the fluids downward in the bulk center region and upward along the side walls at a low Shields number. This circulation has a significant impact on the distribution of the mean streamwise velocity, whose maximum value occurs in the lower half duct, unlike in the plane channel case. The flow resistance is increased and the turbulence intensity is reduced, as the Shields number is decreased. The particles accumulate preferentially at the face center of the bottom wall, due to the effect of the mean secondary flow. It is observed that the collision model has an important effect on the results, but does not change the results qualitatively.

Turbulent duct flow, particle-laden flow, mean secondary flow, fictitious domain method

Introduction

The turbulent flow in a square duct is characterized by the presence of mean cross-stream fluid motions. This kind of secondary flows, induced by turbulence fluctuations, takes the form of eight symmetrical vortices, with two counter-rotating vortices in pairs in each quadrant of the duct (see Fig.2). The mean secondary flows transport the fluid momentum from the bulk region to the corner areas along each corner bisector, and then back to the bulk regions along the wall bisectors. The early experimental measurements of the turbulent flows in a square duct were focused on the Reynolds stresses as the source for the generation of mean secondary flows. Direct numerical simulations of the single-phase turbulent flow in a square duct were performed by Gavrilakis[1], Uhlmann et al.[2]and Pinelli et al.[3].

There are limited studies of the particle-laden turbulent flows in a square duct in the literature. Winkler et al.[4]investigated the preferential concentration of particles in a fully developed turbulent square duct flow, and observed that particles tended to accumulate in regions of high strain-rate and low swirling strength. Sharma and Phares[5]reported that the mean secondary flow enhanced the lateral mixing for passive tracers and low-inertia particles, and higher inertia particles accumulated close to the wall. Winkler et al.[6], Yao and Fairweather[7]and Yao et al.[8]investigated the particle deposition in turbulent square duct flows. The results of Winkler et al.[6]show that the deposition occurs with greater probability near the center of the duct walls than at the corners. On the other hand, Yao and Fairweather[7]and Yao et al.[8]concluded that high-inertia particles tend to deposit close to the corners of the duct floor, while low-inertia particles deposit near the floor center. Yao and Fairweather[9]investi-gated the resuspension of inertial particles in a turbulent square duct flow and demonstrated the important role of the mean secondary flow in the resuspension process.

In these studies of two-phase flows, the pointparticle approximation was employed to deal with the particle motion, which is valid in principle only when the particle size is smaller than the Kolmogorov length scale and the particle volume fraction is low. In recent years, the interface-resolved direct numerical simulation methods were used to study the mechanisms of the interactions between the turbulence and the finitesize particles, to determine the interface between a particle and the fluid and all turbulent structures with the direct numerical simulation method. Such methods were applied to simulations of particle-laden wallbounded turbulent flows such as the pipe flow[10]and the horizontal channel flows[11-15]. We investigated the effects of the finite-size neutrally buoyant particles on the turbulent flows in a square duct with the interfaceresolved DNS method, and observed that the particle addition increased the intensity of the mean secondary flow. In the present paper, we report our results on the particle effects on the turbulent duct flow when the particle sedimentation effect is present. It is noted that in the previous studies based on the point-particle approximation, the modulation of the turbulent duct flow by the particles has not been made clear.

Fig.1 Geometrical model for the duct flow

1. Numerical model

1.1Flow model

A schematic diagram of the geometrical model for the duct flow is shown in Fig.1. Thex-axis is aligned in the streamwise direction. Thez-axis direction is taken as the spanwise direction andy-axis direction as the transverse direction, as in the plane channel flow case. The corresponding velocity components in thedirections arerespectively. The no-slip velocity boundary condition is imposed at the duct walls and the periodic boundary condition is imposed in the streamwise direction. We denote the half width of the duct as. In the present study, the computational domain is

1.2Direct-forcing fictitious domain method

A parallel direct-forcing fictitious domain method (DF/FD) is employed for the simulation of particleladen turbulent duct flows. The fictitious domain (FD) method for the particulate flows was originally proposed by Glowinski et al.[16]. The key idea of this method is that the interior of the particles is filled with the fluids and the inner fictitious fluids are enforced to satisfy the rigid body motion constraint through a pseudo body force, introduced as a distributed Lagrange multiplier in the FD formulation[16]. In what follows, we describe the DF/FD method briefly, and the details can be found in Yu and Shao[17].

For simplicity of description, we will consider only one spherical particle in the following exposition. The particle density, volume and moment of inertia, the translational velocity, the angular velocity and the position are denoted byrespectively. Letrepresent the solid domain andthe entire domain including the interior and the exterior of the solid body. By introducing the following scales for the non-dimensionalization:for length,for velocity,for time,for the pressure, andfor the pseudo body force, the dimensionless FD formulation for the incompressible fluids and the spherical particles can be written as follows:

A fractional-step time scheme is used to decouple the system (1)-(5) into the following two sub-problems.

A finite-difference-based projection method on a homogeneous half-staggered grid is used for the solution of the above fluid sub-problem. All spatial derivatives are discretized with the second-order central difference scheme.

Note that the above equations are reformulated so that all the right-hand side terms are known quantities and consequently the particle velocitiesare obtained without iteration. Then,defined at the Lagrangian nodes are determined from

Finally, the fluid velocitiesat the Eulerian nodes are corrected from

In the above manipulations, the tri-linear function is used to transfer the fluid velocity from the Eulerian nodes to the Lagrangian nodes, and the pseudo body force from the Lagrangian nodes to the Eulerian nodes.

We note that we have extended the fictitious domain method to simulate the fluid-structure interactions[18,19].

1.3Collision model

The collision model is required to prevent the mutual penetration of particles and the penetration of particles into walls. Two soft-sphere collision models are adopted in the present study. One is the following artificial repulsive force collision model (referred to as ARF)

The other collision model adopted is the discrete element model (referred to as DEM) developed originally for the simulation of granular materials. In the DEM collision model, the mechanical elements such as a spring and a dash-dot are employed. Besides the spring-like repulsive force, the viscous damping force in the normal direction and the tangential friction forceare also considered. We adopt the model described in Crowe et al.[20]:

If the following relation is satisfied

The damping coefficients are given by

1.4Parameter settings

The friction Reynolds numberthe particle volume fraction is, the particlefluid density ratio isand the particle radiusnormalized byis, throughout this study. Three different Shields numbers are considered:and 1.0. The Shields number reflects the relative importance of the shear force and the buoyant force on a particle at the wall[21], and is defined as

The value of the Froude number in Eq.(8) can be determined from the Shields number as

The grid number in our computations is 1 024× 128×128, corresponding to the mesh sizeThe time step is. The flow statistics are obtained by averaging the data in the real fluid domain outside the particle boundaries over a period of typically fifty non-dimensional time units after the statistically stationary state is reached.

1.5Validation

Due to lack of benchmark data for the duct flow laden with finite-size particles, we compare our results of the mean and rms velocities in the single-phase case with the previous direct numerical simulations ofGavrilakis[1]and Pinelli et al.[3]to validate the accuracy of our computations, and it is shown that our results are in good agreement with the previous results.

2. Results and discussions

We will present and discuss our results for the mean secondary flow, the mean streamwise velocity, the RMS velocities, and the distribution of the particle concentration in the following sub-sections.

Fig.2 Mean secondary velocity vectors and contours of the corresponding stream-functionsfor the singlephase turbulent duct flow. For the contours, the increment is 0.01, and the dashed lines correspond to negative values and continuous lines to positive ones

2.1Mean secondary motion

The velocity vectors and the stream-function of the mean secondary flow in the single-phase case are shown in Fig.2, in order to make a comparison with the particle-laden flows given below. The flow takes the form of a pair of counter-rotating vortices in each quadrant, transporting the fluid momentum from the bulk region to the corner along the corner bisector. The mean secondary stream-function is determined from the mean streamwise vorticityby using the following equation

The maximum of the stream-function reflects the flow rate in the circulation and thereby the intensity of the mean secondary flow. The magnitude of the mean secondary velocity is about one percent of that of the bulk velocity, and consequently an accurate computation of the mean secondary flow is a challenging task.

Fig.3 Mean secondary velocity vectors for different Shields numbers () and collision models ( DEM and ARF)

Figures 3 and 4 show the velocity vectors and the stream-function of the mean secondary flow forθ= 1.0 and 0.2, respectively. As mentioned earlier, theShields number represents the ratio of the shear force to the buoyant force on a particle. The particle settling effect is stronger at a lower Shields number, and atmost particles settle down to the bottom wall, as shown in Fig.11. Due to the particle sedimentation, the original top-bottom symmetry of the secondary vortices is broken. From the bottom wall to the top wall, the first and the third vortices (i.e., the anticlockwise rotating vortices with a positive stream-function value in the left half duct) are weakened, and the second and the fourth are enhanced by the particles, as seen in Figs.2, 3 and 4. The second and fourth vortices are connected to each other and appear to merge in the particle-laden cases, resulting in a stronger and larger circulation which transports the fluid downward in the center region and upward along the side walls.

Fig.4 Contours of the mean stream-functions1.0 and 0.2 with the DEM collision model. The increment for the contours is 0.01

The suppression of the first vortex in the vicinity of the bottom wall is expected to be primarily caused by the hindering effect of the sediments. In Fig.3, the mean secondary flows obtained from ARF and DEM collision models are also compared, in order to examine the effects of the collision model. We see that the suppression of the first vortex is more significant for the DEM model, which is understandable due to the fact that the tangential friction between the particles and the wall in the DEM model hinders the motion of the particles. The effect of the collision model is important, but does not change the results qualitatively.

Fig.5 Contours of the normal stress difference gradient term in the mean vorticity equation and the normal stress differencefor the particle-free case and the particle-laden case at. The contour increment is 6 for the normal stress difference gradients and 0.1 for the normal stress difference

It seems that there exists a competition between the four vortices. The suppression of the first vortex is expected to be beneficial to the enhancement of the intensity of the second vortex. In addition, the particleeffect may directly enhance the second vortex. In our previous work, it was observed that the particle addition increases the intensity of the mean secondary flow in the neutrally buoyant case and we explained this result by using the mean streamwise vorticity equation. This equation for the fully developed singlephase flow takes a simple form[1]as

Fig.7 Mean streamwise velocity distributions along a line parallel to therespectively

Fig.8 The bulk velocities and mean turbulent kinetic energies as the function of Shields number

The terms in the first bracket on the left-hand side of Eq.(26) represent the convection of the mean vorticity by the secondary flow itself. The last term represents the viscous diffusion of the mean vorticity. The other two terms are related to the Reynolds stresses due to the Reynolds averaging of the Navier-Stokes equation, generally described as the source of vorticity[1]. The first term is associated with the gradients in the Reynolds cross-stream normal stress difference, and the second term is associated with theReynolds secondary shear stress. The flow statistics in the particle-laden case in our simulations are defined in the real fluid domain, thus Eq.(26) should hold approximately true for the particle-laden flows at the relatively low particle volume fractions, as indeed evidenced in our previous simulations in the neutrally buoyant case. It was demonstrated that the normal stress term is more important than the shear stress term for the generation of the mean secondary flow. Therefore, we here only plot the distribution of the normal stress differenceand its gradients in both particle-free and particle-laden cases for the comparison in Fig.5. In the particle-laden case at, we are only concerned with the flow statistics forsince the region forare mainly occupied by the particle sediments, in view of the fact that the particle diameter is.The particle-induced enhancement of the gradients of the normal stress difference in the regions where the second and fourth vortices are located can be observed in Fig.5(b), as compared to Fig.5(a). The enhancement is most pronounced near the side wall at. We then compare the normal stress difference in the particlefree and particle-laden cases in Figs.5(c) and 5(d), respectively. The enhancement of the normal stress difference near the side wall atby the particle effect is clear, which explains the enhancement of its gradients.

Fig.9 (Color online) Contours of the streamwise, transverse and spanwise rms velocity components for

2.2Mean streamwise velocity

We now discuss our results of the mean velocity. The contours of the mean streamwise velocity in the particle-free case and the particle-laden case at0.2 are compared in Fig.6. Although the intensity of the mean secondary flow is smaller by two orders of magnitude than that of the main-stream velocity (i.e., the bulk velocity), it does significantly affect the distribution of the mean streamwise velocity. In the neutrally buoyant case, the mean secondary flow tran-sports the high-speed fluids from the bulk region to the corner region, and the low-speed fluids from the wall region to the bulk region along the center-line, resulting in the distortion of the contours of the streamwise mean velocity in Fig.6(a), namely, with relatively higher streamwise mean velocity along the diagonal line than along the centerline of the cross section at the sameIn the particle-laden case at, the enhanced secondary circulation due to the second and fourth vortices discussed earlier makes the distribution of the streamwise mean velocity more inhomogeneous along the spanwise direction in the upper half duct, and the significant attenuation of the bottom vortex pair makes the distribution of the mean streamwise velocity more homogeneous spanwisely in the lower half duct. In addition, it is interesting to note that the maximum streamwise mean velocity occurs in the lower half duct as a result of the mean secondary flow effect (i.e., the mean secondary flow transports the fluids downwards in the center region), in contrast to the observation in the plane channel flow where there is no mean secondary flow and the maximum mean velocity occurs in the upper half channel[11].

Fig.10 The distributions of the streamwise, transverse and spanwise rms velocity components along a line parallel to thez-axis aty=−0.6, and a line parallel toy-axis atz=−0.5

The profiles of the streamwise mean velocity atfor different Shields numbers are shown in Fig.7, which is a complement to Fig.6 and directly shows that the maximum streamwise mean velocity occurs in the lower half duct due to the particle settling effect. In addition, Fig.7 indicates that theflow rate (or the bulk velocity) is decreased with the decrease of the Shields number, which means that the flow resistance is enhanced as the particle settling effect increases, in view of the fact that the pressure gradient is fixed and the lower flow rate corresponds to the higher flow resistance. Figure 8 shows the bulk velocity as a function of the Shields number for two collision models, and the results again indicate that the effect of the collision model is significant (the relative difference in the bulk velocity atis around 4%), but does not change the results qualitatively. The higher flow resistance for the DEM collision model can be explained by the fact that the friction between the particles and the wall is considered for the DEM model.

2.3rms velocities

We have discussed the role of the Reynolds normal stress difference in the generation of the mean secondary flow. The Reynolds normal stress component is actually the rms velocity. Now we inspect the effects of the particles on the rms velocity components. Figure 9 shows the contours of three rms velocity components forrespectively, with an overview of the rms velocity distribution. As the Shields number is decreased from 1.0 to 0.2, the streamwise rms velocity becomes smaller in most regions, whereas the transverse and spanwise rms velocities are, roughly speaking, enhanced in the lower half duct and weakened in the upper half duct. Figure 10 shows the distribution of the rms velocity components along a line parallel to theand a line parallel towhich further supports the above observation. One reason for the overall decrease of the streamwise rms velocity at a lower Shields number is that the bulk velocity decreases with the decrease of the Shields number. The further decrease of the streamwise mean velocity near the top wall (see Fig.7(a)) due to the effect of the mean secondary flow is responsible for the pronounced attenuation of the streamwise rms velocity near the top wall, and the effect of the particle sediments results in the attenuation of the streamwise rms velocity near the bottom wall (from the comparison between Figs.9(a) and 9(b)). The particle-induced vortices (or vortex shedding) give rise to the enhancement of the transverse and spanwise rms velocity components near the wall. Along the line of, the particle-induced enhancement of the wall-tangential rms velocityis larger than that in the wall-normal rms velocityin the near-wall region, which explains the enhancement of the magnitude of the normal stress differenceobserved earlier.

Shao et al.[11]investigated numerically the effects of the settling particles on the turbulent channel flow at a constant flow rate, and observed that the particle sediments increased significantly the rms velocities in the bulk region due to the vortex-shedding from the particles. In contrast, the present simulations show that the turbulence intensity is reduced when the particle settling effect is increased, as shown in Fig.8(b). The primary reason for the discrepancy may be that in the present study the pressure gradient is fixed and the flow rate decreases with the decrease of the Shields number, whereas the flow rate is fixed in Shao et al.[11].

2.4Particle distribution

It is shown that in the neutrally buoyant case the particle concentration is higher in the near-corner region. The distributions of the particle concentration (or the local volume fraction) forare shown in Fig.11. One can see that the particle concentration is higher at the face center of the bottom wall, which is clearly caused by the mean secondary flow. For, most particles settle down to the bottom wall and form one layer of particle sediments, and the particle can be trapped in the corner for a long time once it is bumped there due to the interaction with other particles, resulting in a high particle concentration there.

Fig.11 (Color online) The spatial concentration distributions of the particles in the cross section for0.2

3. Conclusions

The particle-laden turbulent flows in a squareduct are numerically simulated with a parallel directforcing fictitious domain method. The effects of finitesize heavy particles on the mean secondary flow, the mean streamwise velocity, the root-mean-square of the velocity fluctuation, and the particle concentration distribution are investigated at the friction Reynolds number of 150, the particle volume fraction of 2.36%, the particle diameter of 0.1 duct width, and the Shields number ranging from 1.0 to 0.2. From our results, the following conclusions can be drawn:

(1) The particle sedimentation breaks the updown symmetry of the mean secondary vortices, and results in a stronger circulation which transports the fluids downward in the bulk center region and upward along side walls at a low Shields number. This circulation has a significant impact on the distribution of the mean streamwise velocity, whose maximum value occurs in the lower half duct, unlike the plane channel case.

(2) The flow resistance is increased and the turbulence intensity is reduced, as the Shields number is decreased.

(3) As the Shields number is decreased, the streamwise rms velocity becomes smaller in most regions, whereas the transverse and spanwise rms velocities are, roughly speaking, enhanced in the lower half duct and weakened in the upper half duct. Along the line ofthe particle-induced enhancement of the wall-tangential rms velocityis larger than that in the wall-normal rms velocityin the near-wall region, resulting in the enhancement of the magnitude of the normal stress difference term there, which might be one reason for the enhancement of the second mean secondary vortex from the bottom wall.

(4) The particles accumulate preferentially at the face center of the bottom wall, due to the effect of the mean secondary flow.

(5) The collision model has an important quantitative effect on the results, but does not change the results qualitatively.

Acknowledgements

LPW is supported by the U.S. National Science Foundation (NSF) (Grant No. CBET-1235974) the Air Force Office of Scientific Research (Grant No. FA9550-13-1-0213), the Computing resources at Yellowstone supercomputer are provided by National Center for Atmospheric Research (Grant Nos. CISLP35751014, CISL-UDEL0001 and the University of Delaware through NSF (Grant No. CRI 0958512).

[1] Gavrilakis S. Numerical simulation of low-Reynolds-number turbulent flow through a straight square duct [J].Journal of Fluid Mechanics, 1992, 244: 101-129.

[2] Uhlmann M., Pinelli A., Kawahara G. et al. Marginally turbulent flow in a square duct [J].Journal of Fluid Me-chanics, 2007, 588: 153-162.

[3] Pinelli A., Uhlmann M., Sekimoto A. et al. Reynolds number dependence of mean flow structure in square duct turbulence [J].Journal of Fluid Mechanics, 2010, 644: 107-122.

[4] Winkler C. M., Rani S. L., Vanka S. P. Preferential concentration of particles in a fully developed turbulent square duct flow [J].International Journal of Multiphase Flow, 2004, 30(1): 27-50.

[5] Sharma G., Phares D. J. Turbulent transport of particles in a straight square duct [J].International journal of multiphase flow, 2006, 32(7): 823-837.

[6] Winkler C. M., Rani S. L., Vanka S. P. A numerical study of particle wall-deposition in a turbulent square duct flow [J].Powder Technology, 2006, 170(1): 12-25.

[7] Yao J., Fairweather M. Particle deposition in turbulent duct flows [J].Chemical Engineering Science, 2012, 84: 781-800.

[8] Yao J., Fairweather M., Zhao Y. L. Numerical simulation of particle deposition in turbulent duct flows [J].Industrial and Engineering Chemistry Research, 2014, 53(8): 3329-3341.

[9] Yao J., Fairweather M. Inertial particle resuspension in a turbulent, square duct flow [J].Physics of Fluids, 2010, 22(3): 033303.

[10] Wu T. H., Shao X. M., Yu Z. S. Fully resolved numerical simulation of turbulent pipe flows laden with large neutrally-buoyant particles [J].Journal of Hydrodynamics, 2011, 23(1): 21-25.

[11] Shao X., Wu T., Yu Z. Fully resolved numerical simulation of particle-laden turbulent flow in a horizontal channel at a low Reynolds number [J].Journal of Fluid Mechanics, 2012, 693: 319-344.

[12] Kidanemariam A. G., Chan-Braun C., Doychev T. et al. Direct numerical simulation of horizontal open channel flow with finite-size, heavy particles at low solid volume fraction [J].New Journal of Physics, 2013, 15(2): 025031.

[13] Do-Quang M., Amberg G., Brethouwer G. et al. Simulation of finite-size fibers in turbulent channel flows [J].Physical Review E, 2014, 89(1): 013006.

[14] Picano F., Breugem W. P., Brandt L. Turbulent channel flow of dense suspensions of neutrally buoyant spheres [J].Journal of Fluid Mechanics, 2015, 764: 463-487.

[15] Wang L. P., Peng C., Guo Z. et al. Flow modulation by finite-size neutrally buoyant particles in a turbulent channel flow [J].Journal of Fluids Engineering, 2016,138(4): 041306.

[16] Glowinski R., Pan T. W., Hesla T. I. et al. A distributed Lagrange multiplier/fictitious domain method for particulate flows [J].International Journal of Multiphase Flow, 1999, 25(5): 755-794.

[17] Yu Z., Shao X. A direct-forcing fictitious domain method for particulate flows [J].Journal of Computational Physics, 2007,227(1): 292-314.

[18] Shao X. M., Zhang X. L., Yu Z. S. Numerical studies of the hysteresis in locomotion of a passively pitching foil [J].Journal of Hydrodynamics, 2016, 28(3): 359-368.

[19] Zhang H. Z., Yu Z. S., Shao X. M. Numerical study on the propulsion of a biomimetic jellyfish [J].Chinese Journal of Hydrodynamics, 2016, 31(3): 327-333(in Chinese).

[20] Crowe C. T., Schwarzkopf J. D., Sommerfeld M. et al. Multiphase flows with droplets and particles [M]. Boca Raton, USA: CRC Press, 2011, 124-132.

[21] Peysson Y., Ouriemi M., Medale M. et al. Threshold for sediment erosion in pipe flow [J].International Journal of Multiphase Flow, 2009, 35(6): 597-600.

(Received August 12, 2015, Revised January 27, 2016)

* Project supported by the National Natural Science Foundation of China (Grant Nos. 11372275, 51376162), the Research Fund for the Doctoral Program of Higher Education of China (Grant No. 20130101110035).

Biography: Zhao-wu Lin (1984-), Male, Ph. D. Candidate

Zhao-sheng Yu,

E-mail: yuzhaosheng@zju.edu.cn

猜你喜欢

连平
连平:房地产政策宽松仍需发力
美联储加息与人民币汇率之间的相关性
连平:房地产政策有必要进一步释放维稳信号
连平:对中国房地产市场中长期的发展不应悲观
叶连平 乡村永不熄灭的烛光
叶连平
九旬乡村教师叶连平:“给孩子们上课,我就是享福”
面孔
连平举办精品奇石展
广东省龙川至怀集公路龙川至连平段模块化机房建设