A Comparison Between CFD And Network Models For Predicting Wind-driven .

1y ago
543.50 KB
7 Pages
Last View : 3d ago
Last Download : 4m ago
Upload by : Xander Jaffe

ARTICLE IN PRESSBuilding and Environment 42 (2007) 4079–4085www.elsevier.com/locate/buildenvA comparison between CFD and Network models for predictingwind-driven ventilation in buildingsOmar S. Asfour , Mohamed B. GadiSchool of the Built Environment, The University of Nottingham, University Park, Nottingham NG7 2RD, UKReceived 23 August 2006; received in revised form 12 October 2006; accepted 1 November 2006AbstractThe aim of this paper is to compare the implementation of computational fluid dynamics (CFD) and Network models for airflow rateestimation in buildings. The CFD software used is Fluent 5.5. Comparison between the predicted and simulated airflow rate is suggested as avalidation method of the implemented CFD code, while the common practice is to compare CFD outputs to wind tunnel or full-scalemeasurements. This could be useful for studies that have no access to laboratory or full-scale testing facilities. Results obtained from testing anumber of cases have been compared and analysed, considering normal and oblique wind directions. The comparison held between mathematicaland CFD results generally showed a good agreement, which seems to justify the use of CFD code for predicting natural ventilation in buildings.r 2007 Elsevier Ltd. All rights reserved.Keywords: Natural ventilation; CFD; Network model; Validation1. IntroductionApplication of computational fluid dynamics (CFD) wasoriginally introduced for industrial applications, but todayit has also become a common tool for assessing buildingventilation and environmental performance. As describedby CIBSE [1], CFD is ‘‘a very powerful technique’’ inpredicting air movement and characteristics. CFD model isbased on the concept of dividing the solution domain intosub-zones. Then, for each zone, the mass, momentum, andenergy conservation equations are solved, utilising theprocessing power of computers. This helps to performcalculations more easily and, in comparison with naturalventilation mathematical models, gives more detailedresults. For example, CFD codes are used to predictairflow rate, air velocity and temperature, and airflowpatterns inside and around buildings. Many softwaresbased on CFD codes have been developed like, Fluent,Phoenics, and Flovent.In the last few years, an intensive work has been doneusing CFD. Comparisons of CFD results with wind tunneltests have shown good agreement. However, some studiesare limited to the use of CFD modelling, and have noaccess to experimental testing facilities. In this case, it iscrucial to validate the implemented CFD code prior to theproposed modelling study in order to avoid producing anymisleading results. This is because CFD is considered to bea sophisticated modelling technique. One possible option isto compare CFD results to natural ventilation mathematical models, which is the scope of this paper. These modelsvary from simple to complicated models. The commonlyused Network mathematical model has been found to beappropriate for this study, where calculation process isreasonable.Both implemented CFD and Network models are brieflyexplained below. This followed by the results obtainedfrom both models regarding airflow rate estimation in thecase of normal and oblique wind direction.2. Network mathematical model Corresponding author. Tel.: 44 115 9513118; fax: 44 115 9513159.E-mail addresses: laxosa@nottingham.ac.uk, o.asfour@hotmail.com(O.S. Asfour), mohamed.gadi@nottingham.ac.uk (M.B. Gadi).A flow network is used here to represent different airflowpaths in buildings. This network is represented by a group0360-1323/ - see front matter r 2007 Elsevier Ltd. All rights araa.com

ARTICLE IN PRESSO.S. Asfour, M.B. Gadi / Building and Environment 42 (2007) 4079–40854080of nodes, simulating building openings, and a group oflines, simulating flow paths. Thus, it has more flexibility inanalysing natural ventilation problems. Bahadori andHaghighat [2] explained the implementation of this model.Assuming a multi-zone problem, i.e., a building with Nzones, a number of N equations will be established in theNetwork model. Each equation is dependant on theprevious one.Thus, these equations are solved iteratively, where thesummation of mass flow rates should be zero. In the case ofwind-induced ventilation, the knowledge of internal andexternal pressure coefficients for each opening is requiredto find out pressure difference across this opening. Thesecoefficients are usually obtained experimentally or fromstandard pressure coefficients data. The following equationis commonly used in predicting airflow rate induced bywind pressure �ffiffiffiffiffi Qn ¼ Aeff(1)2Dp r ,where Qn is airflow rate through an opening n (m3/s), Aeff isthe effective area of this opening, DP is pressure differenceacross it (Pa), and r is air density (kg/m3). Pressuredifference mentioned in the previous equation can beestimated using DP ¼ 0:5rV 2 C pn C pi ,(2)where V is wind velocity at datum level (m/s), Cpn ispressure coefficient at opening n, and Cpi is pressurecoefficient inside the space. Substituting Eq. (1) into Eq.(2), we get 1 2Qn ¼ Aeff V C pn C pi C pn C pi .(3)Considering the Law of Mass Conservation, this equationcan be rewritten asNX 1 2Aeff V C pn C pi C pn C pi ¼ 0:(4)n¼1Therefore, it is in possible to estimate internal pressurecoefficient using Eq. (4) and then airflow rate using Eq. (3).To do so, the knowledge of Cpn is essential, as discussed inSection 4.3. CFD codeFluent 5.5 program, one of the most widely usedcommercial CFD codes, has been used in this study.Gambit 1.3 program, which is a pre-processor program,has been used to define the building geometry. In fact,using two-dimensional modelling in room ventilationproblems does not give realistic simulation of airflow, asit does not consider some phenomena that determineairflow characteristics, such as airflow separation overbuilding sharp edges. Thus, the use of the three-dimensional modeling has been chosen in this study. Thedisadvantage, however, is that this method significantlyincreases the required time for solution convergence.Boundary and continuum types have been also defined inGambit program, where velocity inlet and outflowboundary types have been used for the solutiondomain. In order to predict airflow behaviour in themodelled buildings, CFD calculation mesh has beengenerated, which replaces the air inside and around thebuilding. The idea of this mesh is to divide the solutiondomain into small cells, which are used to predict airflowbehaviour using computer processing ability. This can beachieved using different types of mesh, like hexagonal ortetrahedral.For these cells, or volumes in the case of threedimensional modelling, Fluent 5.5 software can numerically solve the three basic conservation equations of mass,momentum and energy in an iterative manner. However,wind-induced ventilation is believed to be more effective inhot climates, when compared with stack-induced one. Thisis because of the relatively lower difference between indoorand outdoor temperatures, which is the main factoraffecting stack ventilation [3]. Therefore, energy settingshave been set off in this study.To process these calculations, an appropriate mathematical model has to be applied on the solution domain. Mostof airflow problems in buildings consider airflow to beturbulent. Thus, definition of turbulent model is requiredto help solving the transport equations. Awbi [4] mentioned that the standard k–e Model is believed to be themost used and developed turbulence model. This model ismost likely to predict reasonable results for airflow studiesin buildings. Turbulence characteristics have been specifiedusing the Turbulence Intensity and Length Scale option,which is recommended in room airflow problems, anddepends on the Reynolds number and the inlet size of thecase. Segregated solver has been used, where the fundamental equations are solved sequentially or segregatedfrom each other. The solution reaches the end when theconvergence criteria are met.The modelled prototype represents a room with an inletand an outlet. In Fluent 5.5 program, it is required todefine air velocity magnitude at this inlet. Air velocitymagnitude varies with height. Thus, it is crucial toestimate air velocity at building height, as will bediscussed in Section 5. Rotating the building modelinside the three-dimensional solution domain, whichrepresents the ambient air, can be used to simulate windangle. It is important to note that Guage pressure shouldbe set to 0. This is because the only pressure acts is theatmospheric pressure. The solution progresses in the formof several iterations, which can be monitored by the user.During this iteration process, the residual information ofthe velocity, continuity, and the turbulence parameters ofthe viscous model are continuously updated. Convergenceis achieved when a sufficient error tolerance, defined by theuser, is reached. Many outputs in different presentationmethods can be obtained from Fluent 5.5 software. Thisincludes airflow rates and contours of air velocitymagnitudes.www.manaraa.com

ARTICLE IN PRESSO.S. Asfour, M.B. Gadi / Building and Environment 42 (2007) 4079–408540814. Prediction of airflow rate using Network mathematicalmodelIt is possible to implement this model for any openingknowing its area and wind static pressure coefficient at thisopening. Prediction process here is divided into four stages: specification of pressure coefficient data;specification of building configurations;calculations of internal pressure coefficient;calculations of airflow rate.Pressure coefficients data that will be used here are thoserecommended by Liddament [5]. These data have beenproduced in wind tunnels for some common buildingconfigurations, and are considered more developed thanthe data mentioned in the British Standards [6]. This is interms of testing more wind angles and considering theeffect of terrain nature. Cases proposed here have nearlythe same volume, but different aspect ratios. Aspect ratiosof 1:1, 1:2 and 2:1 are considered, as depicted in Fig. 1.This figure also shows values of pressure coefficient, with areference number indicating wall number. Initially, calculations have been performed considering a relatively lowand high reference wind velocities, namely 1 and 5 m/s, andtwo wind angles, namely 01 and 451. However, it has beenobserved that discrepancy percentage is not sensitive towind velocity magnitude, as the percentage is nearly thesame in the case of same building geometry and differentapproaching wind speeds. This result supports whatLiddament [5] has mentioned that pressure coefficient isnormally assumed to be independent of wind speed but notdirection. Thus, the following investigation will be confinedto velocity magnitude of 1 m/s.It is important to note that at oblique wind direction,effective area of the windows is less. By knowing that windangle is 451, window effective width can be simplycalculated using Right Triangle Trigonometry. Thus,window effective area can be estimated by multiplyingwindow width by its height. This is 1.4 2 ¼ 2.8 m2.Internal pressure coefficients can be estimated from Eq.(4). As the above-illustrated models have only two zones:indoor and outdoor, this equation can be simplified as C pn C pi 1 2 þ C pðnþ1Þ C pi 1 2 ¼ 0.(5)It is required to know air velocity at building height toimplement Network model for airflow rate estimation. Thisis possible using the following common equation [7]:aV ¼ V r cH ,(6)where V is wind speed at datum level (m/s), Vr is referencewind speed (obtained from meteorological data) (m/s), H isthe height of the building, c is parameter relating windspeed to terrain nature (0.68 in the open country terrain),and a is an exponent relating wind speed to the heightabove the ground (0.17 in the open country terrain).Knowing different building heights, as illustrated inFig. 1. Illustration of the cases tested in this study.Table 1Airflow rate prediction, using the Network modelCaseOpeningAeff .250.125 0.1 0.025 0.2 0.275Qn (m3/s)Qn (kg/s)2.40 2.402.37 2.372.66 2.661.53 1.531.52 1.521.75 1.752.94 2.942.90 2.903.26 3.261.88 1.881.87 1.872.14 2.14Table 1, wind velocity at building height has beenestimated and implemented in Eq. (3) to estimate airflowrate. Results are shown in Table 1.www.manaraa.com

ARTICLE IN PRESS4082O.S. Asfour, M.B. Gadi / Building and Environment 42 (2007) 4079–4085This airflow rate is given in m3/s. It is possible to convertit to mass flow rate, in kg/s, by multiplying it by air densityin order to facilitate the comparison with CFD results.Value of air density should match the one under whichpressure coefficient data were estimated. As there is noindication of this value in the related resource, this valuewill be assumed to be 1.225 kg/m3, which is the defaultvalue used in Fluent 5.5 program and nearly the same ofthe standard air density value, i.e., 1.2 kg/m3.5. Prediction of airflow rate using CFD modellingUsing Fluent 5.5 program, it is possible to estimateairflow rate through building openings. Modelling processincludes the following four stages:Fig. 3. Illustration of the calculated and simulated wind velocity profile(for case 1, as an example). Table 2Air velocities for the different sub-inlets used in CFD modelling, asestimated from Eq. (6) drawing of building models in Gambit program;generating and exporting the calculation mesh to Fluent5.5 program;definition of solution code;computing airflow rate.As discussed in Section 3, building models are simulatedthree-dimensionally in order to obtain more realisticairflow pattern. This explains the significant differencesbetween results obtained using both two and threedimensional modelling in the early stage of this study.Thus, it is important to simulate the ambient air around thebuilding. To do so, the modelled room will be placed insidea three-dimensional box, as shown in Fig. 2.By rotating the room inside this box, it is possible tosimulate oblique wind direction. However, this requires alarger domain to ensure solution convergence. Forexample, in the case of 01 wind direction, the building isplaced in a 30 m 30 m 20 m box. In 451 wind direction,the extra sheer stresses, resulted by the oblique buildingwalls, cause reversed flows to occur. This leads the solutionto diverge, and makes it essential to use a larger domainsize. This size has been gradually increased until anacceptable size of 50 m 50 m 20 m was achieved (Fig. 3).In ventilation modelling, it is important to consider windspeed variation with height, due to the frictional effect ofthe ground. It is possible to do so using Eq. (6). Windvelocity profile can be defined along the velocity-inlet, asillustrated in Fig. 2, using ‘User Defined Function’ optionin Fluent 5.5 software. However, application of this optionFig. 2. Simulation of the ambient air around buildings.Height above the ground (m)Sub-inlet numberWind velocity (m/s)Cases 1 and 4 (room height is 5 m)511021532040.8941.0061.0781.132Cases 2,3,5 and 6 (room height is 4 m)41821231642050.8610.9681.0371.0891.132in three-dimensional simulation required some advancedC programming. It is possible, as an approximationmethod, to divide the velocity inlet into many sub-inlets.Each sub-inlet will have a different air velocity magnitudedepending on its height above the ground. This method hasbeen found to be useful and good results have beenobtained, as has been concluded from this study. To allowfor the estimation of airflow rate, it is recommended to splitthe solution domain into the following zones: room walls,ambient air, room interior, and openings volumes. The lastthree volumes were defined as fluid continuum. The benefitof this arrangement is that Fluent 5.5 program will defineopenings surfaces separately, which facilitates airflowcomputing process (Table 2).The use of three-dimensional simulation leads to thinkabout the resulting mesh sizes and the required processorcapacity. For example, the use of 0.2 m mesh spacing forthe entire domain, 30 cm 30 m 20 m in cases 1–3, wasbeyond computer memory and speed in this study. Acommon solution here is to create a hierarchy in mesh sizeto be fine inside the building and larger around it. In anycase, a trial-and-error process is recommended to find outthe most appropriate mesh configuration. Rough meshescan be suitable and sufficient in many CFD simulationwww.manaraa.com

ARTICLE IN PRESSO.S. Asfour, M.B. Gadi / Building and Environment 42 (2007) 4079–4085cases [8]. Generally, the following mesh configuration hasbeen found to be acceptable: Meshing scheme: hex-map or hex-submap. Mesh elements here are hexahedral, which is applicable to theshape and topological characteristics of the modelledcases. In case of oblique wind, the use of tetrahedralmesh has been found to give more accurate results, asdiscussed in the following section.Mesh node spacing: in case of normal wind direction,this spacing is 0.6 m. In the case of oblique winddirection, solution domain was divided into two zones:room interior, with 0.5 m spacing in the hex-map mesh,and room exterior, with 1 m spacing in the tetrahedralmesh.Once the mesh is ready, it can be exported to Fluent 5.5program in order to perform the calculations, which maytake several hours depending on the size and complexity ofeach individual case. Fig. 4 illustrates the resulting airflowpattern for these cases, presented by contours of velocitymagnitudes over a 2-m height section. Generally, the4083observed airflow pattern, in all the cases, has been found tobe reasonable for such sharp-edge geometries. Forexample, when wind reaches the windward face, a highpressure zone is formed there. This pressure pushes airinside, around, and over the building. Some standardfeatures can also be observed. This includes airflowseparation over building sharp edges. This phenomenonusually occurs when airflow layers hit a sharp edge of thebuilding and thus lose their momentum. After somedistance, the separated airflow joins its original streamagain in a point called the reattachment point.After achieving solution convergence, it is possible toobtain the mass flow rate directly from the software,utilising the Surface Integrals option for the relevant inletsurface. Results obtained for the different cases areillustrated in Table 3.6. Comparison between airflow rate prediction usingNetwork and CFD modelsTable 4 shows a comparison between airflow ratepredicted by the Network and CFD models. The negativeFig. 4. Airflow pattern, presented by contours of velocity magnitudes, for the cases modelled in this study.www.manaraa.com

ARTICLE IN PRESSO.S. Asfour, M.B. Gadi / Building and Environment 42 (2007) 4079–40854084Table 3Airflow rate prediction, using CFD modelCaseOpeningQn (kg/s)CaseOpeningQn (kg/s)11234123.19 3.192.98 2.983.30 3.3041234121.78 1.781.69 1.691.92 1.922356Table 4Discrepancy percentage between estimated and modelled airflow rate forcases 1–6CaseOpeningQn (kg/s)(CFD)Qn (kg/s)(Network)11234121234123.19 3.192.98 2.983.30 3.301.78 1.781.69 1.691.92 1.922.94 2.942.90 2.903.26 3.261.88 1.881.87 1.872.14 2.1423456Discrepancy(%) 5. Skew observed in the hex-map and sub-map meshes used in thecase of 451 wind direction.Table 5Discrepancy percentage between estimated and modelled airflow rate inthe case of oblique wind direction, after changing mesh typeCaseOpeningQn (kg/s)(CFD)Qn (kg/s)(Network)Discrepancy(%)41234121.84 1.841.83 1.832.26 2.261.88 1.881.87 1.872.14 2.14 2.2 5.6 10.75 11.5sign indicates that the modelled airflow rate is less than thecalculated one. In general, it shows that a good agreementhas been achieved. Discrepancy percentage observed isusually acceptable in airflow rate prediction, which is givenas a snapshot and measured in kg/s.However, discrepancy percentage, in general, is higher inthe case of oblique wind direction. CFD code used here isthe same of that one used with the normal wind direction.However, it has been found that the use of hex-map andsub-amp meshes in the oblique wind direction results in aless mesh quality. This is because tilted walls resulted inhigh angular skew between the edges of mesh cells, about0.7, which is considered high. For example, an excellentmesh has a skew less than 0.25 [9] (Fig. 5).Therefore, mesh type of the ambient air zone has beenchanged to the tetrahedral mesh, which has more flexibilityin meshing such geometries. Hex-map mesh is still used forroom interior. Table 5 shows a recalculation of thecomparison held between results obtained from bothmathematical and CFD models in the case of obliquewind direction. Change of mesh type seems to have a goodeffect, as discrepancy percentage has been significantlyreduced in all the three cases.The differences observed between airflow rates predictedby the Network and CFD models at both wind directionscan be justified by many reasons. One of them is theapproximation method used in simulating air velocity6 2.25.3profile, as explained in Section 5. This is because windinduced airflow rate is dependant on wind velocity. As thesquare of air velocity is used in the estimation of thispressure difference, any error in air velocity results in alarger error in airflow rate value.Another reason can be the approximation of themathematical procedure used. This is, on one hand,because the used wind pressure coefficient values areaveraged over the whole specified building face, and nota specified position on it. Liddament [5] highlighted thispoint and told that accurate evaluation of this parameter(i.e. pressure coefficient) is one of the most difficult aspectsof air infiltration modelling. In addition, these data weregenerated in wind tunnel experiments, where air density,and therefore pressure, is affected by air temperature,which is not the case in the isothermal CFD simulationcarried out. On the other hand, air infiltration model usedcontains many assumptions to enable the estimation ofairflow rate through a reasonable mathematical process.On the opposite, CFD considers the different values of airpressure on the opening, and calculate airflow rate as asummation their product.One more reason is related to the pressure coefficientdata used in the case of oblique wind direction. Airpressure distribution around a solid model changes if it isprovided with openings. In the case of 45o wind incidence,www.manaraa.com

ARTICLE IN PRESSO.S. Asfour, M.B. Gadi / Building and Environment 42 (2007) 4079–40854085of the building two windward faces, values of pressurecoefficient for the two windward faces in Fluent 5.5 havebeen found to be different, i.e., 0.28 for the face that hasthe window, and 0.32 for the other one. This can benoted in Fig. 6, showing contours of static pressure.7. ConclusionThis paper has compared the use of CFD and Networkmodels for wind-induced ventilation prediction in buildings. It compares the calculated airflow rate using themathematical Network model and the modelled one usingFluent 5.5 program. Many cases with a variety in buildinggeometries and wind directions have been considered.Results obtained support the use of the proposed CFDcode for wind-induced natural ventilation in buildings, as agood agreement has been achieved. This can be recommended as a validation method for studies that have noaccess to laboratory or full-scale testing facilities. Thisstudy has also revealed that the chose of mesh type andsize, in addition to the domain size are critical parametersin three-dimensional CFD modelling.Fig. 6. Contours of pressure magnitude (Pa) showing different pressuredistributions in the modelled cases and the standard solid ones.building model has two windward faces. In the case of asolid model, average pressure coefficient over these windward faces is the same for the square cases, and has arelatively small difference in the rectangular one. In thecase of placing an opening at one of these two windwardfaces, it is expected that air pressure distribution willchange, and there will be no more balance at itsdistribution on these windward faces. This is becausewindward with solid geometry receives more wind deflection on it. On the opposite, windward with an openingreceives less wind deflection, and therefore wind pressurewill be less too. For example, pressure coefficient values onthe windward faces in case of 451 wind direction and asquare building form are the same, i.e. 0.35, as presentedin the standard data used in this study. This is true for solidmodels. However, in the case of placing an opening in anyReferences[1] CIBSE. Natural ventilation in non-domestic buildings. London:CIBSE; 1997.[2] Bahadori M, Haghighat F. Passive cooling in hot arid regions indeveloping countries by employing domed roofs and reducing thetemperature of internal spaces. Building and Environment1985;20(2):103–13.[3] Chow WK. Wind-induced indoor-air flow in a high-rise buildingadjacent to a vertical wall. Applied Energy 2004;77:225–34.[4] Awbi H. Ventilation of buildings, second ed. London: Spon Press;2003.[5] Liddament MW. Air infiltration calculation techniques—an application guide. Bracknell. Coventry: The Air Infiltration and VentilationCentre; 1986.[6] British Standards Institution. BS 5925: code of practice for ventilationprinciples and designing for natural ventilation. London: BSI; 1991.[7] CIBSE., fifth ed. CIBSE guide, vol. A. London: CIBSE; 1988.[8] Kindamgen J, Krauss G, Depecker P. Effects of roof shapes on windinduced air motion inside buildings. Building and Environment1997;32(1):1–11.[9] Fluent Inc. Gambit 1.3 program help menu; 1988.www.manaraa.com

The CFD software used i s Fluent 5.5. Comparison between the predicted and simulated airflow rate is suggested as a validation method of the implemented CFD code, while the common practice is to compare CFD outputs to wind tunnel or full-scale . Both implemented CFD and Network models are briefly explained below. This followed by the .

Related Documents:

refrigerator & freezer . service manual (cfd units) model: cfd-1rr . cfd-2rr . cfd-3rr . cfd-1ff . cfd-2ff . cfd-3ff . 1 table of contents

430 allocation to elianto cfd o&m 20,577.32 440 allocation to trillium west cfd o&m 27,267.00 450 allocation to west park cfd o&m 70,008.22 460 allocation to festival ranch cfd o&m 177,790.54 480 allocation to tartesso west cfd o&m 27,809.17 481 allocation to anthem sun valley cfd o&

Emphasis is on comparing CFD results, not comparison to experiment CFD Solvers: BCFD, CFD , GGNS Grids: JAXA (D), ANSA (E), VGRID (C) Turbulence Models: Spalart-Allmaras (SA), SA-QCR, SA-RC-QCR Principal results: Different CFD codes on same/similar meshes with same turbulence model generate similar results

The first objective is to make comparison between the three CFD software which consists of ANSYS Fluent, Star-CCM and IESVE Mcroflo according to CFD . In terms of simulation results from the three baseline models, ANSYS Fluent is conclusively recommended for CFD modeling of complicated indoor fluid environment compared with Star-CCM

Trindade, B.C. and J. Vasconcelos. 2011. "Comparison of CFD with Reservoir Routing Model Predictions for Stormwater Ponds." Journal of Water Management Modeling R241-05. doi: 10.14796/JWMM.R241-05. . This may help to explain why such CFD models are seldom used for detention pond analysis. Comparative Study of a CFD Model and a Reservoir .

an opportunity to use Computational Fluid Dynamics (CFD) to supplement experimental data with . and flow across the physics of interest. In order for CFD to be used confidently, these models must be validated against expected literature correlations as well as experimental data. . Wall Nusselt Number Comparison for CFD and Experimental .

EFD, CFD, and UA laboratories. EFD labs were improved and UA was introduced. Complementary CFD labs were also introduced using an advanced research code modified for limited user options. From 1999 to 2002, the research CFD code was replaced by the commercial CFD software (FLUENT) and refinements were made and

The classical approach to public administration, derived from Weber, Wilson and Taylor, largely . Classical and Modern Approaches to Public Administration * Polya Katsamunska is a Ph.D., associate professor at the Public Administration and Regional Development of UNWE, e-mail: polya_katsamunska@yahoo.com. 75 Articles is really impressive and yet "almost no national government would argue .