Final ENSO2124

This article appeared in a journal published by Elsevier. The attached copy is furnished to the author for internal non-commercial research and education use, including for instruction at the authors institution and sharing with colleagues. Other uses, including reproduction and distribution, or selling or licensing copies, or posting to personal, institutional or third party websites are prohibited. In most cases authors are permitted to post their version of the article (e.g. in Word or Tex form) to their personal website or institutional repository. Authors requiring further information regarding Elsevier’s archiving and manuscript policies are encouraged to visit: http://www.elsevier.com/copyright Author's personal copy Environmental Modelling & Software 24 (2009) 389–410 Contents lists available at ScienceDirect Environmental Modelling & Software journal homepage: www.elsevier.com/locate/envsoft Influences of wind flow over heritage sites: A case study of the wind environment over the Giza Plateau in Egypt Ashraf S. Hussein a, *,1, Hisham El-Shishiny b, 2 a b Faculty of Computer and Information Sciences, Ain Shams University, 11566, Cairo, Egypt IBM Center for Advanced Studies in Cairo, IBM Cairo Technology Development Center, IBM WTC, Egypt Branch, P.O.B. 166, El-Ahram, Giza, Egypt a r t i c l e i n f o a b s t r a c t Article history: Received 11 March 2008 Received in revised form 6 August 2008 Accepted 10 August 2008 Available online 8 October 2008 Wind influences are important environmental factors that cause deterioration of historical heritage sites. This work presents a computational framework for investigating the influences of the wind flow over such sites. The wind flow is considered to be fully turbulent, isothermal and incompressible. The present framework employs three-dimensional Reynolds Averaged Navier–Stokes (RANS) equations along with multi-block approach and non-conformal meshes to perform wind flow simulations over such sites with complex geometry. As a case study, the influences of wind flow over the ‘‘Giza Plateau’’, one of the most important Egyptian historical heritage sites, were studied for the Northwest wind (at average wind speed over the year) and the Southwest windstorms. The study addresses the less understood, yet important, influences of the wind flow structure on the site and its famous monuments: the Pyramids and the Great Sphinx. Qualitative and quantitative treatments of the results are carried out for estimating the wind loading on the different monuments within the plateau. Particular attention was paid to the Great Sphinx to investigate its most vulnerable parts to the wind, as one of the critical environmental factors that cause erosion of this colossal statue. Ó 2008 Elsevier Ltd. All rights reserved. Keywords: Environmental modeling Wind engineering Wind over heritage sites Modeling and simulation Software applications Computational fluid dynamics 1. Introduction There are many environmental factors that contribute to the deterioration of heritage sites including wind action, groundwater, humidity, pollution, vibrations and seismic shocks (Conservation of the Sphinx project). Wind flow is one of the most denudation factors that cause erosion of monuments especially when they are located in open areas. There are four major key studies that can be used to investigate the wind structure over complex site terrains: (1) field study; (2) theoretical study; (3) experimental study; (4) numerical study. Both field studies and experiments are usually expensive, and often provide data that are not satisfactorily detailed. On the other hand, the numerical studies often called Computational Fluid Dynamics (CFD) have become definitely one of the most useful tools for studying the atmospheric wind environments in recent years primarily due to the increasing power of the new generations of computers. * Corresponding author. Tel.: þ20106066143; fax: þ20226827574. E-mail addresses: ashrafh@acm.org (A.S. Hussein), shishiny@eg.ibm.com (H. El-Shishiny). 1 Currently visiting scientist at the IBM Center for Advanced Studies (CAS) in Cairo. 2 Tel.: þ20106066078/235361461; fax: þ20235362505. 1364-8152/$ – see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.envsoft.2008.08.002 CFD simulations can provide information on all flow parameters in the entire modeling domain. Moreover, a reliable numerical evaluation of the interaction between wind and the Giza Plateau topography and monuments can be achieved with CFD modeling in an efficient manner (Chattot, 2002). On the other hand, CFD could also be the most superior and cost-effective way over the traditional wind tunnel studies, for the cases of wind flow investigations over complex terrains, in the absence of accurate physical models. To our knowledge, there were no three-dimensional CFD-based simulations performed before to study wind environment over a heritage site. Consequently, our literature survey is limited to the use of CFD for simulating wind environment over complex terrain including the implementation issues and the guidelines for these wind environment predictions. Comprehensive literature reviews on the use of CFD for the simulations of atmospheric environments over complex terrain have been published (Stathopoulos, 1997; Reichrath and Davies, 2002; Blocken and Carmeliet, 2004; Bitsuamlak et al., 2004; Meroney, 2004; Franke et al., 2004). The implementation issues in 3D wind flow prediction over complex terrain have been reviewed and analyzed by Prospathopoulos and Voutsinas (2006) to provide guidelines with respect to the definition of the appropriate boundary conditions and the construction of an adequate and effective computational mesh when a Reynolds Average Navier–Stokes (RANS) solver is implemented. Recently, Author's personal copy 390 A.S. Hussein, H. El-Shishiny / Environmental Modelling & Software 24 (2009) 389–410 Proposed Framework Boundary Surface Conditions Volume Reconstructor Mesh Mesh to CAD Generator Converter Case CFD Simulations Converter Visualization Management Management Parallelization Root User Interface HPC Facility Fig. 1. The proposed framework architecture. Blocken et al. (2007a,b) reviewed and investigated the wall function problems in CFD simulation of the atmospheric boundary layer to provide suggestions to improve the simulations by modifying the wall function roughness. Frank et al. (2007) presented the best practice guidelines for the CFD simulations of flows over microscale obstacle-accommodating meteorological models using RANS equations. On the other hand, Garcia and Boulanger (2006) successfully simulated the low altitude wind flow over the Mount Saint Helens in the United States using the high resolution terrain data provided by the NASA SRTM database by employing the open source CFD solver called ‘‘OpenFoam’’ to solve the Reynolds average Navier–Stokes equations with standard k–3 turbulence model. Bechmann et al. (2007) simulated the neutral atmospheric Site Terrain Level Contour Maps Preprocessing Surface Reconstruction I. Surface Mesh wind over natural terrain employing a hybrid RANS/LES model to reduce the computational cost of the traditional LES model. Finally, Hussein and El-Shishiny (2007) investigated the low speed wind environment around the Great Sphinx, isolated from the surrounding area, in the case of Northwest wind employing parallel Navier–Stokes solver. Although their investigations shed some light on the Sphinx vulnerable parts which experience maximum wind pressure and friction, the study did not consider the Sphinx surrounding topography and the corresponding wind interactions. The goal of this work is to provide an efficient framework to simulate the wind environment over heritage sites taking the Giza Plateau, as a case study, to investigate the wind influences on its monuments with particular attention to the Great Sphinx as the most vulnerable monument to wind effects (Conservation of the Sphinx project). The wind environment around the terrain of the plateau was simulated for the cases of the Northwest wind, the most prevailing wind during the year (Mortensen et al., 2006), at average wind speed in addition to a case representing the Southwest windstorms. The outline of the rest of the paper is as follows: Section 2 describes the present framework. The simulation and the visualization managements are outlined in Section 3. The case study of the wind environment over the Giza Plateau is presented in Section 4. Finally, Section 5 gives the conclusion of this work and an overview of the future work. 2. Framework 2.1. Architecture The present framework comprises several parts as depicted in Fig. 1, which include central CFD simulation and visualization Surface Cleaning and Reduction Surface Mesh to CAD Format Export Primitive Objects & Mesh Parameters 3D Scanned Objects F. Surface Mesh Volume Mesh Generation Definition of Boundaries & Continuum F. Volume Mesh Physical and Numerical Parameters Boundary Conditions Visualization Parameters Images & Animations Visualization System Selected Field Variables Simulation Case File CFD Simulation System Mesh Converter Testing for Grid Independent Solution Simulation Results Fig. 2. CFD simulation and visualization pipeline. Author's personal copy A.S. Hussein, H. El-Shishiny / Environmental Modelling & Software 24 (2009) 389–410 Outflow Y X Inflow Z 391 On the other hand, the visualization toolkit should give users the ability to apply advanced visualization and analysis techniques to their simulation results. These techniques can be applied to help users gain new insights into results. Also, visualization management should provide a full set of tools for manipulating, transforming, processing, realizing, rendering and animating data and allow for visualization and analysis methods based on points, lines, areas, volumes, images or geometric primitives in any combination. These core toolkits should be object-oriented cross platform software packages and have an easy interface for the classes that perform different approaches and algorithms, rendering and interaction techniques. 2.2. Simulation and visualization pipeline Fig. 3. Applications of the boundary conditions. managements. These managements were designed to support uniprocessor or run on based on distributed memory parallel architectures. Also, the components of the present framework include the surface reconstruction tool, surface mesh to CAD format converter, volume grid generator, application of boundary and continuum condition tool, simulation case converter, the user interface with a variety of default functionalities, and a data root for communication with dedicated high performance computing facilities. The proposed framework is based on two primary toolkits; the CFD simulations toolkit for implementation of a diversity of CFD approaches and algorithms, based on the problem under consideration, and the scientific visualization toolkit for the implementation of several visualization algorithms including both scalar and vector data visualization techniques. The CFD toolkit uses the finite volume numerical scheme to solve the equations of conservation on different types of volume meshes. Also, it supports the multi-block approach along with non-conformal unstructured meshes to be suitable for simulating wind environment over sites with complex geometries and diversity between the occupied objects’ sizes. Mapper Abstract Geometric Representations Renderer Selected Technical Data Image Sequences Filter Display/VR Visual Presentation Simulation Results Fig. 4. The visualization model. The solid model of any heritage site is the primary input for the simulation and visualization pipeline as shown in Fig. 2. The solid model of heritage site includes the terrain of the site, the heritage objects, statues and their associated complexes, and primitive objects, if any. Therefore, the pipeline starts by acquiring the heritage site topographic level contour map or terrain data to automatically generate the solid model of the site terrain only, without monuments. Therefore, a post-processing procedure is applied on the topographic map to remove the monuments then the level contour map is converted into a gray scale contour map for the surface altitude. This map is sampled into equidistant nodes and its cloud of points is generated. In this work, Autodesk MAYA modeling software version 8.0 (Maya software) was used to generate the terrain cloud of points. On the other hand, the sophisticated objects like sculptures and statues are reconstructed by acquiring their 3D scanned cloud of points or by obtaining their dimensions and details from contour maps and satellite images. All the 3D clouds are constructed together or individually using (Farouk et al., 2003) or (Hussein and Abdel Aziz, 2007) based on the cloud size. The reconstructed surfaces are efficiently decimated preserving the curvature details by suppressing points at low curvature regions. Depending on the topography of the terrain, about 15–30% (based on our own experience) of the terrain data can be removed, while keeping a low approximation error, by sequentially eliminating triangles that produce the smallest error until an approximation threshold is reached. In order to obtain efficient surface mesh to be useful for further processing, the resulting surface mesh is cleaned by removing very small faces and edges, and sharp angles. In our case, a locally developed tool is used to clean the surface meshes. The clean surface mesh is then converted into CAD format in order to be readable by the volume mesh generator. In this work, the cleaned surface mesh is converted to ACIS standard CAD format. The grid generator imports the converted surface mesh in ACIS format. According to the original site map, the user completes the topology of the site by adding the primitive objects, like tombs, and concatenates them with the imported surface mesh. The computational domain is then constructed according to the guidelines provided by (Prospathopoulos and Voutsinas, 2006; Frank et al., 2007; Bechmann et al., 2007) based on the maximum height of the heritage site ‘‘Hmax’’. Generation of a computational mesh for such geometry is extremely complicated because of the diversity between the dimensions of heritage objects and the terrain geometric features. Consequently, the multi-block approach along with non-conformal unstructured meshes is adopted considering the issues raised in Liu and Shyy (1996). Accordingly, the computational domain is split into multiple subdomains. This technique allows treating the different sub-domains with different volume meshes without preserving identical node locations at the interfaces between these different sub-domains. Author's personal copy 392 A.S. Hussein, H. El-Shishiny / Environmental Modelling & Software 24 (2009) 389–410 Fig. 5. Overview of the Giza Plateau (The Giza archives project). This permits simulating the wind environment over a heritage site for different wind direction angles by just rotating the volume meshes of the inner sub-domains and regenerating only the volume mesh of outer sub-domain. The rotation angle around the vertical axis ‘‘Y-axis’’ is q ¼ (mean wind direction angle (deg.)  360 ). The main factor affecting the minimum mesh spacing near the wall surface is the constraint of placing the first computational node off the wall, at least ks away from the wall, where ks is the sand-grain roughness. This may lead to a contradiction with the concept of good mesh resolution. In this work, the approach of testing several meshes is adopted until obtaining independent mesh solution. Following the generation of the volume meshes of the different sub-domains, the continuum of each mesh is defined and the boundary conditions are applied for each volume mesh. All of the boundaries of the volume-meshes, except the bottom, of the inner sub-domains are marked as ‘‘Interface-Boundary’’. On the other hand, the monuments under particular attention should also be marked to facilitate studying the field variables and loads on their surfaces. In the present framework, Gambit 2.2.3 (Gambit software) is used as a volume mesh generator. Before exporting the simulation case file for the CFD simulation manager, the simulation case file is converted to be imported successfully within the simulation manager. From within the CFD simulation manager, the meshes of different sub-domains are loaded sequentially and the mesh interfaces are defined and the boundary conditions are adjusted for all the boundaries of the global volume mesh. Then, the partial differential equations of conservation are solved. The flow is considered to be fully turbulent, isothermal (neutral atmospheric conditions are assumed) and incompressible. Since this work mainly deals with the prediction of the statistically steady mean flow and turbulence in heritage sites environments for situations with neutral stratification, the Reynolds-averaged Navier–Stokes (RANS) equations are employed considering different types of k–3 turbulence models. CFD simulation requires iterating the solution of the fluid flow equations, starting from initial guess, till it is converged. The flow variables are recalculated in each of the iterations until the equations are solved up to a user-specified threshold. The termination criterion is usually based on the residuals of the corresponding equations. These residuals should tend to zero. Scaling of the residuals is done with the residuals after the first iteration. In many industrial applications, it is confirmed that the residuals are well Fig. 6. Distribution of pressure coefficient along the centerline of the mounted cube. Fig. 7. Velocity profiles at different positions along the bottom centerline. Author's personal copy 393 A.S. Hussein, H. El-Shishiny / Environmental Modelling & Software 24 (2009) 389–410 Table 2 The wake reattachment lengths for the numerical and experimental results Turbulence model Reattachment length/H (from leeward face of the cube) Error% (with reference to the full scale experiment of Hunt, 1982) Standard k–3 RNG k–3 Realizable k–3 Experiment of Castro and Robins (1977) Experiment of Hunt (1982) 1.81 1.8 1.8 3.5–4.5 28 28 28 1.4 Fig. 8. Turbulence intensity profiles at different positions along the bottom centerline. below the general accepted tolerance of 1 103 (Frank et al., 2007), which is in general too high to have a converged solution. In this work, a reduction of the residuals of at least five orders of magnitude is used. For validation purposes of turbulence models, much lower thresholds are used. In addition to the residuals, the target variables are also recorded. If these variables are constant then the solution is confirmed to be converged. Finally, the converged simulation results are stored and the field variables under consideration are exported to the visualization manager. OpenFOAM CFD toolkit version 1.4.1 (OpenFoam software) is served as CFD simulations manager. OpenFOAM is the Open Source Field Operation and Manipulation (OpenFOAM) Cþþ libraries used primarily to create applications on different platforms. OpenFoam code uses the finite volume numerical scheme to solve the equations of conservation on both uniprocessor and it can run on workstation clusters. The IBM OpenDX version 4.4.4 (OpenDX software) is used as the visualization management component which facilitates the development of visual programs to visualize the field variables under consideration via its flexible DX file format. This feature permits the customization of the available visualization capabilities through different scenes to provide the user with high quality images and animations from the visual study. Open Visualization Data Explorer supports uniprocessor and symmetric multiprocessor workstations. It also runs as a distributed application on workstation clusters. 3. Simulation and visualization managements 3.1. Simulation management 3.1.1. Basic equations The wind environment over heritage sites is governed by the conservation laws of mass and momentum. The flow is assumed to be three-dimensional, incompressible Newtonian fluid with Table 1 The coefficients of lift and drag on the mounted cube surfaces obtained with different turbulence models Turbulence model CL Error % CD Error % Standard k–3 RNG k–3 Realizable k–3 Experiment of Castro and Robins (1977) 0.715 0.561 0.641 0.5 43 12.2 28.2 0.727 0.843 0.754 0.78 6.7 8.07 3.3 Fig. 9. Giza Plateau topography. (a) Sample of the Giza Plateau satellite images (The Giza archives project). (b) Level contour map of the Plateau topography (The Giza archives project). Author's personal copy 394 A.S. Hussein, H. El-Shishiny / Environmental Modelling & Software 24 (2009) 389–410 Fig. 10. Details of the Great Sphinx and its surrounding area. (a) Top view of the Great Sphinx and its associated complex (AERAGRAM, 2004). (b) Contour map of the Great Sphinx front view (Lehner and James, 1985). constant density. Based on the above assumptions, the Reynolds average Navier–Stokes (RANS) equations are as (Moshfegh and Nyiredy, 2004), vui ¼ 0 vxi (1) close Eq. (2), a turbulence model is necessary for the turbulent stress tensor. Most turbulence models for RANS follow the eddyviscosity hypothesis by Boussinesq (1877) where the turbulent stresses are assumed to be the product of the fluid strain and an eddy-viscosity, mt. By adopting this concept of Boussinesq, 2 3 r0 u~ i u~ j ¼ 2mt Sij þ dij r0 k;  vui vuj ui 1 v ~i u ~j  pdij þ 2mSij  r0 u ¼ fi  þ r0 vxj vt vxj (2) where the Einstein summation notation is used. ui (i ¼ 1, 2 and 3) denote the velocity components, p is the pressure and m is the dynamic molecular viscosity. fi represents body forces, which act on the mass of the fluid. The pressure, the body forces and the velocity components are decomposed into resolved parts denoted by ,). In order to overbar (,) and unresolved parts denoted by tilde (e   Sij ¼ 0:5 vui =vxj þ vuj =vxi (3) where Sij is the mean rate of strain tensor and the eddy-viscosity, mt, is a positive scalar. The most common way to derive, mt, is to relate it to turbulent kinetic energy, k, and its dissipation rate, 3, by means of the following equation mt ¼ r0 Cm k2 =3. Thus, two addi~ and tional transport equations for the turbulent kinetic energy,k, ~ for the dissipation of turbulent kinetic energy, 3, need to be solved. From these two quantities a length scale (l ¼ k3=2 =3) and a velocity scale (v ¼ k1=2 ) are constructed and the eddy-viscosity computed. Fig. 11. Rendered solid model of the Giza Plateau. (a) The solid model of the Giza Plateau emphasizing the main monuments. (b) Details of the Great Sphinx solid model. Author's personal copy 395 A.S. Hussein, H. El-Shishiny / Environmental Modelling & Software 24 (2009) 389–410 50 Hmax 3.1.2. Turbulence models In this work, three turbulence models were considered; the standard k–3 model (SKE), the renormalization group k–3 model (RNG) and the realizable k–3 model. The standard k–3 model The transport equation for k and 3 after modeling may be expressed as (Launder and Spalding, 1974a,b) Giza Plateau Solid Model      v uj k m vk r0 þ r0 m þ t Vk þ mt S2  r0 3 ¼ V, sk vxj vt 50 Hmax 80 Hmax Giza Plateau Sub-domain    v u3 v3 r0 þ r0 j ¼ V, 50 Hmax vt Giza Plateau Surrounding Sub-domain Hmax 26 Hmax Sphinx Sub-domain Fig. 12. Extensions of the computational domain indicating the different sub-domains. vxj mþ    mt 3 V3 þ C31 mt S2 C32 r0 3 s3 k (4) (5) Typical values for coefficients are ðCm ; sk ; s3 ; C31 ; C32 Þ ¼ ð0:09; 1:00; 1:30; 1:42; 1:92Þ (Launder and Spalding, 1974a). The RNG k–3 model The coefficients of the k–3 are determined from a number of case studies of simple turbulent flows. The k–3 thus has a limited scope of applicability, which yields poor performance in cases with complex flows. This poor performance is suspected to be due to inaccuracies in the 3-equation. The renormalization group k–3 model (RNG) introduces an additional term in the 3-equation, which improves its performance. The basic idea is to systematically filter out the small-scale turbulence to a degree that the remaining scales can be resolved. This is done by the parameter, h, which is the Fig. 13. The computational mesh of the Sphinx sub-domain. (a) Details of Sphinx-complex surface mesh. (b) Details of the Giza Plateau surface mesh. (c) Cross-sectional view in the Sphinx sub-domain volume mesh. (d) Cross-sectional view in the Giza Plateau sub-domain volume mesh. Author's personal copy 396 A.S. Hussein, H. El-Shishiny / Environmental Modelling & Software 24 (2009) 389–410 a 25 22.5 20 % Cells 17.5 15 12.5 10 7.5 5 0 0.00 0.00 0.02 0.05 0.07 0.10 0.12 0.15 0.17 0.19 0.22 0.24 0.27 0.29 0.31 0.34 0.36 0.39 0.41 0.44 0.46 0.48 2.5 Cell Equiangle Skew b 70 60 % Cells 50 40 30 20 516.67 483.33 450.00 383.33 416.67 350.00 283.33 316.67 250.00 183.33 216.67 150.00 83.33 116.67 50.00 16.67 0 0.00 10 Cell Volume (m3) c Cell Volume Fig. 14. Quality of the computational mesh. (a) Histogram of mesh skewness (cell equiangle skew). (b) Histogram of mesh cell volume (m3). (c) Continuity in cell size at the interfaces of the inner sub-domains. ratio between the time scales of the turbulence and the mean flow (Launder and Spalding, 1974a; Rohdin and Moshfegh, 2007). The k– 3 equation is identical with the standard k–3 and 3 equation is given by    v u3 v3 r0 þ r0 j ¼ V, vt vxj mþ    mt 3 V3 þ C31 mt S2 C3*2 r0 3 s3 k (6) Fig. 15. (a) Distribution of the observation points over the Giza Plateau. (b) Distribution of the observation points around the Great Sphinx. (c) Observation points over the chest and the head of the statue. Author's personal copy A.S. Hussein, H. El-Shishiny / Environmental Modelling & Software 24 (2009) 389–410 Table 3 The wall roughness values used in the present computations Zone Wall roughness (m) Sphinx body Giza Plateau Giza Plateau surrounding 0.04 0.2 0.5 where h ¼ Sk=3, h0 ¼ 4:38, b0 ¼ 0:012, Cm ¼ 0.0845, C31 ¼1.42, C32 ¼ 1.68, sk ¼ s3 ¼ 0:7178 and C3*2 ¼ C32 þCm h3 ð1 h=h0 Þ=ð1þ bh3 Þ. Realizable k–3 model ‘‘Realizable’’ means that values of variables, e.g. negative normal stresses, can be removed from the predictions by numerical clipping. To achieve the realizability effect the Cm is no longer constant but a function of the turbulence fields, mean strain and rotation rates. The k equation is identical with the standard k–3 and 3 equation is given by (Rohdin and Moshfegh, 2007) r0      v uj 3 m r 32 v3 m þ t V3 þ C1 Sr0 3  C2 0 pffiffiffiffiffi ¼ V, þ r0 s vxj vt k þ n3 3 (7) where C1 ¼ max½0:43; h=ðh þ 5Þ, C2 ¼ 1.0, sk ¼ 1.0, s3 ¼ 1.2. Generally, the standard k–3 model exhibits very good stability but both the renormalization group k–3 model (RNG) and the realizable k–3 model (RKE) have better characteristics in high Reynolds number flows (Frank et al., 2007). All of the three models were studied in the present work with the modified standard wall functions to consider the wall roughness according to (Bechmann et al., 2007; Blocken et al., 2007a,b; Frank et al., 2007; Prospathopoulos and Voutsinas, 2006). All the discretized equations are solved in a segregated manner with the Semi-Implicit Method for pressure-Linked Equations (SIMPLE) algorithm. The second order Upwind difference scheme is used for spatial discretization. Potential flow solution is used as an initial solution of turbulent flow (Garcia and Boulanger, 2006). 3.1.3. Boundary conditions In the following description of the applied boundary conditions, the computational mesh is oriented so that the X-axis is aligned with the main wind direction while the Y-axis indicates the height above the ground as shown in Fig. 3. Bottom boundary The ground boundary conditions for the velocity are formulated on the basis of the ‘‘wall-function’’ theory (Launder and Spalding, 1974a,b; Blocken et al., 2007b). The basic principles of this theory are the condition of local energy balance and the assumption that the flow near the wall is similar to a one-dimensional Couette flow (Prospathopoulos and Voutsinas, 2006). These standard wall functions were 45 40 35 % Cells 30 25 20 15 10 modified for roughness based on experiments with sand-grain roughness (Blocken et al., 2007b). The sand-grain roughness ks is computed from the relationship with the aerodynamic roughness length y0 as (Prospathopoulos and Voutsinas, 2006) 9:793 y Cs 0 ks ¼ Wall Y plus Fig. 16. Histogram of wall yþ, Northwest wind over the Giza Plateau at 5 m/s. (8) where Cs is the roughness constant and generally set at its default value for sand-grain roughened pipes and channels: 0.5 (Blocken et al., 2007b). The main problem in setting the value of the wall roughness is the contradiction between the good mesh resolutions near the rough wall and the constraint of placing the first computational node off the wall at least ks away from the wall. Inflow boundary The vertical distribution of the inflow velocity follows the logarithmic law of a fully developed boundary layer (Zbigniew, 1989; Prospathopoulos and Voutsinas, 2006). uðx; y; zÞ ¼ u* ðx; zÞ ln k  y y0  (9) where k is the von Karman constant (0.40–0.42) and the friction velocity, u*, and the aerodynamic roughness length, y0, are input data and y is the height above the ground. The friction velocity was adjusted to give a turbulence intensity of 10%. The turbulence kinetic energy and dissipation rate at the inflow boundary are estimated as k ¼ 2 3 Uavg I 2 3 ¼ Cm3=4 k3=2 l (10) (11) where Uavg is the mean velocity at the inlet, I is the turbulence intensity and l is the turbulence length scale. Outflow boundary The open boundary condition based on the constant static pressure boundary condition was used (Hargreaves and Wright, 2007). With this boundary condition the derivatives of all other flow variables are forced to vanish, corresponding to a fully developed flow. Therefore, this boundary was placed far enough from the plateau area to not have any fluid entering into the computational domain through this boundary. Top boundary As the top boundary was placed far enough outside the boundary layer, symmetry boundary conditions were applied to enforce a parallel flow by forcing the velocity component normal to the boundary to vanish and prescribing zero normal derivatives for all other flow variables (Frank et al., 2007; Hargreaves and Wright, 2007). Lateral boundaries In the present framework, symmetry boundary conditions were used at the lateral boundaries as the approach flow direction is parallel to them (Frank et al., 2007; Hargreaves and Wright, 2007). Finally, the surface pressure and friction coefficients are computed as (Chattot, 2002): Cp ¼ p  pN 2 0:5rN UN (12) Cf ¼ sw 2 0:5rN UN (13) 88.94 86.83 84.72 82.61 80.51 78.40 76.29 74.18 72.08 69.97 67.86 65.75 63.65 61.54 59.43 57.32 55.22 53.11 0 51.08 5 397 Author's personal copy 398 A.S. Hussein, H. El-Shishiny / Environmental Modelling & Software 24 (2009) 389–410 Fig. 17. Surface friction coefficient distribution, Northwest wind over the Giza Plateau at 5 m/s. Distribution of the surface friction coefficient over the (a) Giza Plateau, (b) Great Sphinx and its associated complex. where pN , rN and UN are the undisturbed free stream pressure, density and velocity, respectively, and sw is the wall shear stress. 3.2. Visualization management High fidelity, three-dimensional CFD simulations of wind over heritage sites involve complex and severe phenomena especially when these objects include complex details. Therefore, investigation of these phenomena and their effects on the monuments can be augmented by advanced visualization techniques that provide new insights and better understanding. The visualization model, shown in Fig. 4, starts with the source process which reads the simulation results. The filter process either selects or samples data. It is also used to extract spatial or data value domains. In the mapping step, the selected visualization method or superimposed methods convert data into abstract visual representations. The visualization management provides several visualization algorithms that implement different types of mappings; each of them has its specific capabilities to represent simulation results. In most cases, the mapping process leads to a collection of geometric primitives such as triangles, lines, point clouds, etc. This can be combined with textures and material properties of surfaces for improved visual impression. The rendering process is responsible for converting the 3D scene into images to be displayed. In this framework, the interaction activities are separated into multiple feedback loops. The innermost loop feeding back into the rendering process allows modifying user position and orientation thus enabling free roaming within the 3D scene. Also, this loop allows zooming into specific details of the scene. In order to compare different simulation results, the corresponding data files need to be stored in memory enabling quick switching between them. In the next outer feedback loop, a user can interact with the parameters of a selected mapper. Typical interactions with mappers are e.g. the repositioning of a cutting plane, the seeding points of the stream ribbons within the flow field, the definition of a new isovalue of an isosurface, etc. The next outer loop allows interactions with the filtering steps. Filter parameters enable selecting sub-domains with their Author's personal copy A.S. Hussein, H. El-Shishiny / Environmental Modelling & Software 24 (2009) 389–410 399 Fig. 18. Surface pressure coefficient distribution, Northwest wind over the Giza Plateau at 5 m/s. Distribution of the surface pressure coefficient over the (a) Giza Plateau, (b) Great Sphinx and its associated complex. associated values. Finally, the visualized data can be changed via feedback into the input of simulation results. From the inner to the outer loops, the timing requirements for the reaction of the system become less demanding. When moving through a scene or interacting with other animation parameters the system reacts ideally within a 1/30 of a second to give a real time behavior. Table 4 Wind loads (N) and moments (N m) on the main monuments, Northwest wind over the Giza Plateau at 345 and 5 m/s Main monuments X-force Z-force X-moment Y-moment Z-moment Khufu Khafre Menkaure Great Sphinx K1 K2 K3 M1 M2 M3 193,223 155,558 35,962 480 4006 5726 3637 1307 3978 3766 7641 21,458 4412 394 534 2523 2862 1453 265 1415 38,460,149 47,559,910 27,580,231 6097 1,449,757 1,201,768 782,039 11,580,000 13,374,973 9,571,074 37,293,740 1.06e þ 08 35,038,242 9228 68,319 680,923 1,109,849 1,562,676 4,116,282 4,305,740 1.07e þ 08 29,150,139 459,071 9710 8,119,482 9,370,088 9,510,428 1,932,646 1,834,843 1,113,653 4. Case study: Simulation of wind environment over the Giza Plateau in Egypt Giza Plateau, in Egypt, is one of the oldest heritage sites which is mainly known because of the three pyramids and the Great Sphinx that were built there during the 4th Dynasty some 4500 years ago (Fig. 5). The site is located on the west bank of the Nile, opposite Egypt’s modern day capital of Cairo, 17 km to the North of Saqqara and some 8 km to the South of Abu Rawash; it is part of the Northernmost extension of the necropolis of ancient Memphis. The site covers an area of about 2  2 km (Site management of the Giza Plateau). The pyramids and the other monuments of the Giza Plateau are among the great world treasures. The most prevailing wind direction over the year is the Northwest direction at about 345 and wind speed of approximately 5 m/ s (Mortensen et al., 2006). Therefore, simulations were performed for these conditions. In addition, another case for the wind at a direction angle of 225 and 16 m/s was considered to simulate the wind environment resulting from Southwest windstorms. All computations were carried out on 3 GHz Pentium IV PC platform with a 2 GB Ram. Author's personal copy 400 A.S. Hussein, H. El-Shishiny / Environmental Modelling & Software 24 (2009) 389–410 Fig. 19. Streamlines north of the Khafre cause-way seeded from four points at 10 m above the ground level, Northwest wind over the Giza Plateau at 5 m/s. 4.1. Calibration of CFD simulation management Since atmospheric measurements of wind over complex terrain are sparse, the standard test case of a wall mounted cube has been selected to calibrate the simulation setup and to investigate the effect of turbulence model selection. For this case, the flow is highly complex with several recirculation zones and contrary to streamlined bodies (like Askervein hill (Taylor and Teunissen, 1983)). Separation points are clearly defined by the sharp edges of the cube. Simulation of a wall mounted cube is therefore a good test case for validating the turbulence model’s ability to capture the physics of the separated shear layers and complex wake regions behind bluff bodies in atmospheric conditions. The wind tunnel measurements by Castro and Robins (1977) have been selected to validate simulation results. The considered turbulence models were Standard k–3 model (SKE), RNG k–3 model (RNG) and Realizable k–3 model (RKE) described earlier. The computational domain extended 10H vertically, 5H laterally on both sides, 6H upstream and 18H downstream the cube. The dimensions of the cube were the same as the experiments i.e., 0.2 m each side (H). The mean inlet velocity and turbulence intensity profiles were extracted from Castro and Robins (1977). First, the computational mesh independency test Fig. 20. The recirculation zones downwind the pyramids, Northwest wind over the Giza Plateau at 5 m/s. Author's personal copy A.S. Hussein, H. El-Shishiny / Environmental Modelling & Software 24 (2009) 389–410 401 Fig. 21. Wind structure around the Great Sphinx, Northwest wind over the Giza Plateau at 5 m/s. was carried out to obtain the proper mesh for further computations, considering the concept of good mesh (the value of yþ, which is the non-dimensional distance normal to the wall, for the first cell adjacent to the wall is between 30 and 60). Grid adaptation was one of the main steps in validating the model for analysis. The grid was adapted several times based on the values of yþ of the near wall cells until the yþ values were within the limiting range as described before. This was followed by the adaptation based on the change in cell volume and velocity gradients in the flow until a grid independent solution was obtained. To test grid independence, the wall shear stress along the bottom centerline was plotted after each converged solution. After each grid adaptation step, a comparison was made in the shear stress distribution of two successive solutions, and the procedure continued until these plots were matched. By investigating the computed pressure coefficient distributions, along the centerline of the cube using the considered turbulence models, it is found that the pressure coefficient distribution resulting from the three turbulence models are relatively in good agreement with the experimental data (Castro and Robins, Fig. 22. The reverse flow at the right edge of the Sphinx cavity, Northwest wind over the Giza Plateau at 5 m/s. Author's personal copy 402 A.S. Hussein, H. El-Shishiny / Environmental Modelling & Software 24 (2009) 389–410 45 40 35 % Cells 30 25 20 15 10 169.21 165.70 162.19 158.67 155.16 151.65 148.14 144.63 141.13 137.62 134.10 130.59 127.08 123.57 120.05 116.54 113.03 109.51 0 106.13 5 Wall Y plus Fig. 23. Histogram of wall yþ, Southwest windstorm over the Giza Plateau at 16 m/s. 1977) and the numerical results of Thomas and Williams (1999) using Large Eddy Simulation Model. However, the RKE model provided more accurate results than that of the other considered models as shown in Fig. 6. The predicted velocity profiles along the centerline were compared favorably with the experimental data of Castro and Robins (1977) and the numerical results of Thomas and Williams (1999) as shown in Fig. 7. It can be noticed that the experimental results do not show negative velocity components, while the numerical results have negative velocity components in a region near the roof which is 11% of the cube height. Apart from this, the numerical results are in good agreement with the experimental measurements, and interestingly the standard k–3 model results are closer to the experimental observations. Also, the wake produced using RNG model is the thickest among all models. The computed turbulence intensity profiles were compared with the experimental data of Castro and Robins (1977) and the numerical results of Thomas and Williams (1999) as shown in Fig. 8. The computed results using both SKE and RKE are in good agreement with the experimental data while both the computed results using the RNG model and the numerical results (Thomas and Williams, Fig. 24. Surface friction coefficient distribution, Southwest windstorm over the Giza Plateau at 16 m/s. Distribution of the surface friction coefficient over the (a) Giza Plateau, (b) Great Sphinx and its associated complex. Author's personal copy A.S. Hussein, H. El-Shishiny / Environmental Modelling & Software 24 (2009) 389–410 403 Fig. 25. Surface pressure coefficient distribution, Southwest windstorm over the Giza Plateau at 16 m/s. Distribution of surface pressure coefficient over the (a) Giza Plateau, (b) Great Sphinx and its associated complex. 1999) clearly deviate from the experimental data. The computed wind loading is presented in Table 1 in comparison with the experimental data of Castro and Robins (1977). It can be noticed that the RKE model predicts the closest values to the experimental data. Finally, the computed centerline wake reattachment length was compared to the experimental data of Castro and Robins (1977) Table 5 Wind loads (N) and moments (N m) on the main monuments, Southwest windstorm over the Giza Plateau at 225 and 16 m/s Main monuments X-force Khufu Khafre Menkaure Great Sphinx K1 K2 K3 M1 M2 M3 Z-force X-moment Y-moment 724,319 449,967 6.2 þ 08 3.5e þ 08 2,107,901 469,340 4.05e þ 08 9.12e þ 08 379,148 33,899 73,759,856 1.3e þ 08 25,392 23,876 169,766 50,643 29,832 10,632 22,902,854 10,055,479 19,993 5400 38,656,536 7,350,932 24,887 10,113 49,750,200 10,658,871 20,379 7093 8,222,328 304,025 1587 1,423,412 10,221,459 23,297 20,786 557 194,413 8,759,074 Z-moment 48,816,860 4.2e þ 08 1.7e þ 08 440,660 6,828,790 15,475,105 22,598,780 22,635,064 5,490,422 621,799 and Hunt (1982) and all of the considered turbulence models provided almost the same values as shown in Table 2. Based on the results obtained from the above investigation, it is clear that RKE is expected to be the best turbulence model selection for the wind environment simulation, for this work, from both stability and accuracy point of view, which is in accordance with the conclusion given in Frank et al. (2007). 4.2. Simulation of wind environment over the Giza Plateau A solid model of the Giza Plateau including its terrain and monuments has been constructed to provide the simulation setup with a good approximation to the physical domain. The terrain data were acquired from satellite images and contour maps of (AERAGRAM, 2006; The Giza archives project) shown in Fig. 9. Following the pipeline described earlier, the site terrain was represented by 40,000 points and 15% of the terrain data was removed, while keeping a low approximation error, by sequentially eliminating triangles that produce the smallest error until an approximation threshold was reached. Author's personal copy 404 A.S. Hussein, H. El-Shishiny / Environmental Modelling & Software 24 (2009) 389–410 Fig. 26. The recirculation zones downwind the pyramids, Southwest windstorm over the Giza Plateau at 16 m/s. The dimensions of the pyramids were acquired accurately from (Legon, 1979; Cole, 1925; Site management of the Giza Plateau) and constructed considering the topographic maps (Fig. 9). The great pyramid of Khufu was originally 479 ft (146 m) in height, it still stands at an awe-inspiring 449 ft (137 m). Its sides measure 755.76, 755.41, 755.87 and 756.08 ft from west to South in clockwise order. Khafre’s pyramid measures 707 ft (216 m) on each side and was originally 471 ft (143 m) high. For Menkaure’s pyramid, each side measures 356 ft (109 m) and the structure’s completed height was 218 ft (66 m). The dimensions of the surrounding Queens’ pyramids Fig. 27. The structure of recirculation at the Khafre cause-way, Southwest wind over the Giza Plateau at 16 m/s. Author's personal copy A.S. Hussein, H. El-Shishiny / Environmental Modelling & Software 24 (2009) 389–410 405 Fig. 28. The wind structure around the Great Sphinx, Southwest windstorm over the Giza Plateau at 16 m/s. were obtained from AERAGRAM (2006), Legon (1979), and The Giza archives project. The final stage in constructing the Giza Plateau solid model was the construction of the solid model of the Sphinx complex (the sculpture and its associated stone temples). The dimensions of these monuments were obtained from the level contour maps of (AERAGRAM, 2004; Fig. 10a), and the Sphinx map project (Lehner and James, 1985; Fig. 10b). The Sphinx is approximately 240 ft (73 m) long and 66 ft (20 m) high. The details of the Sphinx cavity 1 0.8 were obtained from AERAGRAM (2004) and the real photos acquired from the site. Maya software version 8.0 (Maya software) was used in constructing the solid model of the Sphinx and its surrounding objects. The main geometric features of the Great Sphinx are represented by 12,000 points. Then, all of the constructed monuments were concatenated with the constructed terrain of the plateau. As a final step, a geometry cleanup process was performed to remove short and hard edges, holes, sharp and large angles, and small faces. The rendered model of the plateau is North Northwest Southwest Pressure Coefficient 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 Observation Points Fig. 29. Distribution of the surface friction coefficient at the observation points of the Great Sphinx complex. Author's personal copy 406 A.S. Hussein, H. El-Shishiny / Environmental Modelling & Software 24 (2009) 389–410 0.016 0.014 North Northwest Southwest Friction Coefficient 0.012 0.01 0.008 0.006 0.004 0.002 0 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 Observation Points Fig. 30. Distribution of the surface pressure coefficient at the observation points of the Great Sphinx complex. illustrated in Fig. 11a, indicating the main monuments of the plateau while Fig. 11b shows the rendered model of the Great Sphinx. The Giza Plateau domain of interest extends 1.66 km  1.66 km and the size of the computational domain was selected considering the guidelines provided by Prospathopoulos and Voutsinas (2006), Frank et al. (2007), and Bechmann et al. (2007) based on the maximum height of the plateau ‘‘Hmax’’. The upstream boundary was placed 50Hmax from the upstream edge of the plateau while the downstream boundary was placed 80Hmax. The lateral boundaries were both placed 50Hmax from the sides of the plateau and the ceiling boundary was extended to 25Hmax from the highest point in the plateau. Fig. 12 illustrates the extension of the computational domain from the solid model of the Giza Plateau for the case of North wind. For other wind directions, the solid model of the plateau is rotated around the Y-axis in order to get the situation of inclined wind. Generation of a computational mesh for such geometry is extremely complicated because the surface of the plateau is highly irregular and the dimensions of the Great Sphinx complex are one order of magnitude less than the other monuments of the plateau. Consequently, the multi-block approach along with nonconformal unstructured meshes was adopted considering the issues raised in Liu and Shyy (1996). Accordingly, the computational domain was split into three sub-domains; the Sphinx complex sub-domain, the Giza Plateau sub-domain and the Giza Plateau surrounding sub-domain as shown in Fig. 12. This technique allows treating the different sub-domains with different volume meshes without preserving identical node locations at the interfaces between these different sub-domains. Also, it permits simulating the wind environment over the plateau for different wind direction angles by just rotating the volume meshes of the inner sub-domains and regenerating only the volume mesh of outer sub-domain ‘‘sub-domain 3’’. For each sub-domain, Gambit size function was used (Gambit software) to control the size of mesh intervals for edges and mesh elements for faces or volumes. Other meshing software like Gridgen (Gridgen software) could possibly do the same task of meshing complex domains. For the Sphinx complex sub-domain, as the surface of the Great Sphinx and its complex include complex geometry features, a curvature size function (Gambit software) was attached to the surface of the Great Sphinx and its associated complex to refine the mesh elements near the surface. Three volume meshes have been used to examine mesh independence considering the condition of surface roughness (Blocken et al., 2007b). Mesh 1 had 300,000 tetrahedral cells, mesh 2 had 820,000 tetrahedral cells and mesh 3 had 1,411,375 tetrahedral cells. With the exception of mesh 1 (300,000 cells), the other two meshes presented relevant information in the occupied sub-domain, but mesh 2 showed unsatisfactory results in the proximity of the Sphinx surface. Therefore, results from computations using mesh 3 (1,411,375 cells) are presented in this paper. For this mesh, the minimum mesh size was 0.05 m with a 1.24 growth rate. For the Giza Plateau sub-domain, a curvature size function (Gambit software) is attached to the surface of the Giza Plateau providing a mesh with 700,000 tetrahedral cells. The minimum mesh size was 0.225 m with a 1.26 growth rate. Finally, The Giza Plateau surrounding sub-domain was discretized using 200,000 tetrahedral cells using a fixed size function (Gambit software) with a minimum mesh size of 1.2 m and 1.3 growth-rate. Mesh cell-size continuity was considered at the interfaces between the three subdomains to efficiently minimize the numerical dissipation at these interfaces. Fig. 13a and b shows the surface meshes of the Sphinx and Giza Plateau sub-domains. Fig. 13c and d shows cross sectional views in the volume meshes colored by the cell size value to demonstrate the mesh quality especially near the monuments and their complexes. In order to demonstrate the quality of the entire computational mesh, for the three sub-domains, the mesh skewness (cell equiangle skew) histogram and the mesh cell-size histogram are presented in Fig. 14a and b, respectively. In addition, the mesh cell-size continuity between Sphinx sub-domain and Giza Plateau sub-domain is demonstrated in Fig. 14c. 4.2.1. Selection of observation points A set of observation points was placed on the surface of the plateau, as shown in Fig. 15a, to explore the wind environment features around the monuments. A total of 50 observation points were placed around the different monuments and their associated complexes; 26 of them were concentrated on the Great Sphinx and its associated complex (Fig. 15b) to emphasize the wind effects on the body of the statue. As indicated in Hawass (2004) the Sphinx weak areas are the shoulders, the chest and the top of the hunches Author's personal copy A.S. Hussein, H. El-Shishiny / Environmental Modelling & Software 24 (2009) 389–410 407 (bumps). So, some of the observation points were selected to be placed on these weak areas as shown in Fig. 15b and c. 4.2.2. Simulation of Northwest wind environment around the Giza Plateau The Northwest wind environment around the Giza Plateau was investigated by simulating the Northwest wind blowing over the plateau with a wind direction angle of 345 at 5 m/s (which is the most prevailing wind direction and average wind speed over the year). This simulation was conducted considering the surface roughness, described in Table 3, and employing the Realizable k–3 model (RKE that was selected based on the investigation outlined in Section 4.1) along with the logarithmic wall function with a sandgrain-based roughness modification of Blocken et al. (2007a,b). Computations were conducted until a steady state was reached and took about 20 hours. In order to study the fitness of the computational mesh for boundary layer modeling, the wall yþ values were examined as shown in Fig. 16. The histogram of the wall yþ demonstrates that the wall yþ values are within the acceptable range according to Blocken et al. (2007a). Fig. 17a shows the distribution of the surface friction coefficient on the Giza Plateau. The faces of the pyramids which face the wind, western cemetery and North of Khafre causeway experience relatively high friction coefficient values. For the Great Sphinx, the left shoulder and the top of the hunches experience high wind friction and the surface friction coefficient reaches its maximum value at the head, the face and the top trunk as shown in Fig. 17b. It is interesting to point out that, from this figure, the Sphinx weak areas (as indicated in Hawass, 2004) which are the left shoulder and the top of the hunches are exposed to high wind friction. The resulting pressure coefficient distribution over the surface of the plateau was computed from the predicted wind environment as shown in Fig. 18a. The distribution exhibits maximum values on the North sides of the three great pyramids and high altitude areas of the plateau and the left side of the Sphinx all of which face the wind as shown in Fig. 18b. As shown, the Sphinx vulnerable areas to wind pressure include the left shoulder, the left part of the head, the left part of the trunk and the left thigh. This is critical for the head as it causes a moment on the thorax. The resulting wind loads and moments, as defined in the world axes, acting on the different monuments within the plateau are depicted in Table 4. The tunneling effect North of the Khafre cause-way due to the variation of the plateau altitude was captured using four streamlines seeded from points at 10 m above the surface of the plateau as shown in Fig. 19. The wind structure downwind the pyramids is depicted in Fig. 20 demonstrating generally large recirculation regions with relatively higher momentum downwind Khufu and Khafre pyramids than that of Mankaure’s pyramid. Fig. 21 shows the wind structure around the Sphinx-complex which is very sophisticated especially at the pavement to the left of the Sphinx, the Sphinx cavity and front of the Sphinx’s ‘‘Dream-Stela’’. By investigating the wind structure around the Sphinx-cavity, it is interesting to point out that the Sphinx-cavity edge, to the right of the Sphinx, experiences west flow near the ground as it moves from a relatively high pressure region to a lower one, which is contrary to the main stream around this region, as shown in Fig. 22. Fig. 31. Extracted vortex cores in the Sphinx complex (blue for the extracted using velocity gradient eigenmodes and red for the extracted using vorticity vector). (a) North wind over the Giza Plateau at 5 m/s, (b) Northwest wind over the Giza Plateau at 5 m/s, and (c) Southwest windstorm over the Giza Plateau at 16 m/s. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article). 4.2.3. Simulation of Southwest windstorm environment around the Giza Plateau The Southwest wind environment around the Giza Plateau was investigated by simulating the Southwest windstorms blowing over the plateau with a wind direction angle of 225 at 16 m/s (which is the max wind speed recorded during windstorms, generally blowing from the Southwest, in the area). This simulation was conducted considering the surface roughness, described in Table 3, Author's personal copy 408 A.S. Hussein, H. El-Shishiny / Environmental Modelling & Software 24 (2009) 389–410 Fig. 32. Resemblance between the predicted vulnerable areas to wind effects and the actual week areas of the Great Sphinx. (a) Recent images of the Great Sphinx (current situation). (b) Average shear stress distribution. (c) Average pressure coefficient distribution. and employing the Realizable k–3 model (RKE-that was selected based on the investigations outlined in Section 4.1) along with the logarithmic wall function with a sand-grain-based roughness modification of Blocken et al. (2007a,b). Computations were conducted until a steady state was reached and took about 25 hours. The resulting distribution of the wall yþ was studied as shown in Fig. 23. The wall yþ histogram revealed that the yþ values are within the acceptable range according to Blocken et al. (2007a). Fig. 24a shows the distribution of the surface friction coefficient on the Giza Plateau. Both Mankaure’s complex and Sphinx-complex experience a relatively high wind friction. For the Great Sphinx, the friction coefficient reaches its maximum value at the head, the left shoulder, and the trunk while the top of the hunches, the right shoulder and hand experience a relatively high wind friction as shown in Fig. 24b which is in good agreement with that indicated in Hawass (2004). The resulting pressure coefficient distribution over the surface of the plateau was computed form the predicted wind environment Author's personal copy A.S. Hussein, H. El-Shishiny / Environmental Modelling & Software 24 (2009) 389–410 409 as shown in Fig. 25a. The distribution exhibits maximum values on the pyramids’ sides which face the windstorm and the right side of the Sphinx as shown in Fig. 25b. As shown, the Sphinx vulnerable areas to wind pressure include the right part of the head and thorax, the right part of the trunk, and the right thigh. This is critical for the head only (which experiences pressure coefficient of about 0.9–1) because it causes a moment on the thorax. The resulting wind loads and moments, defined in the world axes, acting on the different monuments within the plateau are depicted in Table 5. As shown the resulting windstorm forces and moments are one order of magnitude larger than that of the normal average wind. The wind structure downwind the pyramids is depicted in Fig. 26 which is dominated by large recirculation regions. Due to the velocity of the windstorm, the wind structure near the ground is very sensitive to the altitude variation as shown downwind Khufu’s complex and to the North of Khafre’s cause-way (Fig. 27). Also, the wind structure around the Sphinx-complex is very sophisticated as shown in Fig. 28. Northwest wind and Southwest windstorms. The Northwest wind and Southwest windstorm cases were simulated using the RKE turbulence model considering surface roughness. Qualitative and quantitative treatments of the results were carried out estimating the wind loads and moments on the different monuments within the plateau. Particular attention was paid to the Great Sphinx to investigate its most vulnerable areas to wind, as one of the critical denudation factors that cause the deterioration of the statue. The present work may give more insight into the effect of the wind around the Giza plateau when developing a global plan for conserving and protecting the site. A similar framework as the one presented in this paper can be considered for studying the effect of wind flow over other heritage sites around the world. As a future work, the mixture of air and dust will be considered as a two phase flow problem for simulating strong dust bearing wind over the Giza plateau. 4.3. Discussion This work is part of an R&D project conducted at the IBM Center for Advanced Studies (CAS) in Cairo. The authors would like to thank Mr. Ahmed Ali for his help in constructing the solid model of the Giza Plateau. For the North and Northwest wind, and Southwest windstorm, the surface friction coefficient at the observation points within the Sphinx complex is shown in Fig. 29. Generally, as the wind changes its direction from the North to the Northwest, the surface friction on the left part of the Sphinx complex increases significantly. This increase reaches an average of 25% at the left shoulder and thigh, and the head of the Great Sphinx. On the other hand, the Southwest windstorm increases the surface friction on the right part of the Sphinx complex especially the right rearward part of the Sphinx. The pressure coefficient at the observation points, for the considered simulation cases is shown in Fig. 30. The left thigh of the statue experiences a significant increase in the surface pressure by more than 100% as the wind changes its direction from the North to the Northwest. The Sphinx top of the hunches is the most vulnerable area to wind pressure in the case of the Southwest windstorm. The vortex cores within the Sphinx-complex were detected as shown in Fig. 31. Both velocity gradient eigenmodes and vorticity vector (Guo et al., 2004) were used to detect these vortex cores. For all the simulated cases, the Sphinx cavity is dominated by recirculation regions that may contribute to sediment accumulation within the cavity. Finally, to examine the average wind influences on the Great Sphinx, as the most vulnerable monument (in this site) to the wind effects, the average pressure and shear stress distributions were evaluated on the surface of the statue due to the most prevailing range of wind and storms directions. The resulting average distributions at the weak areas of the statue are visually compared with the recent images representing the current situation of the statue as shown in Fig. 32. It is interesting to note that all the weak areas indicated in Hawass (2004) experience high pressure or high shear stress or both that may contribute to the denudation of the statue. 5. Conclusions We have described a computational framework for investigating wind flow influences on heritage sites. This framework combines both computational fluid dynamics and visualization capabilities to simulate the wind environment over such sites and visually investigate the resulting environment. The CFD simulation management employs three-dimensional Reynolds Averaged Navier–Stokes (RANS) equations along with multi-block approach and non-conformal meshes to perform wind flow simulations over heritage sites with complex geometry. As a case study, the wind environment around the Giza Plateau was explored for the cases of Acknowledgements References AERAGRAM, 2004. Newsletter of the Ancient Egypt Research Associates 7(1). AERAGRAM, 2006. Newsletter of the Ancient Egypt Research Associates 8(1). Bechmann, A., Sørensen, N.N., Johansen, J., Vinther, S., Nielsen, B.S., Botha, P., 2007. Hybrid RANS/LES method for high Reynolds numbers, applied to atmospheric flow over complex terrain. The science of making torque from wind. Journal of Physics: Conference Series 75, 012054, doi:10.1088/1742-6596/75/1/012054. Bitsuamlak, G.T., Stathopoulos, T., Bedard, C., 2004. Numerical evaluation of wind flow over complex terrain: review. Journal of Aerospace Engineering 17 (4), 135–145, doi:10.1061/(ASCE)0893-1321200417. 4(135). Blocken, B., Carmeliet, J., 2004. A review of wind-driven rain research in building science. Journal of Wind Engineering and Industrial Aerodynamics 92 (13), 1079–1130, doi:10.1016/j.jweia.2004.06.003. Blocken, B., Carmeliet, J., Stathopoulos, T., 2007a. CFD evaluation of wind speed conditions in passages between parallel buildings – effect of wall-function roughness modifications for the atmospheric boundary layer flow. Journal of Wind Engineering and Industrial Aerodynamics 95 (9–11), 941–962. Blocken, B., Stathopoulos, T., Carmeliet, J., 2007b. CFD simulation of the atmospheric boundary layer: wall function problems. Atmospheric Environment 41 (2), 238–252. Boussinesq, J., 1877. Théorie de lécoulement tourbillant. Mem. Présentés par Divers Savants Acad. Science Institution Paris 23, 46–50. Castro, I.P., Robins, A.G., 1977. The flow around a surface-mounted cube in uniform and turbulent streams. Journal of Fluid Mechanics 79 (2), 307–335. Chattot, J.J., 2002. Computational Aerodynamics and Fluid Dynamics: An Introduction. Springer Publications, ISBN 3540434941. Cole, J.H., 1925. Determination of the exact size and orientation of the great pyramid. Published by Government Press. <http://www.theglobaleducationproject.org/ egypt/studyguide/gpmath.php>. Conservation of the Sphinx project: an environmental strategy, 2007. Bibliotheca Alexandrina, Center for Documentation of Cultural and Natural Heritage, Strategic Approach to Egypt’s Cultural Heritage Part IV <www.cultnat.org/ download/Pdfs/part_4/14-Conservation%20of%20the%20Sphinx%20Project.pdf>. Farouk, M., El-Rifai, I., El-Tayar, S., El-Shishiny, H., Hosny, M., El-Rayes, M., Gomes, J., Giordano, F., Rushmeier, H.E., Bernardini, F., Magerlein, K.A., 2003. Scanning and processing 3D objects for web display. In: proceedings of the 4th International Conference on 3-D Digital Imaging and Modeling, 310–317. Frank, J., Hellsten, A., Schlünzen, H., Carissimo, B., 2007. Best practice guideline for the CFD simulation of flows in the urban environment. COST action 732. Quality Assurance and Improvement of Meteorological Models. University of Hamburg, Meteorological Institute, Center of Marine and Atmospheric Sciences. Franke, J., Hirsch, C., Jensen, A.G., Krüs, H.W., Schatzmann, M., Westbury, P.S., Miles, S.D., Wisse, J.A., Wright, N.G., 2004. Recommendations on the use of CFD in wind engineering. COST Action C14, Impact of Wind and Storm on City Life Built Environment. In: van Beeck, J.P.A.J. (Ed.), Proceedings of the International Conference on Urban Wind Engineering and Building Aerodynamics. von Karman Institute, Sint-Genesius-Rode, Belgium. Gambit user’s guide, 2007. Version 2.2.3. <http://www.fluent.com/>. Garcia, M.J., Boulanger, P., 2006. Low altitude wind simulation over Mount Saint Helens using NASA SRTM digital terrain model. 3DPVT’06. Proceedings of the 3rd International Symposium on 3D Data Processing, Visualization, and Transmission, 535–542. Gridgen user’s manual, 2007. <http://www.pointwise.com/products/doc.shtml>. Author's personal copy 410 A.S. Hussein, H. El-Shishiny / Environmental Modelling & Software 24 (2009) 389–410 Guo, D., Evangelinos, C., Patrikalakis, N.M., 2004. Flow feature extraction in oceanographic visualization. In: Proceedings of the Computer Graphics International, doi:10.1109/CGI.2004.1309207 (CGI’04) 162–173. Hargreaves, D.M., Wright, N.G., 2007. On the use of the k–3 model in commercial CFD software to model the neutral atmospheric boundary layer. Journal of Wind Engineering and Industrial Aerodynamics 95 (5), 355–369, doi:10.1016/ j.jweia.2006.08.002. Hawass, Z., 2004. The Secrets of the Sphinx: Restoration Past and Present. The American University in Cairo Press. Hunt, A., 1982. Wind-tunnel measurements of surface pressures on cubic building models at several scales. Journal of Wind Engineering and Industrial Aerodynamics 10, 137–163. Hussein, A.S., Abdel Aziz, A.H., 2007. A framework for implicit surfaces reconstruction from large clouds of points. Proceedings of the 7th IEEE International Symposium on Signal Processing and Information Technology, pp. 611. Hussein, A.S., El-Shishiny, H., 2007. Modeling and simulation of low speed wind over the Great Sphinx. ISCC’07. Proceedings of the 2007 IEEE Symposium of Computers and Communications, 1027–1034, doi:10.1109/ISCC.2007.4381500. Launder, B.E., Spalding, D.B., 1974a. Mathematical Models of Turbulence. Academic Press, London. Launder, B.E., Spalding, D.B., 1974b. The numerical computation of turbulence flows. Computer Methods in Applied Mechanics and Engineering 3, 269–289. Legon, J.A., 1979. The plan of the Giza pyramids. Archaeological Reports 10(1). The Archaeology Society of Staten Island and the Staten Island Society, Archaeological Institute of America. <http://goodfelloweb.com/giza/jlpamp.html>. Lehner, M., James, A., 1985. Sphinx map project. <www.aeraweb.org/sphinx_map. asp>. Liu, J., Shyy, W., 1996. Assessment of grid interface treatments for multi-block incompressible viscous flow computation. Journal of Computers and Fluids 25 (8), 719–740, doi:10.1016/S0045-7930(96)00022-9. Maya user’s guide, 2007. Version 8.0. <www.autodisk.com/>. Meroney, R.N., 2004. Wind tunnel and numerical simulation of pollution dispersion: a hybrid approach. Working paper. Croucher Advanced Study Institute on Wind Tunnel Modeling, Hong Kong University of Science and Technology, 60 pp. Mortensen, N.G., Hansen, J.C., Badger, J., Jørgensen, B.H., Hasager, C.B., Paulsen, U.S., Hansen, O.F., Enevoldsen, K., Youssef, L.G., Said, U.S., Moussa, A.A., Mahmoud, M.A., Yousef, A.E., Awad, A.M., Ahmed, M.A., Sayed, M.A.M., Korany, M.H., Tarad, M.A., 2006. Wind atlas for Egypt: measurements, micro- and mesoscale modelling. In: Proceedings of the 2006 European Wind Energy Conference and Exhibition, Athens, Greece, February 27 to March 2, 10 pp. Moshfegh, B., Nyiredy, R., 2004. Comparing RANS models for flow and thermal analysis of pin fin heat sinks. In: Proceedings of 15th Australasian Fluid Mechanics Conference, Sydney, AustraliaAFMC00253. OpenDx user’s guide, 2007. Version 4.4.4. <http://www.opendx.org/>. OpenFoam user’s guide, 2007. Version 1.4.1. <http://www.opencfd.co.uk/>. Prospathopoulos, J., Voutsinas, S.G., 2006. Implementation issues in 3D wind flow predictions over complex terrain. Journal of Solar Energy Engineering 128 (4), 539–553, doi:10.1115/1.2346702. Reichrath, S., Davies, T.W., 2002. Using CFD to model the internal climate of greenhouses: past, present and future. Agronomie 22 (1), 3–19. Rohdin, P., Moshfegh, B., 2007. Numerical predictions of indoor climate in large industrial premises. A comparison between different k–3 models supported by field measurements. Journal of Building and Environment 42 (11), 3872–3882, doi:10.1016/j.buildenv.2006.11.005. Site management of the Giza Plateau: case study, 2007. Bibliotheca Alexandrina, Center for Documentation of Cultural and Natural Heritage, Strategic Approach to Egypt’s Cultural Heritage Part III, <www.cultnat.org/download/Pdfs/part_3/ 8-Site%20Mgmt%20of%20Giza%20Plat.pdf>. Stathopoulos, T., 1997. Computational wind engineering: past achievements and future challenges. Journal of Wind Engineering and Industrial Aerodynamics 67–68, 509–532, doi:10.1016/S0167-6105(97)00097-4. Taylor, P.A., Teunissen, H.W., 1983. Askervein ’82: Report on the September/October 1982 experiment to study boundary-layer flow over Askervein, south uist. Atmos. Environ. Service, Downsview, Ontario, Technical Report MSRB-83-8. The Giza archives project, 2007. <http://www.gizapyramids.org/code/emuseum. asp?newpage¼index>. Thomas, T.G., Williams, J.J.R., 1999. Simulation of skewed turbulent flow past a surface mounted cube. Journal of Wind Engineering and Industrial Aerodynamics 81 (1–3), 347–360. Zbigniew, S., 1989. Structure of the Atmospheric Boundary Layer. Prentice Hall, Englewood Cliffs, NJ. Ashraf S. Hussein was born in Giza, Egypt, in 1970. He received his B.Sc., M.Sc. and Ph.D. degrees in Aerospace Engineering from the Cairo University, Egypt in 1992, 1996 and 1999, respectively. The major fields of his studies are Modeling and Simulation, and Computer Visualization. During 1992-1996, he worked as a research engineer in the Department of Aerospace Engineering, Cairo University. He joined IBM Cairo Technology Development Centre in 1996 and IBM Cairo Center for Advanced Studies in 2005. In addition, he joined the Faculty of Computer and Information Sciences, Ain Shams University, Egypt in 2001. He is currently Associate Professor in the Department of Scientific Computing, Faculty of Computer and Information Sciences, Ain Shams University. He is also a Visiting Scientist, IBM Center for Advanced Studies in Cairo. His research areas cover Modeling and Simulation, Computational Techniques in Science and Engineering, Computer Graphics and Visualization, and High Performance Computing Applications. He has more than 35 published scientific papers in international conferences and journals. He has participated in more than 15 R&D and incubation projects. Dr. Hussein served as a member in numerous technical committees of international scientific conferences. He is a member of the editorial board of one technical journal and a professional member of the ACM. Hisham El-Shishiny is currently Senior Scientist and Manager of Advanced Technology and Center for Advanced Studies at IBM Technology Development Center in Cairo. His current research interests are Modeling & Simulation, Image Processing & Visualization, HPC, Knowledge Management and Data Mining. During 1981–1989 and 2001 to present, Dr. Shishiny was an adjunct university professor in Computer Science at American University in Cairo (AUC), Alexandria University and Cairo University, and an Independent Consultant for the United Nations. Dr. H. El-Shishiny has authored over 40 published scientific papers and has several issued patents and filed patent applications. He is on the editorial board of two technical journals and a professional member of the ACM. He is the recipient of IBM 1st and 2nd Plateau Invention Achievement Awards and IBM Outstanding Innovation Award. Dr. Shishiny received his Ph.D. and M.Sc. in Computer Science with highest honors from Nancy University in France and his B.Sc. in Electronics and Communications Engineering from Cairo University. He has also completed a two year post-graduate program in Operations Research and conducted post-graduate studies in Solid State Sciences at the American University in Cairo.