Mathematical Modeling and Numerical Analysis of Thermal Distribution in Arch Dams considering Solar Radiation Effect

The effect of solar radiation on thermal distribution in thin high arch dams is investigated. The differential equation governing thermal behavior of mass concrete in three-dimensional space is solved applying appropriate boundary conditions. Solar radiation is implemented considering the dam face direction relative to the sun, the slop relative to horizon, the region cloud cover, and the surrounding topography. It has been observed that solar radiation changes the surface temperature drastically and leads to nonuniform temperature distribution. Solar radiation effects should be considered in thermal transient analysis of thin arch dams.


Introduction
Performance evaluation of high arch dams is of interest in the field of structural engineering due to their complicated design methodology and catastrophic financial and loss consequences due to their failure. Thus, identification, description, and screening of the all parameters affecting the safety of arch dams are important. Dam self-weight, hydrostatic load, silt pressure, thermal loads, and seismic loads (if any) are the most important loads which are usually considered in design and analysis processes. Sources of the thermal loads in concrete dams can be categorized into the three main groups, that is, seasonal temperature and climate variation, mass concrete hydration at early ages, and solar radiation on the exposed faces. It is essential to realize the sensitivity of the thermal distribution within the thin high arch dams with respect to these parameters.
In general there are two approaches in thermal analysis of concrete dam, that is, thermal analysis of thin arch dams in which the assumption of linear temperature distribution through the dam thickness may be acceptable and solar radiation can have significant effect on structural responses of the whole dam body and thermal analysis of the relatively thick sections (i.e., arch-gravity, gravity, and RCC dams) in which nonlinear temperature distribution through the body thickness is expected. In the second case, concrete temperature near the dam face (especially downstream face) changes quickly with respect to the air and water temperature while the temperature inside the body has no meaningful variation. Thus, the solar radiation may affect only the temperature of the exterior areas with no considerable effect on the interior sections [1].
Agullo et al. [2] measured the temperature of each node in numerical model of a concrete dam using explicit finite difference method. They used energy conversion factor for radiation modeling and Newton's law for simulation of the thermal convection. Meyer and Mouvet [3] investigated thermal behavior of a concrete gravity dam in Switzerland utilizing finite element method. They found that the air temperature fluctuation is not a lonely good criterion for estimation of concrete dam's temperature. Solar radiation increases the whole dam's temperature of about 4 ∘ C and the upstream water could regulate the temperature fluctuation amplitude. Castellanos and Marin [4] studied the effect of various parameters on thermal behavior of the concrete arch and gravity dams using a three-dimensional (3D) finite element model. B. Zhu [5] proposed an analytical formula as boundary condition of heat transfer analysis in concrete dams in order to estimate the water temperature through the deep reservoirs. Sheibany and Ghaemian [6] developed 2 The Scientific World Journal a detailed 3D finite element model of an arch dam for estimation of thermal stresses. They found that probable cracks occur in a very narrow region of the downstream face which is mainly due to thermal loads in comparison with self-weight and hydrostatic loads. Ardito et al. [7] studied diagnostic damage analysis of a concrete arch-gravity dam based on seasonal hydrostatic loading considering thermal effects. Jin et al. [8] investigated a practical procedure for predicting nonuniform temperature on the exposed faces of the arch dams due to variation in solar radiation and shading effects. They found that the temperature distribution on the upstream face is mainly due to the reservoir temperature condition, while on the downstream face the temperature distribution is mainly influenced by the solar radiation. They showed that the nonuniform temperature has a significant effect on the surface thermal stresses and crack propagation. K. H. Zheng and B. Y. Zheng [9] studied coupled thermal-creep-stress problem for high arch dams using finite element method. Yang et al. [10] developed a finite element program for thermal analysis of mass concrete embedded with double-layer staggered heterogeneous cooling water pipes based on the equivalent equation of heat conduction including the effect of cooling water pipes and hydration heat of concrete. Comparing calculated results with actual measured data from a monolith of an arch dam in China, their numerical model was proven to be effective in simulating accurate temperature variations of mass concrete. Yangbo et al. [11] proposed two fast algorithms to overcome the huge and time-consuming finite element analysis on temperature control and thermal stress fields of concrete dams' construction. One was overall planning algorithm and the other was incomplete Cholesky conjugate gradient algorithm. By applying these two to the temperature forecasts, it was shown that the calculated results are accurate. Fujun et al. [12] studied the cracking reasons of concrete overflow dam of Hadashan Hydro Project by developing a finite element method to simulate the temperature distribution and thermal stress analysis during the construction period. Results have shown that the crack of the concrete overflow dam is a temperature crack due to combined action of the internal thermal gradient and the external restraints.
In the present paper, impact of the solar radiation on thermal distribution in a high thin arch dam is investigated by solving the 3D differential equation governing the phenomenon within the mass concrete and finally the results are compared with those obtained from recorded data from the dam site.

Mathematical Model for Heat Transfer
2.1. Thermal Conduction. Conduction (or heat conduction) is a mode of transfer of energy within and between bodies of matter due to a temperature gradient. Governing equation of thermal conduction in continuous environments can be achieved according to the conservation principle of thermal energy on constant arbitrary volume, , surrounded by closed surface, . Implementation of heat transfer relation for a system including anisotropic material is as follows: where is density, is specific heat, is temperature of medium, is time, ⃗ is a tensor representing thermal conduction coefficient, and is internal heat generated per volume.

Thermal
Convection. The amount of heat transferred by convection is governed by Newton's cooling law as follows: where ℎ is convective coefficient, sb is temperature of surface boundaries, and is ambient temperature. Convective coefficient is a function of fluid velocity, , and its properties such as type of flow, viscosity, and surface roughness. There are many simplified formulas for calculating this coefficient such as the one proposed by Duffie and Beckman [13]:

Thermal Radiation.
Thermal radiation is an electromagnetic radiation generated by the thermal motion of charged particles in matter. All matter with a temperature greater than absolute zero emits thermal radiation. In the present analysis, radiation involved at the issue is in two forms, that is, radiation energy absorbed by the surfaces and electromagnetic energy released from the surface which is known as thermal radiation. The total amount of the absorbed energy is expressed as [14] = , where is solar absorptivity of surface and is total solar energy on surface. Thermal radiation caused by the temperature difference between surface and its environment can be achieved by Stefan-Boltzmann law: where is emissivity of surface, that is, the ratio of heat released by intended surface to heat released by a black body with the same temperature, and = 5.669 × 10 −8 Wm −2 K −4 is Stefan-Boltzmann constant. It should be noted that the wasted heat between surface and air due to radiation is not significant. Thus, (5) can be rewritten as a pseudolinear equation as follows: where ℎ is defined as Therefore, the radiation effects can be considered by adding ℎ to ℎ and obtaining new convective coefficient.

Boundary Condition for Heat Transfer Equation
The numerical solution of differential equation requires information as known quantities applied to the systems in the form of boundary conditions. Boundary conditions of three-dimensional conduction with convection are shown in Figure 1 schematically including imposed constant temperature on the boundary 1 , imposed heat flux on the boundary 2 , and the edge convection condition on the boundary 3 . Also, the environmental conditions in a typical arch dam that are important in heat transfer analysis and its detailed methodology are shown in Figure 2.  weather station. In the current study, the following data were collected from weather stations: (i) daily mean temperature, , (ii) daily maximum ( max ) and daily minimum temperatures ( min ) , (iii) the maximum monthly mean temperature, max , (iv) the minimum monthly mean temperature, min , (v) annual mean temperature, mean , and (vi) the highest and lowest recorded temperatures at the site.
If there is no information about the air temperature at the dam site, it is possible to use the sinusoidal function with the acceptable accuracy as where is the amplitude, is the annual mean temperature, is time in day, is time in day in which ( ) = , and is period of sine function, which here is 365 days. In this function can be obtained as in which, where is the number of days at month. Figure 3 compares the recorded and estimated air temperatures for one of the stations, in the current research, as a sample.

Reservoir
Temperature. Distribution of the temperature within the dam body is strongly affected by the reservoir water temperature. The best source for this information is the collected statistical data from the similar reservoirs. This type of information is collected by USBR and reported in engineering supplement number 34 [16]. The proposed method by Bofang [5] is another way for estimation of the reservoir temperature which is used in the present study (11). Bofang showed that the proposed method is in good agreement with the recitation values from China reservoirs and can be used for other reservoirs with small changes in the parameters: where is reservoir depth, is time in day, ( , ) is the reservoir temperature at time and depth , ( ) is the mean annual temperature of the reservoir, ( ) is the amplitude of the annual variations in the reservoir temperature, is the frequency of temperature variations, and the parameters and 0 are illustrated in Figure 4. The mean annual temperature of the reservoir, ( ), can be obtained as follows [5]: where is the mean annual temperature of the reservoir at the surface, is the mean annual temperature of the reservoir at the bottom, is the reservoir depth, and the parameter is 0.04. The detailed formulations on how and are calculated can be found in Bofang [5]. Moreover, the amplitude of the annual variations in the reservoir temperature, ( ), is Parameter 0 is a function of reservoir height and is calculated as: The minimum allowed value for 0 is 0.018. Finally the parameter may be calculated by the following equation: where the parameters , , and are 0.085, 2.15, and 1.30, respectively.

Solar Radiation
Effect. Solar radiation, according its nature, leads to increasing of temperature at the exposed surfaces of the dam body and is a function of the slope, direction, and latitude of the surface with respect to sun. The navigation for all parts of the gravity dams is similar and only one value for each side of the dam body would be required, while this value should be calculated for several locations of an arch dam with respect to its shape. The amount of the solar radiation received by an arch dam depends on a series of periodic seasonal changes. This variation is a function of different factors such as the height of the dam site above the sea level, surface directions relative to sun, surface slope relative to horizon, the region cloud cover, surrounding topography of the dam site, and the time of the year. According to the above definitions, the total amount of radiation absorbed by the surface can be calculated from (4). In order to have a clear vision of the total solar energy received on an inclined surface, , Table 1 and Figure 5 summarize various required angles in the solar radiation problem. Consequently, total solar energy received on an inclined surface can be calculated as [16] where is total days that solar is radiated on horizontal surface on the ground, is daily diffused solar radiation The Scientific World Journal 5 Table 1: Summary of the important angles in solar radiation problems.
Symbol Name Definition Latitude The angular position of the given surface to the equator Declination The angular position of the sun at the solar noon to the equator Slop The angle between the given surface to the horizon Surface azimuth angle Illustration diversion of surface normal vector on the horizontal plane to the south latitude

Hour angle
The angular distance of the sun's east and west to origin meridian due to earth rotation in rate of 15 degrees per hour around its axis Angle of incident The angle between the direct radiation on the surface and its normal vector Zenith angle The angle between the vertical line and the line along the sun Solar height angle The angle between sun and horizon Solar azimuth angle The angular distance between the geographic southern and the sun illustration on horizon on horizontal surface on the ground, and is the diffuse reflectance coefficient of the surrounding surface, generally ground surface, and reservoir in issue of arch dams and is given in Table 2 for various surrounding conditions. Parameter is a geometric factor and is defined as the ratio of direct radiation on the inclined surface to the horizontal surface and is derived as where and are calculated as follows: Values of , , and are calculated using the following relations: The total days that solar is radiated, , can be considered as the sum of direct and indirect radiations. Because of different orientation and slope of various locations on the exposed surfaces of arch dams, it is necessary to evaluate the direct and indirect solar radiation. Thus, the Angstrom equation ((20)-2nd part) is used to get the amount of in place sunny hours. Consider where and are experimental constants depending on geographical location of the dam. Their values for Dez Dam, 6 The Scientific World Journal as case study in the present paper, are 0.2515 and 0.446, respectively. Continuously, parameter is mean monthly sunny hour in station, is the sunny day length for each month that can be calculated from (21), and 0 is monthly extraterrestrial radiation on horizontal surface which can be calculated using (22): where sc is solar constant and is equal to 1353 Wm −2 and is number of the day in AD year. Parameter 0 is calculated using (22) on a daily basis while if the mean monthly of 0 is intended, it should be extracted from Table 3. Knowing the values of and 0 , can be calculated using (23) and finally we have the total radiation on the inclined surface:

Finite Element Implementation of Heat Transfer
In the present study, finite element technique is used for numerical implementation of the heat transfer problem within the dam body. The mathematical description of heat transfer relation for a system including anisotropic material in 3D space was explained in (1). Using the appropriate boundary conditions, as discussed in previous sections, it is possible to define the general matrices for heat transfer as [15] [ where { ( ) } is temperature vector, [ ( ) ] is element conspecific heat matrix, [ ( ) ] is element conductivity matrix, { ( ) } is element heat generation vector, and { ( ) } is element convection surface heat flow vector and all are defined as follows: , where , , and are unit normal vectors on the considered surface and , , and are defined as where , , and are conduction coefficients of mass concrete in , , and global directions, respectively, and, finally, is the shape function for the th node in the considered thermal element.

Arch Dam Description. Dez
Dam is located at 20 km NE of Andimeshk, Khuzestan, Iran (Figure 6(a)). This dam is considered as a thin double curvature arch dam. The dam height is 203 m from the foundation and 190 m from the river bottom. Its thickness is 4.5 m at the crest and 21 m at the foundation base level of the crown cantilever [18]. The dam crest length is 212 m and its crest level is 354 asl (above sea level). The maximum operation level of the dam was designed at an altitude of 350 asl. In order to optimize the dam and its hydropower operation, it has been increased to the level of 352 asl [18].
A concrete Pulvino was provided between the concrete body and the foundation rock for suitable load transfer from the dam body to its abutments. There is an overall joint between the dam body blocks and the Pulvino called peripheral joint. Figure 6

Thermal Finite Element Model.
A finite element model of the dam body was provided for thermal transient analysis. The provided model includes 1748 isoparametric linear elements in three layers through the thickness connected through 2120 nodes (Figure 7(a)). Mostly, eight-node solid thermal elements with one thermal degree of freedom were used for modeling while some prism elements were used in vicinity of Pulvino (Figure 7(b)).

Applied Initial Temperature.
In order to perform a valid thermal analysis, it is required to know the initial temperature (reference temperature) at the system. Considering that Dez Dam is about 52 years old with limited instrumentation within the dam body, there is no reliable information about initial temperature distribution within the body [18]. So a heuristic method is used in the current study for determination of the initial temperature. In the first step, a constant temperature of 23 ∘ C is applied to the all nodes which is defined by tracking the similar dams and the engineering judgment. Then, thermal transient analysis is performed for one year considering the imposed boundary conditions. Results of the thermal distribution on the nodes at the end of this step are applied to the model again as the new initial conditions for the next step and this loop is repeated for several years until approaching the stable responses. The results of the thermal distribution at the end of the last step can be interpreted as real initial temperature of the model after several years. Figure 8 shows the proposed algorithm for calculating the initial temperature distribution on the dam body. Figure 9 compares the results obtained from numerical analysis with the data collected from thermometers installed in block number 9 corresponding to year 1973. Results indicated that the responses from thermal analysis are in good compliance with the actual values. It should be mentioned that there is no middle point in numerical model and the results for this point are the average of two side nodes and some differences between the actual and theoretical data are expected.

Air Temperature.
As mentioned before both the recorded data and proposed sinusoidal function in (8) could be used for determination of the air temperature at the dam site. Considering that the nearest meteorological station to the dam site is station number 21-144, so these results are used as reference. Air temperature is recorded monthly in this station, so it is required to use the linear interpolation in order to extract daily air temperature. Figure 10 shows the sample of air temperature variation for the period of six years in this station.

Reservoir Water Temperature. The empirical-analytical
Bofang's relations are used to calculate the temperature of the reservoir water [5]. Using (11)-(15) and considering = 2.15 and = 1.30, it is possible to simulate the water temperature at different elevations. A sample of computational results in comparison with values recorded by the water temperature thermometers at the upstream face in block number 9 is shown in Figure 11. According to these figures, the calculated values show good agreement with the real values. The differences between analysis and thermometer data for level 328 m asl is due to decreasing reservoir water in some months to lower than this level and so this part of the dam is in expose of the air temperature instead of the water.
The Scientific World Journal

Effect of Solar Radiation.
In order to implement the solar radiation, the direct and indirect values of daily radiation on horizontal plane should be determined. Due to the lack of the required instruments for measuring radiation near the dam site, its amount is determined according to the recorded daily sunshine hours at the Dez station for several years and calculated using (12)- (16). Direct, indirect, and total radiations calculated for the six years are shown in Figure 12.
The reflection coefficient is assumed to be 0.2 for the region surrounding the site and 0.3 for water surface. Estimated coefficient of the surface reflectance values and variation of reflection coefficient from water surface are shown in Table 3 and Figure 13, respectively. Also the absorption coefficient for the body concrete is considered to be 0.7.

Reservoir Water Level Changes.
In order to properly investigate the environmental temperature effects on the thermal distribution of the dam body, the reservoir water level was tracked carefully corresponding to reservoir operation regime. To do so, the upper part of the upstream face, which is not in contact with water, is detected based on reservoir daily water level. At each step of thermal transient analysis, the lower part of the upstream face, immersed in water, is affected only by reservoir temperature while the part out of the water is affected by the air temperature and solar radiation. Figure 14 shows the variation in water level for the period of six years.
In the thermal analysis process, the first step is to determine daily thermal parameters such as air temperature surrounding the dam, water reservoir temperature, water level variations, and radiation energy. The second step is implementation of these values in the thermal transient analysis. In the current study, thermal transient analysis is performed for a time period of six years and the air and water temperatures are applied to the model as a convective flow while the solar radiation is applied to respective nodes as heat flow. Also, the foundation temperature was considered to be constant. The overview of applied thermal boundary conditions to the finite element model is shown in Figure 15.

Numerical Results versus Recorded Data
6.1. Thermal Calibration of Model. The model calibration was performed comparing the calculated data with the actual values recorded by thermometers at block number 9, which is the only instrumented block in the Dez Dam. The calibration process was conducted for six years, 1972 to 1977, because only the data for these years were available. To do so, the initial thermal condition was determined for January 1972 and thermal transient analysis was performed using a daily time step. In other words, there were 2190 load steps in thermal transient analysis. Figure 16 shows the location of block number 9 in the dam and also the position of the thermometers at various levels and sections. Table 4 represents the initial and final values of important parameters which the calibration of the finite element model is based on. All these parameters were obtained after almost wide range sensitivity analysis on the responses. Finally, Figure 17 shows the comparison of the results of calibrated model at various levels and sections with those recorded by thermometers for considered time period at reference block. Results confirmed that the provided thermal model considering the solar radiation effects is able to model the real conditions and, therefore, can be used for predicting the thermal behavior of the dam body.

Thermal Analysis considering Solar Radiation Effects.
According to Figure 18, solar radiation has an undeniable role at creation of the nonuniform temperature distribution on the dam surfaces. It means that the temperature differences between some coincident nodes highly affected by solar radiation are about 6 ∘ C to 10 ∘ C for the most hot and cold month, respectively. Also, the reservoir water fluctuation is effective in temperature controlling of the parts under water level on the upstream face. In the arch dams, their own especial geometry leads to heat concentration at some areas even at the same level. Due to location of Dez Dam which is in northern of the earth and considering that the axis of   0  5  10  15  20  25  30  35  40  45   0  5  10  15  20  25  30  35  40  45   0  5  10  15  20  25  30  35  40  45   0  5  10  15  20  25  30  35  40  45   0  5  10  15  20  25  30  35  40  45   0  5  10  15  20  25  30  35  40  45   0  5  10  15  20  25  30  35  40

12
The Scientific World Journal Jan. Feb.

May
Jun.

May
Jun.
Nov. Dec. the dam has a slight angle (about 6 ∘ ) to the geographic south, the sun which is located in the southern glows to downstream face and leads to more heat concentration near the abutments and middle elevations. The effect of direct solar radiation is almost negligible on the upstream face because it is fully behind the sun. This distribution is affected by indirect radiation reflected from the environmental surfaces, which is almost insignificant. The contours depicted in Figure 18 also show the effects of reservoir water level and temperature on thermal distribution of the upstream face. Temperature variation is not considerable in the lower part of the dam, especially along the same level, because the concrete in these regions is affected mostly by the reservoir temperature instead of the air temperature and solar radiation. Figure 19 shows the temperature distribution in the dam body neglecting the solar radiation effect. It can be clearly seen that the solar radiation had a lot of effects on temperature distribution. Neglecting solar radiation effects at the downstream face leads to completely uniform temperature distribution because only the air temperature acts as the thermal boundary condition. On the other hand the same trend is shown on the upstream face in which temperature distribution changes only by changing reservoir levels and there is no nonuniform temperature at the nodes located at the same level.

Conclusion
Thermal transient analysis of a double curvature high thin arch dam was studied in the present paper considering solar radiation effect. Air temperature, water temperature in various levels of the reservoir, and also its fluctuations were considered in the present study. Solar radiation effects were implemented in finite element model considering factors such as the height of the dam site above sea level, surface directions relative to the sun, surface slope relative to horizon, the region cloud cover, surrounding topography of the site, and the time of the year. Also the effects of direct and indirect solar energy absorption were investigated.
It is obvious that the temperature distribution on the upstream in two cases is equal except where the water level is below the thermometer's level. In this case, the temperature distribution on the upstream surface is changed due to the effect of solar radiation. According to the results of thermal analysis, the following can be concluded.
(i) The used boundary conditions provide realistic responses from the numerical models according to the results obtained from calibration of the finite element model.
(ii) The estimated reservoir water temperature from empirical-analytical Bofang's relation shows a good agreement with recorded values.
(iii) Comparison of temperature distribution on the upstream face represents that the parts lower the water level are mostly affected by reservoir temperature and its fluctuations and have uniform temperature distribution along the same levels, while the part above the water level is affected mainly by solar radiation and air temperature which leads to nonuniform thermal distribution along the various levels.
(iv) Comparing the temperature distribution on the downstream face shows that it is affected mostly by solar radiation and air temperature, which leads to nonuniform distribution of temperature along different levels. Due to double curvature shape of the considered dam, the areas in vicinity of the abutments have more variation in temperature than the central parts.
(v) Radiation effect reveals up to 10 ∘ C temperature difference between coincident nodes especially on the downstream face when the effects of sunlight are considered.
(vi) Comparison of temperature distribution which resulted from finite element method and recorded data through the thickness of the dam at different levels shows that the analytical results can estimate results of thin arch dams with good accuracy.