**Study of Wind-Induced Interference Effects on the Fujian Earth-Buildings**

Peng Xingqian, Liu Chunyan and Chen Yanhong *College of Civil Engineering, Huaqiao University, Quanzhou China* 

#### **1. Introduction**

As the only large-scale rammed-earth dwelling worldwide,Fujian earth-building gets much attention for its unique style, grand scale,ingenious structure, abundant cultural connotation, reasonable layout and the concept of keeping harmony with nature. In July 2008, Chuxi earthbuilding cluster, Hongkeng earth-building cluster, Gaobei earth-building cluster, Yangxiang Lou and Zhenfu Lou were listed among world heritage. They are important parts of Fujian earth-building with a long history, vast distribution, various types and rich connotation. Earthbuilding culture roots in oriental ethical relations and provides specific historical witness to traditional style of living by clansman, It is a unique achievement by employing rammed raw earth in large scale with "outstanding universal value".

Because of the high frequency of typhoon between summer and autumn in mountainous areas of the western Fujian, buildings in high and open areas often get serious damages, as shown in figure 1. In 2006, the 4th cyclone "BiLiSi" brought heavy damage to Daoyun Lou which is 400 years old. Six rooms in it collapsed, several tiles were blew off and the total number of damaged rooms reached more than 10. As one of the world cultural heritages, the protection, utilization and development of Fujian earth-building is the major issue to be deal with. Presently, the theory study for wind-resistant of low buildings is still not enough, the failure mechanism hasn't been studied thoroughly. For low buildings often appear in the form of groups, the related studies are even less. So research of wind interference effect in earth-building groups can not only fill the blank of research studies but also put forward some corresponding measures for protection of the world cultural heritage.

Fig. 1. Storm damage to the roof of earth-buildings

Study of Wind-Induced Interference Effects on the Fujian Earth-Buildings 5

interference effects between groups, and ignoring the interference of the building which less

In previous tests of wind interference, we remain in the interference effect between two buildings for a long time, and rarely consider the interference effects between more than three buildings. Professor Xiezhuangning[6] studied the wind-induced disturbance response between the three buildings, and analysis the interference of the characteristics and mechanism by neural networks, spectral analysis and statistics. The results show that the combined effects of the two buildings would be stronger than a single building. Under the landscape of Class B, Interference factor of two buildings would be increased more than 79%

Holmes[7] study the wind characteristics of the street on both sides of the building, and found upstream of the shadowing effect and the building construction the distance between a great relationship. Zhao qingchun[8] have studied the low gable roof wind-induced interference effect, found group effect on the windward roof pressure front and rear degree of influence. When the workshop's distance was 2b the interference obviously. The wind tunnel experiments show that: when buildings adjacent cross-wind side by side, the gap flow effect presence in the region when S / D ≤ 2,(S for the building spacing, D for the side of building);

Tsutsumi,J.[9]conducted a model test in different wind direction of wind load characteristics of the group, received the average wind pressure coefficient of the windward and leeward of buildings' surface. Compare and analyze the model's average wind pressure coefficient under different architectural layout, Get the average wind pressure coefficient varies with the change of wind direction. Generally speaking, the flow separation zone will increased by the skew

Fujian earth-building is ring-shaped with one-ring building or more. Here we simplify the model with ignoring the hallway, ancestral temple and such subsidiary structures. In the study of wind-induced interference, we select the typical circular earth-building to do the numerical simulation. The diameter of the biggest circular earth-building is 82 meters and the smallest is 17 meters while the common number of stories is 2 ~ 4 [10]. The research model this article selects is 28 meters in diameter, 3 layers, 11.2 meters high and under conditions as shown in fig. 3. This chapter mainly studies the characteristics of wind load and the air flowing field of circular earth-buildings,the change rule under different spacing of which are also explored. Here we major change the windward spacing between two earth-buildings and wind direction, as shown in fig.03. S stands for spacing between two earth-buildings, D stands for horizontal scale (the diameter of the bigger circular earth-building), n is valued respectively by 0.15, 0.5, 0.75, 0.25, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5 and 4.0. The characteristics of Fujian earth-buildings are: clay wall, general 1-m ~ 1.5 m of wall thick, chines-style tile roof and big pick eaves. When

When Buildings are along the windward, shielding effect exists in S / D ≤ 3 regions.

wind, air disturbance will more severe, and the flow will become more complex.

**3. Analysis of wind interference effect between two Earth-buildings** 

than half the height of buildings[5].

of a single building.

**2.3 The influence of number of buildings** 

**2.4 The influence of spacing of building group** 

**2.5 The impact of the wind stream** 

**3.1 Calculation model** 

## **2. The influencing factors of wind interference effect**

Fujian earth-buildings are often located in the form of groups, as shown in figure 2. Surface wind load is heavily influenced by the surrounding buildings and the main influencing factors include the height of the building, the relative position between buildings, section size and shape, the wind speed and wind direction, the type of wind field, etc.

Fig. 2. Tianluokeng earth-building cluster

#### **2.1 The influence of landscape**

Roughness of the landscape has a great influence on the structure wind loads, And under different wind, the interference effects between the buildings are quite different from the wind. Compared to the isolated building at the open area, Walker and Roy[1] found that the average load, peak load and bending moment are increased in the urban area of wind load. Under the open countryside and suburban areas of different topography, Case. P.C [2] study on the transient external pressure of the gable roof building experimental. He pointed out that wind load at buildings in the city suburbs is lower than that in the open landscape. And the arrangement of groups help reduced the load on a single. Blessmann[3] studied variety of landscape effects of wind interference, The results show that the moderating effect of the open landscape is most evident. Because of the turbulence is relatively low in open landscape, The pulse of wake in the upstream building has a strong correlation, Therefore, wind loads on downstream buildings caused by increased.

#### **2.2 The influence of building's width and height**

The width of the windward side of the housing has great influence on eddy size behind the leeward side, And the size of the upstream building construction also affect the downstream response of wind interference. Taniike[4] study the Wind-induced interference effects under low turbulence contour and different section size in square columns, He pointed out that the average wind load to the along wind will decline with the increase size of upper building's section, and that dynamic response to the along wind will increase with increasing section width. Under normal circumstances, when the height of adjacent buildings is equal to or greater more than half of the height of the building, we should take into account the mutual interference effects between groups, and ignoring the interference of the building which less than half the height of buildings[5].

#### **2.3 The influence of number of buildings**

4 Fluid Dynamics, Computational Modeling and Applications

Fujian earth-buildings are often located in the form of groups, as shown in figure 2. Surface wind load is heavily influenced by the surrounding buildings and the main influencing factors include the height of the building, the relative position between buildings, section

Roughness of the landscape has a great influence on the structure wind loads, And under different wind, the interference effects between the buildings are quite different from the wind. Compared to the isolated building at the open area, Walker and Roy[1] found that the average load, peak load and bending moment are increased in the urban area of wind load. Under the open countryside and suburban areas of different topography, Case. P.C [2] study on the transient external pressure of the gable roof building experimental. He pointed out that wind load at buildings in the city suburbs is lower than that in the open landscape. And the arrangement of groups help reduced the load on a single. Blessmann[3] studied variety of landscape effects of wind interference, The results show that the moderating effect of the open landscape is most evident. Because of the turbulence is relatively low in open landscape, The pulse of wake in the upstream building has a strong correlation, Therefore,

The width of the windward side of the housing has great influence on eddy size behind the leeward side, And the size of the upstream building construction also affect the downstream response of wind interference. Taniike[4] study the Wind-induced interference effects under low turbulence contour and different section size in square columns, He pointed out that the average wind load to the along wind will decline with the increase size of upper building's section, and that dynamic response to the along wind will increase with increasing section width. Under normal circumstances, when the height of adjacent buildings is equal to or greater more than half of the height of the building, we should take into account the mutual

size and shape, the wind speed and wind direction, the type of wind field, etc.

**2. The influencing factors of wind interference effect** 

Fig. 2. Tianluokeng earth-building cluster

wind loads on downstream buildings caused by increased.

**2.2 The influence of building's width and height** 

**2.1 The influence of landscape** 

In previous tests of wind interference, we remain in the interference effect between two buildings for a long time, and rarely consider the interference effects between more than three buildings. Professor Xiezhuangning[6] studied the wind-induced disturbance response between the three buildings, and analysis the interference of the characteristics and mechanism by neural networks, spectral analysis and statistics. The results show that the combined effects of the two buildings would be stronger than a single building. Under the landscape of Class B, Interference factor of two buildings would be increased more than 79% of a single building.

#### **2.4 The influence of spacing of building group**

Holmes[7] study the wind characteristics of the street on both sides of the building, and found upstream of the shadowing effect and the building construction the distance between a great relationship. Zhao qingchun[8] have studied the low gable roof wind-induced interference effect, found group effect on the windward roof pressure front and rear degree of influence. When the workshop's distance was 2b the interference obviously. The wind tunnel experiments show that: when buildings adjacent cross-wind side by side, the gap flow effect presence in the region when S / D ≤ 2,(S for the building spacing, D for the side of building); When Buildings are along the windward, shielding effect exists in S / D ≤ 3 regions.

#### **2.5 The impact of the wind stream**

Tsutsumi,J.[9]conducted a model test in different wind direction of wind load characteristics of the group, received the average wind pressure coefficient of the windward and leeward of buildings' surface. Compare and analyze the model's average wind pressure coefficient under different architectural layout, Get the average wind pressure coefficient varies with the change of wind direction. Generally speaking, the flow separation zone will increased by the skew wind, air disturbance will more severe, and the flow will become more complex.

#### **3. Analysis of wind interference effect between two Earth-buildings**

#### **3.1 Calculation model**

Fujian earth-building is ring-shaped with one-ring building or more. Here we simplify the model with ignoring the hallway, ancestral temple and such subsidiary structures. In the study of wind-induced interference, we select the typical circular earth-building to do the numerical simulation. The diameter of the biggest circular earth-building is 82 meters and the smallest is 17 meters while the common number of stories is 2 ~ 4 [10]. The research model this article selects is 28 meters in diameter, 3 layers, 11.2 meters high and under conditions as shown in fig. 3. This chapter mainly studies the characteristics of wind load and the air flowing field of circular earth-buildings,the change rule under different spacing of which are also explored.

Here we major change the windward spacing between two earth-buildings and wind direction, as shown in fig.03. S stands for spacing between two earth-buildings, D stands for horizontal scale (the diameter of the bigger circular earth-building), n is valued respectively by 0.15, 0.5, 0.75, 0.25, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5 and 4.0. The characteristics of Fujian earth-buildings are: clay wall, general 1-m ~ 1.5 m of wall thick, chines-style tile roof and big pick eaves. When

Study of Wind-Induced Interference Effects on the Fujian Earth-Buildings 7

while domain 2 adopts convergence higher structured hexahedral meshes. Fujian Earthbuilding is located in rural mountain areas , h0 =10m,v0 =5.35m/s, Fujian Earth-building

I(z)=0.194,turbulence integral scale Lu =60.55m,kinetic energy k(z) and dissipation rate

4 () ( )

(a) Meshing of internal domain 1 (b) Meshing of internal domain 2

This paper adopt every 45° wind direction to do wind interference analysis, for symmetry of the structure, three conditions were simulated in this paper under the same spacing. This paper analyzed wind pressure coefficient at local wind vector at Earth-building 2/3 highly level profile and centre vertical profile, and contrasted wind pressure coefficient between

At the windward area, flow has a positive stagnation point at 2/3 highly level profile, from the stagnation point airflow radiate outward[9]. At the area above the point, the current rise upward and beyond Earth-building roof top; at the area below the point, airflow downward and flow to the ground. So this paper choose 2/3 highly level profile to discuss the wind field characteristics. Meanwhile, this paper select of center vertical profile as features

*z*

 

*kz Iz uz C kz*

3/4 3/2 ( ) 0.5 [ ( ) ( )]

*u*

*KL* 

The surface of buildings use non-slip wall, the two sides and top surface of the numerical wind tunnel use free gliding wall, the outlet of the numerical wind tunnel use open pressure export. This paper argues turbulence is fully development (The static pressure is zero).

model).

2

=0.16. Turbulence intensity

(1)

area belongs to the class B landform, roughness index

Turbulence model adopt shear stress transport model (SST *k*

ε(z) are adopted as the following form:

Fig. 5. Meshing of domain

**3.2.1 0° wind direction** 

**3.2 Analysis of wind characteristic** 

monomer Earth-building and group Earth-building

typhoon comes, the big pick eaves are most easily swept away, leading to the damage of whole roof structure. Therefore, the roof zoning plan of the research object (i.e. the disturbed body) is shown in figure 4. The dividing is in a counterclockwise direction and the roof is divided into eaves part and ridge part. The upper surface of outside carry eaves are signed respectively by WTS1 ~ WTS8, the lower surface by WTX1 ~ WTX8. The upper surface of inside carry eaves are signed respectively by NTS1~ NTS8, the lower surface by NTX1 ~ NTX8. The ridge part is divided into the inside part and outside part, and they are signed respectively by NJ1~ NJ8 and WJ1~ WJ8.

Fig. 3. Plan of Earth-building and wind direction

#### Fig. 4. Roofing zoning

Settings of basic parameters in numerical simulation: according to reference literatures[11], domains' size can be set as following: BLH=600m500m100m, its blocking rate is 0.6%, meets requirements. As shown in figure 5, the whole calculation domain is divided into two parts: internal area and external area. The cylinder with 380m diameter is the internal areadomain 1, the other part is external area-domain2. Domain 1 use the tetrahedron meshes while domain 2 adopts convergence higher structured hexahedral meshes. Fujian Earthbuilding is located in rural mountain areas , h0 =10m,v0 =5.35m/s, Fujian Earth-building area belongs to the class B landform, roughness index =0.16. Turbulence intensity I(z)=0.194,turbulence integral scale Lu =60.55m,kinetic energy k(z) and dissipation rate ε(z) are adopted as the following form:

$$\begin{cases} k(z) = 0.5 \times [I(z) \times \overline{u}(z)]^2 \\ \varepsilon(z) = \frac{4C\_{\mu}^{3/4} k(z)^{3/2}}{KL\_u} \end{cases} \tag{1}$$

The surface of buildings use non-slip wall, the two sides and top surface of the numerical wind tunnel use free gliding wall, the outlet of the numerical wind tunnel use open pressure export. This paper argues turbulence is fully development (The static pressure is zero). Turbulence model adopt shear stress transport model (SST *k* model).

(a) Meshing of internal domain 1 (b) Meshing of internal domain 2

Fig. 5. Meshing of domain

6 Fluid Dynamics, Computational Modeling and Applications

typhoon comes, the big pick eaves are most easily swept away, leading to the damage of whole roof structure. Therefore, the roof zoning plan of the research object (i.e. the disturbed body) is shown in figure 4. The dividing is in a counterclockwise direction and the roof is divided into eaves part and ridge part. The upper surface of outside carry eaves are signed respectively by WTS1 ~ WTS8, the lower surface by WTX1 ~ WTX8. The upper surface of inside carry eaves are signed respectively by NTS1~ NTS8, the lower surface by NTX1 ~ NTX8. The ridge part is divided into the inside part and outside part, and they are signed

Settings of basic parameters in numerical simulation: according to reference literatures[11], domains' size can be set as following: BLH=600m500m100m, its blocking rate is 0.6%, meets requirements. As shown in figure 5, the whole calculation domain is divided into two parts: internal area and external area. The cylinder with 380m diameter is the internal areadomain 1, the other part is external area-domain2. Domain 1 use the tetrahedron meshes

respectively by NJ1~ NJ8 and WJ1~ WJ8.

Fig. 3. Plan of Earth-building and wind direction

Fig. 4. Roofing zoning

#### **3.2 Analysis of wind characteristic**

This paper adopt every 45° wind direction to do wind interference analysis, for symmetry of the structure, three conditions were simulated in this paper under the same spacing. This paper analyzed wind pressure coefficient at local wind vector at Earth-building 2/3 highly level profile and centre vertical profile, and contrasted wind pressure coefficient between monomer Earth-building and group Earth-building

#### **3.2.1 0° wind direction**

At the windward area, flow has a positive stagnation point at 2/3 highly level profile, from the stagnation point airflow radiate outward[9]. At the area above the point, the current rise upward and beyond Earth-building roof top; at the area below the point, airflow downward and flow to the ground. So this paper choose 2/3 highly level profile to discuss the wind field characteristics. Meanwhile, this paper select of center vertical profile as features

Study of Wind-Induced Interference Effects on the Fujian Earth-Buildings 9

buildings is more complex. Air pressure coefficient isocline mutual surrounded relatively intense, and the value has reduce trend. Which is especially noteworthy is the wind pressure coefficient isocline of perturbation building is quite different to monomer at windward direction; it appears two isocline large regions. The interfered building is affected by two vortexes at the tail of the front Earth-building. When the spacing is 0.75 D, vortex is gradually developed, air pressure coefficient isocline between two Earth-buildings is linked together, and mutual interference is still evident. When the spacing is 1.5 D, the whirlpool basically develops fully and the isocline is tending to independence. When the spacing continues to increase to 3.0 D, development of whirlpool is fully, air pressure coefficient isocline around two Earth-buildings is full independence and tend to monomer conditions. At 0 ° wind direction, generally speaking, flow field of downstream Earth-building changes

Figure 7 gives wind pressure coefficients isocline of center vertical profile in different

(b) Wind pressure coefficient of center vertical

 (d) Wind pressure coefficient of center vertical profile when S = 0.75D

(f) Wind pressure coefficient of center vertical

profile when S =3D

Fig. 7. Wind pressure coefficient of central vertical profile of 0 ° wind direction

profile when S = 0.15D

greatly, downstream Earth-building under the more obvious influence.

(2) Center vertical profile of 0 °wind direction

(a) Wind pressure coefficient of center

(c) Wind pressure coefficient of center

(e) Wind pressure coefficient of center

vertical profile when S = 0.75D

vertical profile when S = 2D

vertical vertical profile

spacing.

surface, analyze flow field characteristics between Earth groups Building through the observation of wind pressure coefficient graph of this vertical profile.

(1) Level cross section at 0° wind direction

Fig. 6 shows the isocline of the air pressure coefficient at 2/3 highly level profile. From Fig 6(a) we can see the isocline wind pressure coefficient is very plump at the windward area and the two sides of single building. Wind pressure coefficient is positive in the windward,

 (a) Wind pressure coefficient of cross section of single building

(c) Wind pressure coefficient of cross section

 (e) Wind pressure coefficient of cross section when S=2.0D

 (d) Wind pressure coefficient of when S=1.5D

 (f) Wind pressure coefficient of cross when S=3.0D

Fig. 6. Wind pressure coefficient of 2/3highly level profile at 0 ° wind direction

and the closer to the building, the bigger it is. While this coefficient is negative in the side, and the closer to the building, the bigger the absolute value is. But in the leeward surface, we can see two air pressure coefficient equivalent envelope for the two vortexes formed at leeward. Figure 6 (b) is air pressure coefficient graph of two spacing is 0.15 D circular Earthbuilding. Due to the distance between the two buildings smaller, flow between the two

 (b) Wind pressure coefficient of cross when S=0.15D

buildings is more complex. Air pressure coefficient isocline mutual surrounded relatively intense, and the value has reduce trend. Which is especially noteworthy is the wind pressure coefficient isocline of perturbation building is quite different to monomer at windward direction; it appears two isocline large regions. The interfered building is affected by two vortexes at the tail of the front Earth-building. When the spacing is 0.75 D, vortex is gradually developed, air pressure coefficient isocline between two Earth-buildings is linked together, and mutual interference is still evident. When the spacing is 1.5 D, the whirlpool basically develops fully and the isocline is tending to independence. When the spacing continues to increase to 3.0 D, development of whirlpool is fully, air pressure coefficient isocline around two Earth-buildings is full independence and tend to monomer conditions. At 0 ° wind direction, generally speaking, flow field of downstream Earth-building changes greatly, downstream Earth-building under the more obvious influence.

(2) Center vertical profile of 0 °wind direction

8 Fluid Dynamics, Computational Modeling and Applications

surface, analyze flow field characteristics between Earth groups Building through the

Fig. 6 shows the isocline of the air pressure coefficient at 2/3 highly level profile. From Fig 6(a) we can see the isocline wind pressure coefficient is very plump at the windward area and the two sides of single building. Wind pressure coefficient is positive in the windward,

S=0.15D

when S=1.5D

when S=3.0D

Fig. 6. Wind pressure coefficient of 2/3highly level profile at 0 ° wind direction

and the closer to the building, the bigger it is. While this coefficient is negative in the side, and the closer to the building, the bigger the absolute value is. But in the leeward surface, we can see two air pressure coefficient equivalent envelope for the two vortexes formed at leeward. Figure 6 (b) is air pressure coefficient graph of two spacing is 0.15 D circular Earthbuilding. Due to the distance between the two buildings smaller, flow between the two

(b) Wind pressure coefficient of cross when

(d) Wind pressure coefficient of

(f) Wind pressure coefficient of cross

observation of wind pressure coefficient graph of this vertical profile.

(1) Level cross section at 0° wind direction

(a) Wind pressure coefficient of cross

(c) Wind pressure coefficient of cross section

(e) Wind pressure coefficient of cross section

section of single building

when S=0.75D

when S=2.0D

Figure 7 gives wind pressure coefficients isocline of center vertical profile in different spacing.

 (a) Wind pressure coefficient of center vertical vertical profile

vertical profile when S = 0.75D

 (e) Wind pressure coefficient of center vertical profile when S = 2D

 (f) Wind pressure coefficient of center vertical profile when S =3D

Fig. 7. Wind pressure coefficient of central vertical profile of 0 ° wind direction

 (b) Wind pressure coefficient of center vertical profile when S = 0.15D

Study of Wind-Induced Interference Effects on the Fujian Earth-Buildings 11

Figure 8 45°wind direction is given level of the wind pressure coefficients cross section in 2/3 housing height place, we can see that Earth-buildings wind field changes significantly around in different wind direction for monomer Earth-building, the situation is similar with above, here is not to say much. For two Earth-buildings speaking, when spacing is 0.15 D, oblique flow fields makes Earth-buildings both sides have larger wind speed, but it is affected behind Earth-buildings, wind pressure coefficients isocline have inter-permeation by each other, and is very strong between two Earth-buildings wind pressure isocline with monomer markedly different in two Earth-buildings adjacent area and leeward surface for the front of Earthbuildings, drafting produces whirlpool is impeded, which leading to wind pressure coefficients reduce, and most regional present negative, interference phenomenon is seriously in the leeward surface Earth-buildings also wind pressure isocline with monomer markedly different in two Earth-buildings adjacent area and leeward surface for Behind Earth-buildings,

When the spacing becomes larger between two Earth-buildings, and from spacing 0.75D to 2.0D, with drafting place whirlpool developed slowly in the front of Earth-buildings, wind pressure coefficients isocline tend to be independent, interference become weak. When spacing for 3.0 D, interference has not obvious, the flow fields around the Earth-buildings is

In figure 9, 45° wind direction are given under different spacing vertical section center air pressure coefficient isocline map, From figure 9 (a) which can be seen, air pressure coefficient value of Earth-buildings in the windward side is lesser and negative, near the Earth-buildings metopic air pressure coefficient absolute value increases, in the outside carry eaves, separated phenomenon of air pressure coefficient appeared .It's all negative value in fluctuation pick eaves. Both internal and external roof are affected by negative pressure, and the internal roof endure a bigger negative pressure. The leeward side is in negative pressure area, because of the blocking by Earth-buildings windward surface, wind pressure reduced. Fluctuation pick eaves pressure coefficients of the inside carry eaves are all negative value which offset each other. Outside carry eaves fluctuation surface wind pressure coefficient size differ not quite, the most air pressure coefficient negative value appeared in the leeward side metopic place. If both earth-buildings exist together, the mutual influence is obvious. In figure 9 (b) the spacing is 0.15D, Between two Earthbuildings regional wind pressure isocline showed great difference when monomer, between two Earth-buildings it has even been inter-permeation phenomenon, Behind Earth-buildings air pressure coefficient value in the windward side is negative and its absolute value increases of the monomer Earth-building. It shows that when two Earth-buildings interact with each other, the buildings in the downstream are in the area of architectural drafting upstream effect . The influence of its surface is opposite bigger. As spacing increase, center profile around two earth-buildings distributions of air pressure mutual interference gradually decreased, and change contour tend to monomer condition. When spacing is 3.0D, wind pressure coefficient changes contour line in center vertical of Earth-buildings is similar

**3.2.2 45° wind direction** 

similar with monomer.

with monomer condition.

(1) Level cross section of 45°wind direction

but wind pressure changes very little in lateral area.

(2) Center vertical profile of 45 ° wind direction

From figure 7 we can see that wind pressure coefficients of Earth-buildings center vertical profile are similar with one, when spacing for 3.0 D. Wind pressure coefficients isocline appear separation phenomenon is quite serious in the external roofs, where the separation point expose many isocline. Under the surface of wind pressure coefficients significantly greater than upper one, which above is negative, the other is positive in the external roofs. Wind pressure coefficients of upper and under surface is close in the external roofs, wind pressure coefficients of upper and under surface almost to zero in the internal roofs, when they are in the leeward flow fields. When both ones spacing is 0.15 D, The prevailing wind direction of wind field and leeward are significantly different, the prevailing flow fields is not affected, and drafting leeward surface whirlpool didn't develop completely Because of the stop function behind Earth-buildings. Whirlpool gradually development, Earthbuildings mutual interference slowly reduce, as spacing is increasing, finally wind field becomes into a monomer.

 (a) Wind pressure coefficient of cross section of single building

 (b) Wind pressure coefficient of cross section when S=0.15D

 (c) Wind pressure coefficient of cross when S=0.75Dsection

 (e) Wind pressure coefficient of cross section when S=2.0D

 (d) Wind pressure coefficient of cross section when S=1.5D

 (f) Wind pressure coefficient of cross section when S=2.0D

Fig. 8. Wind pressure coefficient of 2/3highly level profile at 45 ° wind direction

#### **3.2.2 45° wind direction**

10 Fluid Dynamics, Computational Modeling and Applications

From figure 7 we can see that wind pressure coefficients of Earth-buildings center vertical profile are similar with one, when spacing for 3.0 D. Wind pressure coefficients isocline appear separation phenomenon is quite serious in the external roofs, where the separation point expose many isocline. Under the surface of wind pressure coefficients significantly greater than upper one, which above is negative, the other is positive in the external roofs. Wind pressure coefficients of upper and under surface is close in the external roofs, wind pressure coefficients of upper and under surface almost to zero in the internal roofs, when they are in the leeward flow fields. When both ones spacing is 0.15 D, The prevailing wind direction of wind field and leeward are significantly different, the prevailing flow fields is not affected, and drafting leeward surface whirlpool didn't develop completely Because of the stop function behind Earth-buildings. Whirlpool gradually development, Earthbuildings mutual interference slowly reduce, as spacing is increasing, finally wind field

(b) Wind pressure coefficient of cross section

(d) Wind pressure coefficient of cross section

(f) Wind pressure coefficient of cross section

when S=0.15D

when S=1.5D

when S=2.0D

Fig. 8. Wind pressure coefficient of 2/3highly level profile at 45 ° wind direction

becomes into a monomer.

section of single building

S=0.75Dsection

when S=2.0D

(a) Wind pressure coefficient of cross

(c) Wind pressure coefficient of cross when

(e) Wind pressure coefficient of cross section

#### (1) Level cross section of 45°wind direction

Figure 8 45°wind direction is given level of the wind pressure coefficients cross section in 2/3 housing height place, we can see that Earth-buildings wind field changes significantly around in different wind direction for monomer Earth-building, the situation is similar with above, here is not to say much. For two Earth-buildings speaking, when spacing is 0.15 D, oblique flow fields makes Earth-buildings both sides have larger wind speed, but it is affected behind Earth-buildings, wind pressure coefficients isocline have inter-permeation by each other, and is very strong between two Earth-buildings wind pressure isocline with monomer markedly different in two Earth-buildings adjacent area and leeward surface for the front of Earthbuildings, drafting produces whirlpool is impeded, which leading to wind pressure coefficients reduce, and most regional present negative, interference phenomenon is seriously in the leeward surface Earth-buildings also wind pressure isocline with monomer markedly different in two Earth-buildings adjacent area and leeward surface for Behind Earth-buildings, but wind pressure changes very little in lateral area.

When the spacing becomes larger between two Earth-buildings, and from spacing 0.75D to 2.0D, with drafting place whirlpool developed slowly in the front of Earth-buildings, wind pressure coefficients isocline tend to be independent, interference become weak. When spacing for 3.0 D, interference has not obvious, the flow fields around the Earth-buildings is similar with monomer.

#### (2) Center vertical profile of 45 ° wind direction

In figure 9, 45° wind direction are given under different spacing vertical section center air pressure coefficient isocline map, From figure 9 (a) which can be seen, air pressure coefficient value of Earth-buildings in the windward side is lesser and negative, near the Earth-buildings metopic air pressure coefficient absolute value increases, in the outside carry eaves, separated phenomenon of air pressure coefficient appeared .It's all negative value in fluctuation pick eaves. Both internal and external roof are affected by negative pressure, and the internal roof endure a bigger negative pressure. The leeward side is in negative pressure area, because of the blocking by Earth-buildings windward surface, wind pressure reduced. Fluctuation pick eaves pressure coefficients of the inside carry eaves are all negative value which offset each other. Outside carry eaves fluctuation surface wind pressure coefficient size differ not quite, the most air pressure coefficient negative value appeared in the leeward side metopic place. If both earth-buildings exist together, the mutual influence is obvious. In figure 9 (b) the spacing is 0.15D, Between two Earthbuildings regional wind pressure isocline showed great difference when monomer, between two Earth-buildings it has even been inter-permeation phenomenon, Behind Earth-buildings air pressure coefficient value in the windward side is negative and its absolute value increases of the monomer Earth-building. It shows that when two Earth-buildings interact with each other, the buildings in the downstream are in the area of architectural drafting upstream effect . The influence of its surface is opposite bigger. As spacing increase, center profile around two earth-buildings distributions of air pressure mutual interference gradually decreased, and change contour tend to monomer condition. When spacing is 3.0D, wind pressure coefficient changes contour line in center vertical of Earth-buildings is similar with monomer condition.

Study of Wind-Induced Interference Effects on the Fujian Earth-Buildings 13

when S=0.15D

when S=1.5D

when S=3.0D

Fig. 10. Wind pressure coefficient of 2/3highly level profile at 45 ° wind direction

conditions are 1.578 already smaller maximum absolute value.

coefficient in the adobe negative side, the greater the closer the absolute value of the earthbuilding wall, or even 1.5.But in the leeward surface, due to the drafting place formed two swirls, we can see two wind pressure coefficient equivalent envelope. When there are two circular Earth-buildings, the spacing is 0.15 D according to figure 10 (b), air flow are the prevailing wind direction, due to shunt bypass side collision between the smaller ones, adobe air spacing interaction, air pressure coefficient negative, and absolute 2.239, at maximum achieve isocline wind pressure coefficient, and mutual surrounded relatively intense numerical more monomer adobe has the tendency of increase. Along with the increasing of the spacing distance, two 0.75 D adobe air pressure coefficient between each other 1.811, to an absolute value of interference still obvious, spacing for 1.5 D, isocline tend to independence, air pressure coefficient absolute 1.712, when spacing continue to increase to 3.0 D, two adobe air pressure coefficient isocline around almost completely independent, air pressure coefficient for 1.599, with monomer absolute adobe air pressure coefficient

(b) Wind pressure coefficient of cross-section

(d) Wind pressure coefficient of cross-section

(f) Wind pressure coefficient of cross-section

(a) Wind pressure coefficient of single

(c) Wind pressure coefficient of cros-section

(e) Wind pressure coefficient of cross-

building

when S=0.75D

section when S=2.0D

 (a) Wind pressure coefficient of center profile of single building

 (c) Wind pressure coefficient of center vertical profile when S = 0.75D

vertical profile when S = 2D

 (b) Wind pressure coefficient of center vertical vertical profile when S = 0.15D

 (d) Wind pressure coefficient of center vertical profile when S = 1.5D

 (f) T wind pressure coefficient of center vertical profile when S =3D

Fig. 9. 45° direction Angle of wind pressure coefficient vertical profile central figure

#### **3.2.3 90° direction angle**

Suppose define 0° direction angle as serial, adobe layout 45 ° direction angle is inclined column, then 90° direction angle, think adobe arrangement as coordination. Through the previous analysis, we know that , with the increasing distance, mutual interference will gradually decrease, buildings' field which surround by also tend to flow around the single Earth-building conditions.

(1) Level cross section of 90°wind direction

From Figure 10 (a),we can see that in 2/3 single adobe houses at the height of the level of cross-section of the wind pressure coefficient contour maps is same with 0° wind direction, on the windward side and side lines are full, the wind pressure coefficient is positive for integrity, and the more close from the adobe metopic walls ,the greater the pressure

(b) Wind pressure coefficient of center vertical

vertical profile when S = 0.15D

 (d) Wind pressure coefficient of center vertical profile when S = 1.5D

(f) T wind pressure coefficient of center

vertical profile when S =3D

Fig. 9. 45° direction Angle of wind pressure coefficient vertical profile central figure

Suppose define 0° direction angle as serial, adobe layout 45 ° direction angle is inclined column, then 90° direction angle, think adobe arrangement as coordination. Through the previous analysis, we know that , with the increasing distance, mutual interference will gradually decrease, buildings' field which surround by also tend to flow around the single

From Figure 10 (a),we can see that in 2/3 single adobe houses at the height of the level of cross-section of the wind pressure coefficient contour maps is same with 0° wind direction, on the windward side and side lines are full, the wind pressure coefficient is positive for integrity, and the more close from the adobe metopic walls ,the greater the pressure

(a) Wind pressure coefficient of center

(c) Wind pressure coefficient of center

(e) Wind pressure coefficient of center

(1) Level cross section of 90°wind direction

vertical profile when S = 0.75D

vertical profile when S = 2D

**3.2.3 90° direction angle** 

Earth-building conditions.

profile of single building

 (a) Wind pressure coefficient of single building

 (c) Wind pressure coefficient of cros-section when S=0.75D

 (b) Wind pressure coefficient of cross-section when S=0.15D

 (d) Wind pressure coefficient of cross-section when S=1.5D

section when S=2.0D

 (f) Wind pressure coefficient of cross-section when S=3.0D

Fig. 10. Wind pressure coefficient of 2/3highly level profile at 45 ° wind direction

coefficient in the adobe negative side, the greater the closer the absolute value of the earthbuilding wall, or even 1.5.But in the leeward surface, due to the drafting place formed two swirls, we can see two wind pressure coefficient equivalent envelope. When there are two circular Earth-buildings, the spacing is 0.15 D according to figure 10 (b), air flow are the prevailing wind direction, due to shunt bypass side collision between the smaller ones, adobe air spacing interaction, air pressure coefficient negative, and absolute 2.239, at maximum achieve isocline wind pressure coefficient, and mutual surrounded relatively intense numerical more monomer adobe has the tendency of increase. Along with the increasing of the spacing distance, two 0.75 D adobe air pressure coefficient between each other 1.811, to an absolute value of interference still obvious, spacing for 1.5 D, isocline tend to independence, air pressure coefficient absolute 1.712, when spacing continue to increase to 3.0 D, two adobe air pressure coefficient isocline around almost completely independent, air pressure coefficient for 1.599, with monomer absolute adobe air pressure coefficient conditions are 1.578 already smaller maximum absolute value.

Study of Wind-Induced Interference Effects on the Fujian Earth-Buildings 15

and pick eaves is just alike. With adobe spacing 0.75 D increases, spacing, adobe air pressure coefficient between areas surrounded by abate, and absolute phenomenon of wind pressure has reduce and decrease. With increased, when spacing distance, two Earth-buildings center 2.0 D air pressure changes around the isocline section together with monomer Earth-

This paper through interference factor to quantitative description of interference effect

C C *pI IF pA*

Figure 12 shows that under wind direction 0°, the average wind pressure coefficient interference factors of each zone of the interfered Earth-building roof is changed with the change of distance. In the Figure 12, abscissa S denotes for distance, D denotes for diameter.

> WTS1 WTS2 WTS3 WTS4 WTS5 WTS6 WTS7 WTS8

WJ1 WJ2 WJ3 WJ4 WJ5 WJ6 WJ7 WJ8 0.0 0.2 0.4 0.6 0.8 1.0 1.2

IF


outside carry eaves

inside carry eaves

IF

*pA* are separately average wind pressure coefficients after and without wind

(2)

外挑檐下表面平均风压系数干扰因子

内挑檐上表面风压系数干扰因子

 (d) Interference factors of average wind pressure coefficients on upper surface of

0.15 0.25 0.50 0.75 1.00 1.50 2.00 2.50 3.00 3.50 4.00 S/D

 (b) Interference factors of average wind preasure coefficients on under surface of

0.15 0.25 0.50 0.75 1.00 1.50 2.00 2.50 3.00 3.50 4.00 S/D

WTX1 WTX2 WTX3 WTX4 WTX5 WTX6 WTX7 WTX8

NTS1 NTS2 NTS3 NTS4 NTS5 NTS6 NTS7 NTS8

**3.3 The change rule of average wind pressure coefficient disturbances** 

building working outline similar.

adobe residences groups:

**3.3.1 0°wind direction** 

外挑檐上表面平均风压系数干扰因子

外屋脊平均风压系数干扰因子

(c) Interference factors of average wind pressure coefficients of on external roof

0.15 0.25 0.50 0.75 1.00 1.50 2.00 2.50 3.00 3.50 4.00 S/D

(a) Interference factors of average wind pressure coefficients of on upper surface

0.15 0.25 0.50 0.75 1.00 1.50 2.00 2.50 3.00 3.50 4.00 S/D

C

0.0 0.2 0.4 0.6 0.8 1.0 1.2

> 0.0 0.2 0.4 0.6 0.8 1.0 1.2

ridge

IF

outside carry eaves

IF

*<sup>p</sup>*I And C

interference.

#### (2) Center vertical profile of 90 ° wind direction

Figure 11 is a different spacing center vertical section, air pressure coefficient isocline can see from figure 11, vertical center section on either side of the air pressure coefficient about isocline obvious symmetry, airflow around side of the surface wind pressure coefficient for exterior wall lateral negative, and the farther from metopic, the small wind pressure coefficient absolute YanXia outside carry with external surface wind pressure coefficient of the measured wind pressure coefficient on the side, almost the same, pick up within the surface wind pressure coefficients were coping negative, but under the surface for the absolute value is opposite bigger. Adobe, adobe has two air around the isocline except in two adjacent area changes remarkably, adobe big changes in other areas. When spacing is 0.15 D, air pressure coefficient isocline surrounded very intense, adjacent area outside carry eaves surface wind pressure coefficient negative, and more monomer when absolute value change, adobe air next pick eaves coefficient is bigger, near the Earth-buildings absolute value change big trend, maximum achieve 2.24, metopic isocline under the changing trends

(a) Wind pressure coefficient of center vertical of single building

(c) Wind pressure coefficient of center vertical profile when S = 0.75D

 (e) Wind pressure coefficient of centerprofile when S = 2D

 (b) Wind pressure coefficient of profile of center vertical profile when S = 0.15D

 (d) Wind pressure coefficient of center vertical profile when S = 0.75D

 (f) Wind pressure coefficient of center vertical vertical profile when S =3D

and pick eaves is just alike. With adobe spacing 0.75 D increases, spacing, adobe air pressure coefficient between areas surrounded by abate, and absolute phenomenon of wind pressure has reduce and decrease. With increased, when spacing distance, two Earth-buildings center 2.0 D air pressure changes around the isocline section together with monomer Earthbuilding working outline similar.

#### **3.3 The change rule of average wind pressure coefficient disturbances**

This paper through interference factor to quantitative description of interference effect adobe residences groups:

$$IF = \frac{\mathbf{C}\_{pI}}{\mathbf{C}\_{pA}} \tag{2}$$

C *<sup>p</sup>*I And C *pA* are separately average wind pressure coefficients after and without wind interference.

#### **3.3.1 0°wind direction**

14 Fluid Dynamics, Computational Modeling and Applications

Figure 11 is a different spacing center vertical section, air pressure coefficient isocline can see from figure 11, vertical center section on either side of the air pressure coefficient about isocline obvious symmetry, airflow around side of the surface wind pressure coefficient for exterior wall lateral negative, and the farther from metopic, the small wind pressure coefficient absolute YanXia outside carry with external surface wind pressure coefficient of the measured wind pressure coefficient on the side, almost the same, pick up within the surface wind pressure coefficients were coping negative, but under the surface for the absolute value is opposite bigger. Adobe, adobe has two air around the isocline except in two adjacent area changes remarkably, adobe big changes in other areas. When spacing is 0.15 D, air pressure coefficient isocline surrounded very intense, adjacent area outside carry eaves surface wind pressure coefficient negative, and more monomer when absolute value change, adobe air next pick eaves coefficient is bigger, near the Earth-buildings absolute value change big trend, maximum achieve 2.24, metopic isocline under the changing trends

.

.

Fig. 11. Wind pressure coefficient of 90° wind direction of central vertical profile

 (b) Wind pressure coefficient of profile of center vertical profile when S = 0.15D

 (d) Wind pressure coefficient of center vertical profile when S = 0.75D

(f) Wind pressure coefficient of center vertical

vertical profile when S =3D

(2) Center vertical profile of 90 ° wind direction

(a) Wind pressure coefficient of center

(c) Wind pressure coefficient of center vertical profile when S = 0.75D

(e) Wind pressure coefficient of

centerprofile when S = 2D

vertical of single building

Figure 12 shows that under wind direction 0°, the average wind pressure coefficient interference factors of each zone of the interfered Earth-building roof is changed with the change of distance. In the Figure 12, abscissa S denotes for distance, D denotes for diameter.

(a) Interference factors of average wind pressure coefficients of on upper surface outside carry eaves

(c) Interference factors of average wind pressure coefficients of on external roof ridge

 (b) Interference factors of average wind preasure coefficients on under surface of outside carry eaves

 (d) Interference factors of average wind pressure coefficients on upper surface of inside carry eaves

Study of Wind-Induced Interference Effects on the Fujian Earth-Buildings 17

3. External roof ridge: The change regulation of the wind pressure coefficient interference factors in external roof ridge is basic same as the wind pressure coefficient interference factors in upper surface of inside carry eaves. The range ability of interference factor WJ1 and WJ8 is minimum, which is about 10%. The range ability of interference factor WJ4 and WJ5 is maximum, which is up to 70%, WJ3 and WJ6 interference factors

reduced margin around 30%, WJ2 and WJ7 are same as WJ1 and WJ8 changes. 4. Upper surface of inside carry eaves: From figure 12 (d), in addition to see surface wind pressure coefficient interference factor NTS4 and NTS5 have obvious change, other various surface marked change are small coefficient of wind pressure reduction. Judging from the numerical simulation results, surface wind pressure coefficient NTS4 and NTS5 in smaller values, on change, although magnitude of 0.01 small changes, but the wind pressure coefficient embodied in the disturbances have changed greatly. Overall, wind pressure coefficient interference factor upper surface of inside carry eaves

5. Under surface of inside carry eaves: from figure 12 (e), the each zoning of the eaves changes consistently, the wind pressure coefficient decreases, it is favorable to stand up the wind. The amplitude of NTX3 and NTX6 is relatively larger,20%,it is the premises where backflow happens inside the Earthen ring. Although it is obstructed by the windward, the wind pressure has decreased; however, airflow returns violently, the wind pressure coefficient of the premises changes more than other premises. When the distance of the Earth-building is 1.5D, the influential factors have approached to 1.0, we

6. Inner roof ridge: The wind pressure coefficient disturbances of inner roof ridge and inside carry eaves have the similar variation tendency, but the variation amplitude of inner roof ridge is larger. The variation amplitude of NJ4 and NJ5 are still the biggest, they reach the 80%, the disturbances of NJ3 and NJ6 decrease about 30%, other each

7. Bare wind pressure coefficient interference factors in outside carry eaves: For characteristic of Earth-building suction's 2.5 meter large carry eaves, we consider the up and down surface wind pressure coefficient of carry eaves respectively, then compose them, so we can obtain the bare wind pressure coefficient value that it will be used in design. Increase the interval, the disturbances in previous analysis will increase from the value that less than 1 to 1.0,in considering the bare wind pressure coefficient disturbances, the surfaces disturbances of WT4 and WT5 are both more than 1.0,the reason is that surface wind pressure absorb up and press down, it is disadvantage of structure to withstand wind. With the change of interval, the surface disturbances of WT1 and Wt8 are almost no impact, the reason is that up and down surface offset each other. The surface wind pressure of WT2 and WT7 decrease more than 70%, it is advantageous to withstand wind. The surface wind pressure of WT3 and WT6 also

8. Bare wind pressure coefficient interference factors in inside carry eaves: The surface wind pressure coefficient disturbances of NT2 and NT7 increase to a little range, it is about

does not change significantly disturbances, the maximum 16%.

surface has the small variation amplitude of about 10%.

decrease to a certain degree, it is about 40%~50%.

3.0D,it is close to 1.0.

can neglect the influences.

WTX5 at the leeward side, wind pressure coefficient changing, reducing up to 150%,when the distance increases, the pressure coefficient become negative from positive, interference factor become positive from negative, and when the spacing

(e) Interference factors of average wind pressure coefficients on under surface of inside carry eaves

(g) Interference factors of net wind pressure coefficient on outside carry eaves

 (h) Interference factors of net wind pressure coefficient on inside carry eaves

Fig. 12. The interference factors of wind pressure coefficient on each roof partition of 0° wind direction

What can be obtained by figure 12:


 (f) Interference factors of average wind pressure coefficients on inner roof ridge

0.0 0.2 0.4 0.6 0.8 1.0 1.2

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4

IF

IF

内屋脊风压系数干扰因子

 (f) Interference factors of average wind pressure coefficients on inner roof ridge

内挑檐净风压系数干扰因子

coefficient on inside carry eaves

0.15 0.25 0.50 0.75 1.00 1.50 2.00 2.50 3.00 3.50 4.00 S/D

0.15 0.25 0.50 0.75 1.00 1.50 2.00 2.50 3.00 3.50 4.00 S/D

(h) Interference factors of net wind pressure

NJ1 NJ2 NJ3 NJ4 NJ5 NJ6 NJ7 NJ8

> NT1 NT2 NT3 NT4 NT5 NT6 NT7 NT8

NTX1 NTX2 NTX3 NTX4 NTX5 NTX6 NTX7 NTX8

WT1 WT2 WT3 WT4 WT5 WT6 WT7 WT8

Fig. 12. The interference factors of wind pressure coefficient on each roof partition of 0°

1. Upper surface of outside carry eaves: to sum up , the wind pressure coefficient in upper surface of outside carry eaves is reduced compared with single Earth-building, along with the increase of distance between, this trend weakened gradually, and when the distance reached 3D, interference factor approached to 1.0, interference affect basically can be ignored. Interference factor WTS1 and WTS8 have minimum amplitude, about10%, the two surfaces are the farthest from the interfered Earthbuilding. Interference factor WTS4 and WTS5 have maximum amplitude, it reached 70 %, it has great influence with scrambling Earth-building, to the benefit of windresistant. Interference factor WTS3 and WTS6 decrease amplitude is about 30%, and interference factor WTS2 and WTS7 change amplitude is less ,is basically similar to

2. Under surface of outside carry eaves: Roofing partition is symmetrical, the wind pressure coefficient interference factors under wind direction 0°change have obvious symmetry, and wind pressure coefficient are obviously reduce compared with single Earth-building. Interference factor WTS1 and WTS8 change amplitude is less changed with distance increases, in the 1.0 external floating up and down 5%, surface wind pressure coefficient interference factor WTS2 and WTS7, WTS3 and WTS6 increase with instance increase, among them , interference factor WTS2 and WTS7 are reduce mostly 20%, interference factor WTS3 and WTS6 reduce 40%. It's worth noting that WTX4 and

内挑檐下表面风压系数干扰因子

外挑檐净风压系数干扰因子

(g) Interference factors of net wind pressure

coefficient on outside carry eaves

What can be obtained by figure 12:

WTS1 and WTS8.

0.15 0.25 0.50 0.75 1.00 1.50 2.00 2.50 3.00 3.50 4.00 S/D

(e) Interference factors of average wind pressure coefficients on under surface of

0.15 0.25 0.50 0.75 1.00 1.50 2.00 2.50 3.00 3.50 4.00 S/D

0.0 0.2 0.4 0.6 0.8 1.0 1.2

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6

wind direction

IF

inside carry eaves

IF

WTX5 at the leeward side, wind pressure coefficient changing, reducing up to 150%,when the distance increases, the pressure coefficient become negative from positive, interference factor become positive from negative, and when the spacing 3.0D,it is close to 1.0.


20%.Because of airflow in the Earth-building is backflow, the wind pressure coefficient disturbances of NT1、NT8、NT3、NT6 decrease to a little range, it is about 30%.

Overall, under 0°wind direction, the disturbed is in the upstream, then the disturbing effect is not so obvious, in addition to the individual; each partition of roof is advantageous to withstand wind. If downstream Earth-building have closer distance, the outer cornice of WT4 and WT5 are disadvantage to withstand wind, but value of its surface wind pressure coefficient is less, so it has little influence on withstand wind design. When interval of Earthbuilding reaches 3D, disturbances will approach to 1.0, the interference effect almost can be ignored.

The disturbances of outer roof ridge and outer carry eaves have the similar change rules, the disturbances variation amplitude of inner roof are all smaller than outer roof, and outer carry eaves' absolute value of wind pressure coefficient is bigger than inner carry eaves, in real life situation, destroy the roof mainly begins from lifting tile of outer cornice roof, outer carry eaves is in a very disadvantage condition, so we mainly consider the outer cornice's change rule of wind pressure coefficient disturbances below.

#### **3.3.2 45° wind direction**

Figure 13 is a 45 ° wind direction under the pressure of the partition coefficient of confounding factors change with the pitch curve.

(a) Interference factors of average wind pressure coefficients on the upper surface of outside overhangs.

(c) Interference factors of net wind pressure coefficients of onside overhangs.

 (b) Interference factors of average wind pressure coefficients on the lower surface of outside overhangs.

WTX1 WTX2 WTX3 WTX4 WTX5 WTX6 WTX7 WTX8 Study of Wind-Induced Interference Effects on the Fujian Earth-Buildings 19

1. External overhangs on the surface: WTS5 in the leeward surface, and adjacent interference Earth-building, where up to 1.5.There are increasing rapidly wind pressure and very adverse wind resistance. WTS6 interference factor is less than 1.0. WTS5 in a relatively inclined all along flow position, and wind pressure coefficient disturbances is asymmetry, because of the other Earth-building influence. There are WTS1, WTS2 WTS4 and WTS8 about 1.0, which interference is not obvious. WTS3, WTS6 and WTS7 interference factor has a small decrease. When the spacing reach to 2.5 D, interference

2. The lower surface of outer overhangs: WTX5 maximum interference factor of 1.6, WTX3 and WTX7 interference factor of Leeward surface decreases above 30%, Range of other surface pressure coefficients interference factors varies by less, remain in the

3. The interference factor of net air pressure coefficient of outside carry eaves: it reaches 2.0 or above on WT6, interference effect is serious, but the wind pressure coefficient of this surface is numerical small, 0.1 orders of magnitude, the wind resistant design does not control the surface. The interference factor of WT4 is 1.4 ,which is adverse for wind resistance. The interference factor of WT5 is negative in small spacing, and for WT5 is in the leeward surface where air pressure coefficient is small, the interference effect is beneficial although it is serious. The minimum interference factor of WT3 is 0.1, and it is getting bigger with increases of spacing (the interference factor turns to 1.0 when the spacing is 2.5D). Other interference factors are floating near 1.0; the interference

4. The interference factor of net wind pressure coefficient of inside pick eaves: the variation interference factors of each roof partition are small, interference phenomenon

Overall, under the 45 °wind direction, except the interference factors of some roof partitions are near 1.0, the variation of interference factors of WT4, WT5 and WT6 are large, which

Figure 14 is a 90° wind direction under the pressure of the partition coefficient of confounding

1. External overhangs on the surface: WTS6 interference factor reached a maximum of 2.1, WTS5 interference factor is close to 2.0, WTS4 interference factor 1.6, WTS3 interference factor is also 1.53, with control effect in the gorge, the wind pressure coefficient increases more, is not conducive to Wind. Range of other surface confounding factors varies by less, WTS2 surface disturbance factor of 1.3, WTS7 interference factor of 0.9, WTS1 and WTS8 fluctuations are around 1.0. When the distance increased to 1.5D, external overhangs are reaching the district interference factor of 1.0, as the spacing

2. The lower surface of outer overhangs: WTX5 maximum interference factor of 2.1, when wind pressure coefficient of -1.23, air flow around the earth-buildings in the region seriously affected with each other. WTX4 interference factor is also higher, at 1.4. WTX6 and WTX7 interference factor decreases above 100%, large amplitude, but at the leeward wind pressure coefficient values smaller. Range of other surface pressure

coefficients interference factors varies by less, remain in the vicinity of 1.0.

factor affect will be ignore.

phenomenon is not obvious.

should be taken into consideration of design.

increases, interference effects gradually weakened.

factors is changed with the pitch curve. What can be drawn from figure 14:

vicinity of 1.0.

is not obvious.

**3.3.3 90° wind direction** 

 (d) Interference factors of net wind pressure coefficients of inside overhangs.

Fig. 13. The fact interference ors of average wind pressure coefficient on each roof partition of 450 wind direction

What can be drawn from figure 13:


Overall, under the 45 °wind direction, except the interference factors of some roof partitions are near 1.0, the variation of interference factors of WT4, WT5 and WT6 are large, which should be taken into consideration of design.

#### **3.3.3 90° wind direction**

18 Fluid Dynamics, Computational Modeling and Applications

The disturbances of outer roof ridge and outer carry eaves have the similar change rules, the disturbances variation amplitude of inner roof are all smaller than outer roof, and outer carry eaves' absolute value of wind pressure coefficient is bigger than inner carry eaves, in real life situation, destroy the roof mainly begins from lifting tile of outer cornice roof, outer carry eaves is in a very disadvantage condition, so we mainly consider the outer cornice's

Figure 13 is a 45 ° wind direction under the pressure of the partition coefficient of

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8

> 0.6 0.7 0.8 0.9 1.0 1.1 1.2

IF

Fig. 13. The fact interference ors of average wind pressure coefficient on each roof partition

outside overhangs.

IF

外挑檐下表面平均风压系数干扰因子

内挑檐净风压系数干扰因子

(d) Interference factors of net wind pressure

coefficients of inside overhangs.

0.15 0.25 0.50 0.75 1.00 1.50 2.00 2.50 3.00 3.50 4.00 S/D

 (b) Interference factors of average wind pressure coefficients on the lower surface of

0.15 0.25 0.50 0.75 1.00 1.50 2.00 2.50 3.00 3.50 4.00 S/D

WTX1 WTX2 WTX3 WTX4 WTX5 WTX6 WTX7 WTX8

NT1 NT2 NT3 NT4 NT5 NT6 NT7 NT8

WTS1 WTS2 WTS3 WTS4 WTS5 WTS6 WTS7 WTS8

> WT1 WT2 WT3 WT4 WT5 WT6 WT7 WT8

change rule of wind pressure coefficient disturbances below.

confounding factors change with the pitch curve.

外挑檐上表面平均风压系数干扰因子

外挑檐净风压系数干扰因子

(c) Interference factors of net wind pressure

coefficients of onside overhangs.

What can be drawn from figure 13:

of 450 wind direction

0.15 0.25 0.50 0.75 1.00 1.50 2.00 2.50 3.00 3.50 4.00 S/D

(a) Interference factors of average wind pressure coefficients on the upper surface of

0.15 0.25 0.50 0.75 1.00 1.50 2.00 2.50 3.00 3.50 4.00 S/D

disturbances of NT1、NT8、NT3、NT6 decrease to a little range, it is about 30%. Overall, under 0°wind direction, the disturbed is in the upstream, then the disturbing effect is not so obvious, in addition to the individual; each partition of roof is advantageous to withstand wind. If downstream Earth-building have closer distance, the outer cornice of WT4 and WT5 are disadvantage to withstand wind, but value of its surface wind pressure coefficient is less, so it has little influence on withstand wind design. When interval of Earthbuilding reaches 3D, disturbances will approach to 1.0, the interference effect almost can be

ignored.

0.4 0.6 0.8 1.0 1.2 1.4 1.6

> -0.5 0.0 0.5 1.0 1.5 2.0 2.5

IF

outside overhangs.

IF

**3.3.2 45° wind direction** 

20%.Because of airflow in the Earth-building is backflow, the wind pressure coefficient

Figure 14 is a 90° wind direction under the pressure of the partition coefficient of confounding factors is changed with the pitch curve.

What can be drawn from figure 14:


Study of Wind-Induced Interference Effects on the Fujian Earth-Buildings 21

Through the exhaustive numerical simulation about two typical circular earth-buildings in different wind directions and different spacing, we obtained the distribution characteristics of wind pressure both on the single earth-building and the groups under 0°, 45°, 90°wind direction, the interference effects of wind flow and the change rules that interference factor of average wind pressure coefficient changing with spacing on each roof partition. From the

1. Numerical simulation can accurately simulate the air flowing field in different wind directions and wind pressure distribution on each roof partition of low buildings. Through the study of wind pressure coefficient and wind velocity contour on two feature faces, we can make qualitative analyses about the wind characteristics in buildings. Under the same wind direction, air flowing field around the earth-building in groups gets close to that of single building along with the increase of building spacing. Under different wind directions, the air flowing field is not the same when

2. Because of the characteristics of earth-buildings: thick walls, long cornices and tile roof, we mainly analyses the variation of interference factors on each roof partition. The actual situation is that damage of tile roof starts with the lifting of roof tiles from the outside carry eaves. After analysis of the variation of interference factors under different wind directions, we found the change laws of interference factors on the ridge and outside carry eaves are consistent and there is less change of interference factors on the inside roof than the outside. This paper focuses on analyzing the variation of interference factors of wind

3. Under different wind directions, the variation range of interference factors of average wind pressure coefficient on both the upper and lower surface are smaller than that of net wind pressure coefficient of the outside carry eaves. Interference factors on the upper surface and that on the lower surface mutually reinforce the effects on the

4. Compared with a single earth-building, when two earth-buildings are in a line the variation of interference factors is not obvious. When the spacing between two earthbuildings reach 3 times of the bigger radius,the interference factor is close to 1.0 and

5. when the wind direction is 45° with the line of two earth-buildings,the wind interference is very different to the situation when the wind direction is 0°. When the spacing between two earth-buildings reach 2.5 times of the bigger radius,the interference factor is close to 1.0 and the interference affect can be ignored in the wind resistant design. Under the 45° wind direction, maximum interference factor reach 2.1 on part WT6 which is unfavorable disturbance, but interference factors decrease on part

6. Under the 90° wind direction, most roof partitions are severely disturbed. It shows obvious effect of narrow and wind pressure coefficient increases greatly. The interference factors on part WT5 reach up to 2.5 which are extremely unfavorable to wind resistant and should be paid attention to. When the spacing between two earthbuildings reach 1.5 times of the bigger radius,the interference effect is too weak to be

pressure coefficient on the outside carry eaves under various conditions.

windward side and mutually reduce the effects on the leeward side.

the interference affect can be ignored basically.

considered in the wind resistant design.

WT5 which is positive.

**4. Conclusion** 

research above we can draw the conclusion that:

building spacing is different.

(a) Interference factors of average wind pressure coefficients on the upper surface of outside overhangs

(c) Interference factors of net wind pressure coefficients of outside overhangs.

 (b) Interference factors of average wind pressure coefficients on the lower surface of outside overhangs

 (d) Interference factors of net wind pressure coefficients of inside overhangs.

Fig. 14. The interference factors of average wind pressure coefficient on each roof partition of 900 wind direction


Overall, in 90 ° wind direction, the interference is obvious, considering the surface pressure coefficients, most of the partition surface disturbance factor greater than 1.0, but considering the superimposed effect of the upper and lower surfaces, the interference is not so prominent, but WT5, NT4 disadvantaged status of wind, earth-building roof to the attention of conservation measures.

## **4. Conclusion**

20 Fluid Dynamics, Computational Modeling and Applications


0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4

IF

Fig. 14. The interference factors of average wind pressure coefficient on each roof partition

3. Net wind pressure coefficients interference factor outside the overhangs: WT5 surface disturbance factors up to 2.5, the pressure coefficient of 0.53, interference is obvious. WT7 interference factor of -1.7, but the surface is in the leeward wake region, pressure coefficient value is small. WT4 and WT6 interference factor decreases by about 80%, favorable wind. The surface of the other confounding factors the district did not change significantly, at 1.0 fluctuate. When the spacing of 1.5D in the Earth-building. District

4. Net wind pressure coefficients interference factor in the overhangs: the overhangs at the wind pressure coefficient absolute value of the partition is generally small, between 0.2 and 0.4. NT4 interference factor 1.32, NT5 minimum interference factor of 0.64, Earthen Ring back airflow significantly. NT8 interference factor of 1.22, NT7 interference factor of 1.18, a result of disturbed earth-building wake interference earthen interference

Overall, in 90 ° wind direction, the interference is obvious, considering the surface pressure coefficients, most of the partition surface disturbance factor greater than 1.0, but considering the superimposed effect of the upper and lower surfaces, the interference is not so prominent, but WT5, NT4 disadvantaged status of wind, earth-building roof to the attention

confounding factors close to 1.0, interference effects can be ignored.

outside overhangs

IF

外挑檐下表面平均风压系数干扰因子

内挑檐净风压系数干扰因子

(d) Interference factors of net wind pressure

coefficients of inside overhangs.

0.15 0.25 0.50 0.75 1.00 1.50 2.00 2.50 3.00 3.50 4.00 S/D

 (b) Interference factors of average wind pressure coefficients on the lower surface of

0.15 0.25 0.50 0.75 1.00 1.50 2.00 2.50 3.00 3.50 4.00 S/D

WTX1 WTX2 WTX3 WTX4 WTX5 WTX6 WTX7 WTX8

NT1 NT2 NT3 NT4 NT5 NT6 NT7 NT8

WTS1 WTS2 WTS3 WTS4 WTS5 WTS6 WTS7 WTS8

WT1 WT2 WT3 WT4 WT5 WT6 WT7 WT8

外挑檐上表面平均风压系数干扰因子

外挑檐净风压系数干扰因子

(c) Interference factors of net wind pressure

coefficients of outside overhangs.

facilities was greatly changed.

of conservation measures.

of 900 wind direction

0.15 0.25 0.50 0.75 1.00 1.50 2.00 2.50 3.00 3.50 4.00 S/D

(a) Interference factors of average wind pressure coefficients on the upper surface of

0.15 0.25 0.50 0.75 1.00 1.50 2.00 2.50 3.00 3.50 4.00 S/D

0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2


IF

outside overhangs

IF

Through the exhaustive numerical simulation about two typical circular earth-buildings in different wind directions and different spacing, we obtained the distribution characteristics of wind pressure both on the single earth-building and the groups under 0°, 45°, 90°wind direction, the interference effects of wind flow and the change rules that interference factor of average wind pressure coefficient changing with spacing on each roof partition. From the research above we can draw the conclusion that:


**0**

**2**

and Rafael Reséndiz1

*México City*

*México*

**Mass–Consistent Wind Field Models: Numerical**

<sup>1</sup>*Departamento de Matemáticas, Universidad Autónoma Metropolitana Iztapalapa,*

For several meteorological problems and a large number of applications, the knowledge of the 3–D wind field over a region is required. Examples include prediction of the transport, diffusion and dispersion of air pollutants in the atmosphere (Finardi et al., 2010; Sherman, 1978), realization of wind maps for the design of different urban and general projects (Castino et al., 2003), and the effect of wind on structures and fire spreading (Potter & Butler, 2009), among others. Moreover, meteorological wind fields are also required inputs for air quality models. In practice, usually limited horizontal wind field measurements are available, and therefore the calculation of the vertical motion must be predicted or calculated. Several methods and strategies, with various levels of complexity, have been proposed to address this problem. They can be included into two general model types: prognostic models and diagnostic models. Prognostic models are complex time–dependent hydrodynamic models governing air flow, including thermal effects, density variation and turbulent interaction. While these models are "realistic", they are expensive to operate, need extensive computer facilities, and require specialized training for their operation. On the other hand, diagnostic wind models do not require the integration of the non–linear hydrodynamic equations. Instead, available interpolated data is used to generate wind fields, which satisfy some physical or dynamical constraints. For instance, to assure mass conservation, a simplified steady–state version of the continuity equation is imposed, and the resulting model is then called a mass–consistent model. A review of these models is available in Ratto et al. (1994)

We focus in a variational mass–consistent model which is based in the original formulation by Sasaki (Sasaki, 1958). This approach has been used for a variety of meteorological problems (Castino et al., 2003; Pennel, 1983; Sherman, 1978; Wang et al., 2005). Mass–consistent models are attractive because of their simplicity, and because they are easy and economical to operate. In some applications, these models outperform the more sophisticated and expensive dynamical models (Ratto et al., 1994). However, mass–consistent models have some disadvantages, because they are based on incomplete or idealized models and have difficulty

**1. Introduction**

and Ratto (1996).

**Techniques by L2–Projection Methods**

<sup>2</sup>*División de Ciencias Básicas, Universidad Juárez Autónoma de Tabasco*

L. Héctor Juárez1, María Luisa Sandoval1, Jorge López2

#### **5. References**


## **Mass–Consistent Wind Field Models: Numerical Techniques by L2–Projection Methods**

L. Héctor Juárez1, María Luisa Sandoval1, Jorge López2 and Rafael Reséndiz1 <sup>1</sup>*Departamento de Matemáticas, Universidad Autónoma Metropolitana Iztapalapa, México City* <sup>2</sup>*División de Ciencias Básicas, Universidad Juárez Autónoma de Tabasco México*

#### **1. Introduction**

22 Fluid Dynamics, Computational Modeling and Applications

[1] Walker.G; Roy.R. Wind loads on houses in an urban environment[R]. University of Roorkee, India: Asia Pacic Symposium on wind engineering, 1985. [2] Case.P.C; Isyumov.N. Wind loads on low buildings with 4:12 gable roofs in open

[3] Blessmann.J. Buffeting effects of neighboring tall buildings [J]. Journal of Journal of wind

[4] Taniike,Yoshihito. Turbulence effect on mutual interference of buildings [J]. Journal of

[5] Zhang Xiangting. Engineering wind resistance design and calculation manual [M].

[6] Xie Zhuangning. Research of Interference Effects of Wind Loads of a Cluster of Tall

[7] Holmes.J.D. Wind pressures on tropical housing[J]. Journal of Wind Engineering and

[8] Zhao Qingchun; Peng Xingqian; Zhou Xianpeng; Qiao Changgui. Numerical simulation

[9] Tsutsumi.J; Katayama.T; Nishida.M. Wind tunnel tests of wind pressure on regularly

[10] Huang Hanmin. Fujian Earth-building [M]. Beijing: Sanglian Bookstore, 2003. (in

[11] Shao Kun; Peng Xingqian; Liu Chunyan; Xu Gang. Computational Domain Setting

ZhengZhou institute of light industry. 2010, 25(4):55-58. (in Chinese)

analysis of wind interference effects on the roof of low-rise gable-roofed buildings

aligned buildings [J]. Journal of Wind Engineering and Industrial Aerodynamics.

About Numerical Wind Tunnel Simulation of Earth-building [J].Journal of

engineering and industrial aerodynamics. 1985, 18(1): 100-105.

Buildings[D]. shanghai, Tongji university , 2003. (in Chinese)

[J]. Journal of Fuzhou university.2008, 36(6): 863-867. (in Chinese)

Aerodynamics. 1998(77-78): 107-118.

Engineering Mechanics. 1991, 117(3): 443-456.

Industrial Aerodynamics. 1994,53(1-2): 105-123.

1992, 43(3): 1799-1810.

Chinese)

China architecture &building press, 1998.(in Chinese)

country and suburban exposures [J]. Journal of Wind Engineering and Industrial

**5. References** 

For several meteorological problems and a large number of applications, the knowledge of the 3–D wind field over a region is required. Examples include prediction of the transport, diffusion and dispersion of air pollutants in the atmosphere (Finardi et al., 2010; Sherman, 1978), realization of wind maps for the design of different urban and general projects (Castino et al., 2003), and the effect of wind on structures and fire spreading (Potter & Butler, 2009), among others. Moreover, meteorological wind fields are also required inputs for air quality models. In practice, usually limited horizontal wind field measurements are available, and therefore the calculation of the vertical motion must be predicted or calculated. Several methods and strategies, with various levels of complexity, have been proposed to address this problem. They can be included into two general model types: prognostic models and diagnostic models. Prognostic models are complex time–dependent hydrodynamic models governing air flow, including thermal effects, density variation and turbulent interaction. While these models are "realistic", they are expensive to operate, need extensive computer facilities, and require specialized training for their operation. On the other hand, diagnostic wind models do not require the integration of the non–linear hydrodynamic equations. Instead, available interpolated data is used to generate wind fields, which satisfy some physical or dynamical constraints. For instance, to assure mass conservation, a simplified steady–state version of the continuity equation is imposed, and the resulting model is then called a mass–consistent model. A review of these models is available in Ratto et al. (1994) and Ratto (1996).

We focus in a variational mass–consistent model which is based in the original formulation by Sasaki (Sasaki, 1958). This approach has been used for a variety of meteorological problems (Castino et al., 2003; Pennel, 1983; Sherman, 1978; Wang et al., 2005). Mass–consistent models are attractive because of their simplicity, and because they are easy and economical to operate. In some applications, these models outperform the more sophisticated and expensive dynamical models (Ratto et al., 1994). However, mass–consistent models have some disadvantages, because they are based on incomplete or idealized models and have difficulty

**2. Mathematical formulation of the problem**

Fig. 1. Bounded region Ω.

such that **u** · **n** = 0 on Γ*N*.

closed space

∑*d*

point of *J*:

Let Ω be an open, simply connected and bounded region in **R***<sup>d</sup>* (*d* = 2 or 3) with Lipschitz boundary *∂*Ω = Γ*<sup>N</sup>* ∪ Γ*D*, where Γ*<sup>N</sup>* is the part of the boundary associated to the surface terrain (topography), and Γ*<sup>D</sup>* is the rest of the boundary (artificial vertical boundaries and top boundary), as shown in Figure 1. Given an initial vector field **u**<sup>I</sup> in Ω (which can be obtained

Mass–Consistent Wind Field Models: Numerical Techniques by L2–Projection Methods 25

by interpolating atmospheric data, or by other means), our goal is to generate a solenoidal field **u** –called adjusted field– as close to **u**<sup>I</sup> as possible in a sense that will be clarified below,

We define the following vector function spaces: **L**2(Ω) = *L*2(Ω)*<sup>d</sup>* and **H**(*div*; Ω) = { **v** ∈ **L**2(Ω) : ∇ · **v** ∈ *L*2(Ω) }. Then, the adjusted wind field **u** must belong to the normed

<sup>1</sup> *vi wi* is the usual scalar product in **<sup>R</sup>***d*. We can now formulate the problem as a least squares projection problem. For this purpose, we define a convex quadratic functional *J* : **V** → **R** as

Due to the properties of this functional, **u** ∈ **V** is a minimizer of *J* if and only if it is a stationary

*<sup>S</sup>*(**<sup>u</sup>** <sup>−</sup> **<sup>u</sup>**<sup>I</sup>

 Ω

The Lax–Milgram theorem guaranties that this equation has a unique solution.

*<sup>S</sup>*,Ω<sup>=</sup> <sup>1</sup> 2 Ω

with the norm �·�*S*,<sup>Ω</sup> associated to the inner product �**u**, **<sup>v</sup>**�*<sup>S</sup>* <sup>=</sup>

<sup>2</sup> � **<sup>v</sup>** <sup>−</sup> **<sup>u</sup>**<sup>I</sup> �<sup>2</sup>

*<sup>J</sup>*(**v**) = <sup>1</sup>

Then, our problem can be stated as follows:

*∂*

*∂� <sup>J</sup>*(**<sup>u</sup>** <sup>+</sup> *�* **<sup>v</sup>**)|*�*=<sup>0</sup> <sup>=</sup>

**V** = { **v** ∈ **H**(*div*; Ω) : ∇ · **v** = 0 and **v** · **n** = 0 on Γ*<sup>N</sup>* }, (3)

*<sup>S</sup>*(**<sup>v</sup>** <sup>−</sup> **<sup>u</sup>**<sup>I</sup>

Given **<sup>u</sup>**<sup>I</sup> <sup>∈</sup> **<sup>L</sup>**2(Ω), find **<sup>u</sup>** <sup>∈</sup> **<sup>V</sup>** such that *<sup>J</sup>*(**u**) <sup>≤</sup> *<sup>J</sup>*(**v**), <sup>∀</sup> **<sup>v</sup>** <sup>∈</sup> **<sup>V</sup>**. (5)

) · (**<sup>v</sup>** <sup>−</sup> **<sup>u</sup>**<sup>I</sup>

) · **v** *d***x** = 0, ∀ **v** ∈ **V**. (6)

<sup>Ω</sup> *S***u** · **v** *d***x**, where **v** · **w** =

) *d***x**. (4)

representing flows accurately in data–sparse regions as mountains or oceans. Despite these limitations, mass–consistent models are a valuable tool for air quality applications and consequently several developments have taken place over last decades (Ferragut et al., 2010; Ratto, 1996; Ratto et al., 1994; Ross et al., 1988; Wang et al., 2005). Most of the results presented in this chapter has been published in the last few years (Flores et al., 2010; Núñez et al., 2007; 2006), but we also include some additional ideas and recent results.

The variational method proposed by Sasaki uses the continuity equation ∇ · **u** = 0, where **u** is the wind velocity vector field on a given domain Ω. The method is based on the minimization of the functional *L* defined by

$$L\left(\mathbf{u},\lambda\right) = \frac{1}{2} \int\_{\Omega} \left\{ S\left(\mathbf{u} - \mathbf{u}^{\mathrm{I}}\right) \cdot \left(\mathbf{u} - \mathbf{u}^{\mathrm{I}}\right) + \lambda \left[\nabla \cdot \mathbf{u}\right] \right\} \, dV \, \tag{1}$$

where **u**<sup>I</sup> is an initial observed wind field, *λ* is a Lagrange multiplier and *S* is a diagonal matrix with weighting parameters *α<sup>i</sup>* > 0, *i* = 1, 2, 3, called Gaussian precision moduli, related to the scales of the respective components of the velocity field. The vertical component of the initial wind field is taken as zero because meteorological stations usually do not measure this component. The Euler–Lagrange equations of (1) are:

$$\mathbf{u} = \mathbf{u}^{\mathrm{l}} + \mathbf{S}^{-1} \nabla \lambda \,\tag{2}$$

Usually **u** is obtained from (2), after *λ* is computed. Since ∇ · **u** = 0, then from (2) we obtain the elliptic equation −∇ · *<sup>S</sup>*−1∇*<sup>λ</sup>* <sup>=</sup> ∇ · **<sup>u</sup>**<sup>I</sup> , from which *λ* is obtained. To complement (close) this equation, two types of boundary conditions are commonly used: homogeneous Dirichlet boundary conditions, *λ* = 0, for open or "flow through" boundaries (like truncated boundaries), and Neumann boundary conditions, *∂λ*/*∂***n** = 0, for closed or "no flow through" boundaries (like the surface terrain or topography). Many authors have been used and recommend these boundary conditions (Kitada et al., 1986; 1983; Ratto et al., 1994; Sherman, 1978). However they are physically and mathematically inconsistent as we will show in this work. Even though, there have been several sophisticated developments in the numerical simulations of this model as, for instance, the application of multigrid methods (Wang et al., 2005), and the application of genetic algorithms to estimate parameters (Montero et al., 2005), it seems that the analysis of boundary conditions has not attracted the attention of the community in meteorology.

In this work we study how boundary conditions affect solutions of the elliptic equation for *λ*. We show that the application of incorrect boundary conditions may degrade the solutions several orders of magnitude, and we propose some strategies to overcome this problem. In particular, we introduce a new approach based on the saddle–point formulation of the constrained least squares formulation of the problem, which allows the introduction of successful techniques from computational fluid dynamics. This new approach does not require boundary conditions for the multiplier. It produces much better results, and it also helps us to establish more consistent boundary conditions on truncated nonphysical boundaries. We also explore other boundary conditions for the multiplier better suited for artificial truncated boundaries. Furthermore, we present some preliminary numerical results using a meshfree method based on a radial basis function collocation method.

#### **2. Mathematical formulation of the problem**

Let Ω be an open, simply connected and bounded region in **R***<sup>d</sup>* (*d* = 2 or 3) with Lipschitz boundary *∂*Ω = Γ*<sup>N</sup>* ∪ Γ*D*, where Γ*<sup>N</sup>* is the part of the boundary associated to the surface terrain (topography), and Γ*<sup>D</sup>* is the rest of the boundary (artificial vertical boundaries and top boundary), as shown in Figure 1. Given an initial vector field **u**<sup>I</sup> in Ω (which can be obtained

Fig. 1. Bounded region Ω.

2 Will-be-set-by-IN-TECH

representing flows accurately in data–sparse regions as mountains or oceans. Despite these limitations, mass–consistent models are a valuable tool for air quality applications and consequently several developments have taken place over last decades (Ferragut et al., 2010; Ratto, 1996; Ratto et al., 1994; Ross et al., 1988; Wang et al., 2005). Most of the results presented in this chapter has been published in the last few years (Flores et al., 2010; Núñez et al., 2007;

The variational method proposed by Sasaki uses the continuity equation ∇ · **u** = 0, where **u** is the wind velocity vector field on a given domain Ω. The method is based on the minimization

where **u**<sup>I</sup> is an initial observed wind field, *λ* is a Lagrange multiplier and *S* is a diagonal matrix with weighting parameters *α<sup>i</sup>* > 0, *i* = 1, 2, 3, called Gaussian precision moduli, related to the scales of the respective components of the velocity field. The vertical component of the initial wind field is taken as zero because meteorological stations usually do not measure this

Usually **u** is obtained from (2), after *λ* is computed. Since ∇ · **u** = 0, then from (2) we obtain

(close) this equation, two types of boundary conditions are commonly used: homogeneous Dirichlet boundary conditions, *λ* = 0, for open or "flow through" boundaries (like truncated boundaries), and Neumann boundary conditions, *∂λ*/*∂***n** = 0, for closed or "no flow through" boundaries (like the surface terrain or topography). Many authors have been used and recommend these boundary conditions (Kitada et al., 1986; 1983; Ratto et al., 1994; Sherman, 1978). However they are physically and mathematically inconsistent as we will show in this work. Even though, there have been several sophisticated developments in the numerical simulations of this model as, for instance, the application of multigrid methods (Wang et al., 2005), and the application of genetic algorithms to estimate parameters (Montero et al., 2005), it seems that the analysis of boundary conditions has not attracted the attention of the

In this work we study how boundary conditions affect solutions of the elliptic equation for *λ*. We show that the application of incorrect boundary conditions may degrade the solutions several orders of magnitude, and we propose some strategies to overcome this problem. In particular, we introduce a new approach based on the saddle–point formulation of the constrained least squares formulation of the problem, which allows the introduction of successful techniques from computational fluid dynamics. This new approach does not require boundary conditions for the multiplier. It produces much better results, and it also helps us to establish more consistent boundary conditions on truncated nonphysical boundaries. We also explore other boundary conditions for the multiplier better suited for artificial truncated boundaries. Furthermore, we present some preliminary numerical results

using a meshfree method based on a radial basis function collocation method.

<sup>=</sup> ∇ · **<sup>u</sup>**<sup>I</sup>

+ *λ* [∇ · **u**]

**<sup>u</sup>** <sup>=</sup> **<sup>u</sup>**<sup>I</sup> <sup>+</sup> *<sup>S</sup>*−1∇*λ*, (2)

, from which *λ* is obtained. To complement

*dV* , (1)

2006), but we also include some additional ideas and recent results.

of the functional *L* defined by

the elliptic equation −∇ ·

community in meteorology.

*<sup>L</sup>* (**u**, *<sup>λ</sup>*) <sup>=</sup> <sup>1</sup>

component. The Euler–Lagrange equations of (1) are:

*<sup>S</sup>*−1∇*<sup>λ</sup>*

2 Ω *S* **<sup>u</sup>** <sup>−</sup> **<sup>u</sup>**<sup>I</sup> · **<sup>u</sup>** <sup>−</sup> **<sup>u</sup>**<sup>I</sup> 

> by interpolating atmospheric data, or by other means), our goal is to generate a solenoidal field **u** –called adjusted field– as close to **u**<sup>I</sup> as possible in a sense that will be clarified below, such that **u** · **n** = 0 on Γ*N*.

> We define the following vector function spaces: **L**2(Ω) = *L*2(Ω)*<sup>d</sup>* and **H**(*div*; Ω) = { **v** ∈ **L**2(Ω) : ∇ · **v** ∈ *L*2(Ω) }. Then, the adjusted wind field **u** must belong to the normed closed space

$$\mathbf{V} = \{ \mathbf{v} \in \mathbf{H}(div; \Omega) \; : \; \nabla \cdot \mathbf{v} = 0 \quad \text{and} \quad \mathbf{v} \cdot \mathbf{n} = 0 \quad \text{on} \quad \Gamma\_N \} \,\tag{3}$$

with the norm �·�*S*,<sup>Ω</sup> associated to the inner product �**u**, **<sup>v</sup>**�*<sup>S</sup>* <sup>=</sup> <sup>Ω</sup> *S***u** · **v** *d***x**, where **v** · **w** = ∑*d* <sup>1</sup> *vi wi* is the usual scalar product in **<sup>R</sup>***d*. We can now formulate the problem as a least squares projection problem. For this purpose, we define a convex quadratic functional *J* : **V** → **R** as

$$J(\mathbf{v}) = \frac{1}{2} \left\| \mathbf{v} - \mathbf{u}^{\mathrm{I}} \right\|\_{S, \Omega}^2 = \frac{1}{2} \int\_{\Omega} S(\mathbf{v} - \mathbf{u}^{\mathrm{I}}) \cdot (\mathbf{v} - \mathbf{u}^{\mathrm{I}}) \, d\mathbf{x}.\tag{4}$$

Then, our problem can be stated as follows:

Given **<sup>u</sup>**<sup>I</sup> <sup>∈</sup> **<sup>L</sup>**2(Ω), find **<sup>u</sup>** <sup>∈</sup> **<sup>V</sup>** such that *<sup>J</sup>*(**u**) <sup>≤</sup> *<sup>J</sup>*(**v**), <sup>∀</sup> **<sup>v</sup>** <sup>∈</sup> **<sup>V</sup>**. (5)

Due to the properties of this functional, **u** ∈ **V** is a minimizer of *J* if and only if it is a stationary point of *J*:

$$\frac{\partial}{\partial \boldsymbol{\varepsilon}} J(\mathbf{u} + \boldsymbol{\varepsilon} \cdot \mathbf{v})|\_{\boldsymbol{\varepsilon}=\boldsymbol{0}} = \int\_{\Omega} S(\mathbf{u} - \mathbf{u}^{\mathrm{I}}) \cdot \mathbf{v} \, d\mathbf{x} = \mathbf{0}, \quad \forall \, \mathbf{v} \in \mathbf{V}. \tag{6}$$

The Lax–Milgram theorem guaranties that this equation has a unique solution.

respectively. Thus, the finite element algorithm is: Given **u**<sup>I</sup>

*<sup>S</sup>*−1∇*λ<sup>h</sup>* · ∇*q d***<sup>x</sup>** <sup>=</sup> <sup>−</sup>

*<sup>h</sup>*) · **v** *d***x** −

<sup>1</sup> 1.2 1.4 1.6 1.8 <sup>2</sup> 2.2 2.4 -0.2

example, we obtain *er* <sup>=</sup> 1.9 <sup>×</sup> <sup>10</sup>−<sup>2</sup> and *mdiv* <sup>=</sup> 4.1 <sup>×</sup> <sup>10</sup>−2.

(1, 2) <sup>×</sup> (0, 1), so that **<sup>u</sup>** <sup>∈</sup> **<sup>V</sup>**. Assuming that we have **<sup>u</sup>**<sup>I</sup>

From now on, we identify the algorithm (14)–(15) as the *E1–algorithm*.

*<sup>h</sup>* <sup>∈</sup> **<sup>L</sup>***<sup>h</sup>* is the interpolant of the given initial velocity field **<sup>u</sup>**<sup>I</sup>

 Ω  Ω **u**I

Mass–Consistent Wind Field Models: Numerical Techniques by L2–Projection Methods 27

solving the resulting system of linear equations, and the numerical approximation **u***<sup>h</sup>* of **u** is computed by the weak version of (2) as follows: Find **u***<sup>h</sup>* ∈ **L***<sup>h</sup>* with **u***<sup>h</sup>* · **n** = 0 on Γ*<sup>N</sup>* such that

**Example 1**. We consider the solenoidal vector field **u**(*x*, *z*)=(*x*, −*z*) defined in Ω =

wind field, we want to apply the *E1–algorithm* to see how much we can recover of the vertical component of **u**. For this numerical calculation, Ω is divided into a 80 × 80 triangular mesh, and we choose the following values for the Gaussian Precision moduli: *α*<sup>1</sup> = 1 and *α*<sup>3</sup> = 0.001. Figure 2 shows the exact field in red and the computed adjusted field in blue. Both fields agree fairly well almost everywhere, except on the vertical artificial boundaries *x* = 1 and *x* = 2.

 Ω

 Ω (*S***u**<sup>I</sup>

where **u**<sup>I</sup>

 Ω

(*S***u***h*) · **v** *d***x** =

0 0.2 0.4 0.6 0.8 1 1.2

*er* <sup>=</sup> ||**<sup>u</sup>** <sup>−</sup> **<sup>u</sup>***h*||<sup>2</sup> ||**u**||<sup>2</sup>

*<sup>h</sup>* ∈ **L***h*, find *λh*∈*Hh* such that

*<sup>h</sup>* · ∇*q d***x**, ∀ *q* ∈ *Hh*, (14)

(*x*, *z*)=(*x*, 0) as an initial horizontal

*λh*∇ · **v** *d***x**, ∀ **v** ∈ **L***h*, **v** · **n** = 0 on Γ*N*. (15)

1.85 1.9 1.95 2 2.05

{∇· **u***h*(**x***i*) | **x***<sup>i</sup>* is a interior vertex}, (16)

**u***<sup>h</sup>* · ∇*φ<sup>i</sup> d***x** , (17)

0.8 0.82 0.84 0.86 0.88 0.9 0.92 0.94 0.96 0.98 1

Fig. 2. Exact field **u** = (*x*, −*z*) in red, adjusted field obtained by the *E1–algorithm* in blue.

 Ω

where *φ<sup>i</sup>* is the piece–wise linear base function associated to vertex node **x***i*. For the present

The values for the Gaussian precision moduli were chosen based on numerical performance. Table 1 shows the behavior of *er* and *mdiv* for different values of *α*<sup>3</sup> when *α*<sup>1</sup> is kept constant

The relative error and the mean divergence of the computed solution are defined as

respectively. The point–wise divergence is computed in a weak sense, as follows

∇ · **u***h*(**x***i*) = −

, and *mdiv* <sup>=</sup> mean **<sup>x</sup>***<sup>i</sup>*

. We obtain *λ<sup>h</sup>* after

#### **3. The traditional approach. Advantages and difficulties**

#### **3.1 Derivation of the elliptic problem**

The first approach is based on a Helmholtz–type decomposition of the Hilbert vector space **L**2(Ω), and it reduces to the traditional approach used by meteorologists. **Proposition 1** *The orthogonal complement in* **L**2(Ω) *of the closed subspace* **V** *is*

$$\mathbf{V}^{\perp} = \{ \ \nabla \boldsymbol{\eta} \; : \; \boldsymbol{\eta} \in H^{1}(\Omega), \quad \boldsymbol{\eta} = 0 \quad \text{on} \quad \Gamma\_{D} \}.$$

An argument very similar to that given by Girault and Raviart (Girault & Raviart, 1986), shows that this decomposition is valid (details are given in (Núñez et al., 2007)). Therefore we get from (6) that *S* **<sup>u</sup>** <sup>−</sup> **<sup>u</sup>**<sup>I</sup> = ∇*λ*, with *λ* in

$$H\_D^1(\Omega) \equiv \{ \: q \in H^1(\Omega) \; : \; q = 0 \,\text{on} \,\Gamma\_D \}. \tag{7}$$

With the above properties, we obtain a saddle–point problem for **u** and *λ* (left), as well as the correspondent elliptic problem for *λ* (right):

$$\mathbf{S} \cdot \mathbf{S} \mathbf{u} - \nabla \lambda = \mathbf{S} \mathbf{u}^{\mathrm{I}} \quad \text{and} \quad \nabla \cdot \mathbf{u} = 0 \quad \text{in } \Omega, \qquad -\nabla \cdot \left( \mathbf{S}^{-1} \nabla \lambda \right) = \nabla \cdot \mathbf{u}^{\mathrm{I}} \quad \text{in } \Omega, \tag{8}$$

$$
\lambda = 0 \quad \text{on } \Gamma\_{D\prime} \tag{9}
$$

$$\mathbf{u} \cdot \mathbf{n} = 0 \quad \text{on } \Gamma\_N. \qquad \qquad \qquad -\mathbb{S}^{-1} \nabla \lambda \cdot \mathbf{n} = \mathbf{u}^\mathrm{I} \cdot \mathbf{n} \quad \text{on} \quad \Gamma\_N. \tag{10}$$

To obtain the elliptic problem, we eliminate **u** from the saddle–point problem using that **u** belongs to **V**. Once *λ* is calculated from (8)–(10), the adjusted field is recovered from (2).

Equation (8) has traditionally been used by meteorologists. However, this equation is generally introduced from a discussion in which it is not clear how to establish the proper boundary conditions for *λ*. The crucial argument in our study is the decomposition of **L**2(Ω) in orthogonal subspaces **V** and **V**⊥, from which the boundary conditions for *λ* arises in a natural way, from the mathematical point of view. We would like to mention that the boundary condition (10) has already been used in recent research (Ferragut et al., 2010).

#### **3.2 Finite element solution of the elliptic problem**

The variational formulation of the elliptic problem (8)–(10) is

$$\int\_{\Omega} \mathbf{S}^{-1} \nabla \lambda \cdot \nabla q \, d\mathbf{x} = -\int\_{\Omega} \mathbf{u}^{\mathrm{I}} \cdot \nabla q \, d\mathbf{x}, \quad \forall q \in H\_D^{1}(\Omega). \tag{11}$$

Here, we consider the two–dimensional case. Let T*<sup>h</sup>* be a finite element triangulation of <sup>Ω</sup> <sup>⊂</sup> **<sup>R</sup>**<sup>2</sup> (Ciarlet, 2002), where *<sup>h</sup>* is taken as the space discretization step. Let's denote by *P*<sup>1</sup> the space of polynomials of degree less or equal than 1. Then, **L**2(Ω) and *H*<sup>1</sup> *<sup>D</sup>*(Ω) are approximated by the finite dimensional spaces

$$\mathbf{L}\_{\hbar} = \left\{ \mathbf{v}\_{\hbar} \in \mathcal{C}^{0}(\tilde{\Omega})^{2} \; : \; \mathbf{v}\_{\hbar}|\_{T} \in P\_{1} \times P\_{1\prime} \; \forall \; T \in \mathcal{T}\_{\hbar} \right\},\tag{12}$$

$$H\_{\hbar} = \left\{ \begin{array}{l} q \in \mathcal{C}^{0}(\vec{\Omega}) \;:\; q \vert\_{T} \in P\_{1\prime} \; \forall \; T \in \mathcal{T}\_{\hbar \prime} \; q = 0 \text{ on } \Gamma\_{D} \right\} \end{array} \tag{13}$$

4 Will-be-set-by-IN-TECH

The first approach is based on a Helmholtz–type decomposition of the Hilbert vector space

**<sup>V</sup>**<sup>⊥</sup> <sup>=</sup> { ∇*<sup>q</sup>* : *<sup>q</sup>* <sup>∈</sup> *<sup>H</sup>*1(Ω), *<sup>q</sup>* <sup>=</sup> 0 on <sup>Γ</sup>*<sup>D</sup>* }.

An argument very similar to that given by Girault and Raviart (Girault & Raviart, 1986), shows that this decomposition is valid (details are given in (Núñez et al., 2007)). Therefore

With the above properties, we obtain a saddle–point problem for **u** and *λ* (left), as well as the

To obtain the elliptic problem, we eliminate **u** from the saddle–point problem using that **u** belongs to **V**. Once *λ* is calculated from (8)–(10), the adjusted field is recovered from (2). Equation (8) has traditionally been used by meteorologists. However, this equation is generally introduced from a discussion in which it is not clear how to establish the proper boundary conditions for *λ*. The crucial argument in our study is the decomposition of **L**2(Ω) in orthogonal subspaces **V** and **V**⊥, from which the boundary conditions for *λ* arises in a natural way, from the mathematical point of view. We would like to mention that the boundary condition (10) has already been used in recent research (Ferragut et al., 2010).

> Ω

*P*<sup>1</sup> the space of polynomials of degree less or equal than 1. Then, **L**2(Ω) and *H*<sup>1</sup>

Here, we consider the two–dimensional case. Let T*<sup>h</sup>* be a finite element triangulation of <sup>Ω</sup> <sup>⊂</sup> **<sup>R</sup>**<sup>2</sup> (Ciarlet, 2002), where *<sup>h</sup>* is taken as the space discretization step. Let's denote by

<sup>2</sup> : **<sup>v</sup>***h*|*<sup>T</sup>* <sup>∈</sup> *<sup>P</sup>*<sup>1</sup> <sup>×</sup> *<sup>P</sup>*1, <sup>∀</sup> *<sup>T</sup>* ∈ T*<sup>h</sup>*

*<sup>q</sup>* ∈ C0(Ω¯ ) : *<sup>q</sup>*|*<sup>T</sup>* <sup>∈</sup> *<sup>P</sup>*1, <sup>∀</sup> *<sup>T</sup>* ∈ T*h*, *<sup>q</sup>* <sup>=</sup> 0 on <sup>Γ</sup>*<sup>D</sup>*

*<sup>D</sup>*(Ω) ≡ { *<sup>q</sup>* <sup>∈</sup> *<sup>H</sup>*1(Ω) : *<sup>q</sup>* <sup>=</sup> 0 on <sup>Γ</sup>*<sup>D</sup>* }. (7)

<sup>=</sup> ∇ · **<sup>u</sup>**<sup>I</sup> in <sup>Ω</sup>, (8)

*<sup>D</sup>*(Ω). (11)

, (12)

, (13)

*<sup>D</sup>*(Ω) are

 *<sup>S</sup>*−1∇*<sup>λ</sup>* 

**<sup>u</sup>**<sup>I</sup> · ∇*q d***x**, <sup>∀</sup> *<sup>q</sup>* <sup>∈</sup> *<sup>H</sup>*<sup>1</sup>

*λ* = 0 on Γ*D*, *λ* = 0 on Γ*D*, (9) **<sup>u</sup>** · **<sup>n</sup>** <sup>=</sup> 0 on <sup>Γ</sup>*N*. <sup>−</sup>*S*−1∇*<sup>λ</sup>* · **<sup>n</sup>** <sup>=</sup> **<sup>u</sup>**<sup>I</sup> · **<sup>n</sup>** on <sup>Γ</sup>*N*. (10)

**3. The traditional approach. Advantages and difficulties**

**L**2(Ω), and it reduces to the traditional approach used by meteorologists. **Proposition 1** *The orthogonal complement in* **L**2(Ω) *of the closed subspace* **V** *is*

= ∇*λ*, with *λ* in

, and ∇ · **u** = 0 in Ω, −∇ ·

**3.1 Derivation of the elliptic problem**

 **<sup>u</sup>** <sup>−</sup> **<sup>u</sup>**<sup>I</sup> 

correspondent elliptic problem for *λ* (right):

**3.2 Finite element solution of the elliptic problem**

 Ω

approximated by the finite dimensional spaces

**L***<sup>h</sup>* = 

*Hh* = 

The variational formulation of the elliptic problem (8)–(10) is

*<sup>S</sup>*−1∇*<sup>λ</sup>* · ∇*q d***<sup>x</sup>** <sup>=</sup> <sup>−</sup>

**<sup>v</sup>***<sup>h</sup>* ∈ C0(Ω¯ )

*H*1

we get from (6) that *S*

*<sup>S</sup>***<sup>u</sup>** − ∇*<sup>λ</sup>* <sup>=</sup> *<sup>S</sup>***u**<sup>I</sup>

respectively. Thus, the finite element algorithm is: Given **u**<sup>I</sup> *<sup>h</sup>* ∈ **L***h*, find *λh*∈*Hh* such that

$$\int\_{\Omega} S^{-1} \nabla \lambda\_{\hbar} \cdot \nabla \eta \, d\mathbf{x} = - \int\_{\Omega} \mathbf{u}\_{\hbar}^{\mathrm{I}} \cdot \nabla \eta \, d\mathbf{x}, \quad \forall \, \eta \in H\_{\hbar \nu} \tag{14}$$

where **u**<sup>I</sup> *<sup>h</sup>* <sup>∈</sup> **<sup>L</sup>***<sup>h</sup>* is the interpolant of the given initial velocity field **<sup>u</sup>**<sup>I</sup> . We obtain *λ<sup>h</sup>* after solving the resulting system of linear equations, and the numerical approximation **u***<sup>h</sup>* of **u** is computed by the weak version of (2) as follows: Find **u***<sup>h</sup>* ∈ **L***<sup>h</sup>* with **u***<sup>h</sup>* · **n** = 0 on Γ*<sup>N</sup>* such that

$$\int\_{\Omega} (\mathbf{S} \mathbf{u}\_h) \cdot \mathbf{v} \, d\mathbf{x} = \int\_{\Omega} (\mathbf{S} \mathbf{u}\_h^\mathrm{I}) \cdot \mathbf{v} \, d\mathbf{x} - \int\_{\Omega} \lambda\_h \nabla \cdot \mathbf{v} \, d\mathbf{x}, \quad \forall \, \mathbf{v} \in \mathbf{L}\_h, \quad \mathbf{v} \cdot \mathbf{n} = 0 \quad \text{on} \quad \Gamma\_N. \tag{15}$$

From now on, we identify the algorithm (14)–(15) as the *E1–algorithm*.

**Example 1**. We consider the solenoidal vector field **u**(*x*, *z*)=(*x*, −*z*) defined in Ω = (1, 2) <sup>×</sup> (0, 1), so that **<sup>u</sup>** <sup>∈</sup> **<sup>V</sup>**. Assuming that we have **<sup>u</sup>**<sup>I</sup> (*x*, *z*)=(*x*, 0) as an initial horizontal wind field, we want to apply the *E1–algorithm* to see how much we can recover of the vertical component of **u**. For this numerical calculation, Ω is divided into a 80 × 80 triangular mesh, and we choose the following values for the Gaussian Precision moduli: *α*<sup>1</sup> = 1 and *α*<sup>3</sup> = 0.001. Figure 2 shows the exact field in red and the computed adjusted field in blue. Both fields agree fairly well almost everywhere, except on the vertical artificial boundaries *x* = 1 and *x* = 2.

Fig. 2. Exact field **u** = (*x*, −*z*) in red, adjusted field obtained by the *E1–algorithm* in blue. The relative error and the mean divergence of the computed solution are defined as

$$\varepsilon\_{r} = \frac{||\mathbf{u} - \mathbf{u}\_{h}||\_{2}}{|||\mathbf{u}|||\_{2}}, \quad \text{and} \quad m d\dot{v} = \underset{\mathbf{x}\_{i}}{\text{mean}} \{ |\nabla \cdot \mathbf{u}\_{h}(\mathbf{x}\_{i})| \, | \, \mathbf{x}\_{i} \text{ is a interior vertex} \}, \tag{16}$$

respectively. The point–wise divergence is computed in a weak sense, as follows

$$\nabla \cdot \mathbf{u}\_h(\mathbf{x}\_i) = -\int\_{\Omega} \mathbf{u}\_h \cdot \nabla \phi\_i \, d\mathbf{x}\_\prime \tag{17}$$

where *φ<sup>i</sup>* is the piece–wise linear base function associated to vertex node **x***i*. For the present example, we obtain *er* <sup>=</sup> 1.9 <sup>×</sup> <sup>10</sup>−<sup>2</sup> and *mdiv* <sup>=</sup> 4.1 <sup>×</sup> <sup>10</sup>−2.

The values for the Gaussian precision moduli were chosen based on numerical performance. Table 1 shows the behavior of *er* and *mdiv* for different values of *α*<sup>3</sup> when *α*<sup>1</sup> is kept constant

Furthermore, **u***<sup>λ</sup>* must satisfy (20) which has the following equivalent strong version

introduce the linear operator *A* from *L*2(Ω) into *L*2(Ω) defined by

*S***u***<sup>q</sup>* · **v** *d***x** = −

<sup>Ω</sup> *<sup>q</sup>*� ∇·**u***<sup>q</sup>* <sup>=</sup>

*<sup>S</sup>* **<sup>u</sup>***<sup>q</sup>* · **<sup>u</sup>***<sup>q</sup>* <sup>&</sup>gt; *<sup>c</sup>* �**u***q*�<sup>2</sup>

 Ω

*<sup>S</sup>* **<sup>u</sup>***<sup>k</sup>* · **<sup>v</sup>** *<sup>d</sup>***<sup>x</sup>** <sup>=</sup> <sup>−</sup>

Compute *λk*+<sup>1</sup> = *λ<sup>k</sup>* + *α<sup>k</sup> dk*, **u***k*+<sup>1</sup> = **u***<sup>k</sup>* + *α<sup>k</sup>* **u***k*, *gk*+<sup>1</sup> = *g<sup>k</sup>* + *α<sup>k</sup> wk*.

3. If �*gk*, *<sup>g</sup>k*� ≤ *<sup>ε</sup>*�*g*0, *<sup>g</sup>*0�, take *<sup>λ</sup>* <sup>=</sup> *<sup>λ</sup>k*+<sup>1</sup> and **<sup>u</sup>** <sup>=</sup> **<sup>u</sup>***k*+<sup>1</sup> and stop. Otherwise, compute

where **u***<sup>q</sup>* ∈ **V***<sup>N</sup>* is the solution of

**4.2 Conjugate gradient algorithm**

<sup>Ω</sup> *<sup>q</sup>*� *Aq d***<sup>x</sup>** <sup>=</sup> <sup>−</sup>

*qAq d***x** =

1. Given *<sup>λ</sup>*<sup>0</sup> <sup>∈</sup> *<sup>L</sup>*2(Ω), solve for **<sup>u</sup>**<sup>0</sup> <sup>∈</sup> **<sup>V</sup>***<sup>N</sup>* Ω

Set *<sup>g</sup>*<sup>0</sup> <sup>=</sup> −∇ · **<sup>u</sup>**<sup>0</sup> and *<sup>d</sup>*<sup>0</sup> <sup>=</sup> <sup>−</sup>*g*0.

infinite dimensional problem (25):

following: Solve for **<sup>u</sup>***k*<sup>∈</sup> **<sup>V</sup>***<sup>N</sup>*

 Ω

equation

 Ω

*<sup>S</sup>* **<sup>u</sup>**<sup>0</sup> · **<sup>v</sup>** *<sup>d</sup>***<sup>x</sup>** <sup>=</sup>

 Ω

Set *<sup>w</sup><sup>k</sup>* <sup>=</sup> −∇ · **<sup>u</sup>***<sup>k</sup>* and compute *<sup>α</sup><sup>k</sup>* <sup>=</sup> �*gk*, *<sup>g</sup>k*�

 Ω −∇· **<sup>u</sup>***<sup>λ</sup>* <sup>=</sup> ∇ · **<sup>u</sup>**<sup>I</sup>

Mass–Consistent Wind Field Models: Numerical Techniques by L2–Projection Methods 29

A key point is that problem (21)–(22) can be formulated as a functional equation. For this we

 Ω

With this definition, it is clear, from (21)–(22), that the multiplier *λ* satisfies the functional

*<sup>A</sup><sup>λ</sup>* <sup>=</sup> ∇ · **<sup>u</sup>**<sup>I</sup>

Therefore, the following iterative conjugate gradient algorithm may be used to solve the

 Ω

*<sup>S</sup>* **<sup>u</sup>***<sup>I</sup>* · **<sup>v</sup>** *<sup>d</sup>***<sup>x</sup>** <sup>−</sup>

2. For *<sup>k</sup>* <sup>≥</sup> 0, assuming we know *<sup>λ</sup>k*, *<sup>g</sup>k*, *<sup>d</sup>k*, **<sup>u</sup>***k*, find *<sup>λ</sup>k*<sup>+</sup>1, *<sup>g</sup>k*<sup>+</sup>1, *<sup>d</sup>k*<sup>+</sup>1, **<sup>u</sup>***k*<sup>+</sup>1, doing the

�*dk*, *<sup>w</sup>k*�

*<sup>d</sup>k*+<sup>1</sup> <sup>=</sup> <sup>−</sup>*gk*+<sup>1</sup> <sup>+</sup> *<sup>β</sup><sup>k</sup> <sup>d</sup><sup>k</sup>* where *<sup>β</sup><sup>k</sup>* <sup>=</sup> �*gk*<sup>+</sup>1, *<sup>g</sup>k*+1�

.

 Ω

Operator *A* is selfadjoint, and strongly elliptic, since from (23) and (24) we have

 Ω . (22)

*A q* = −∇ · **u***<sup>q</sup>* , (23)

*S* **u***q*� · **u***<sup>q</sup> d***x** ∀ *q*, *q*� ∈ *L*2(Ω),

*<sup>d</sup>k*∇ · **<sup>v</sup>** *<sup>d</sup>***x**, <sup>∀</sup> **<sup>v</sup>** <sup>∈</sup> **<sup>V</sup>***N*.

*<sup>L</sup>*2(Ω) ∀ *q* �= 0 (0 < *c* < min{*αi*})

*<sup>λ</sup>*0∇ · **<sup>v</sup>** *<sup>d</sup>***x**, <sup>∀</sup> **<sup>v</sup>** <sup>∈</sup> **<sup>V</sup>***N*.

�*gk*, *<sup>g</sup>k*� .

*q*∇ · **v** *d***x**, ∀ **v** ∈ **V***N*. (24)

. (25)

and equal to one. Clearly the best results were obtained with *α*<sup>3</sup> = 0.001. We will explain this behavior later on. But, for the moment, we want to emphasize that this algorithm produces satisfactory results almost everywhere, except on the boundary Γ*D*, where homogeneous Dirchlet boundary conditions were imposed.


Table 1. Numerical performance of *E1–algorithm* for different values of *α*3.

We can say that the main advantage of this traditional way to solve the problem is its simplicity, since it only involves the solutions of an elliptic partial differential equation (PDE). On the other hand, one of its major drawbacks is that inconsistent or incorrect boundary conditions, on truncated artificial boundaries, degrade the accuracy of the solution. In the rest of the chapter, we introduce some alternatives to overcome these problems.

#### **4. A saddle–point formulation and the conjugate gradient algorithm**

#### **4.1 Derivation of the formulation**

The second approach to solve the problem (5), or equivalently problem (6), is based on the usual methodology to solve constrained optimization problems. That is, we introduce the space of vector functions

$$\mathbf{V}\_N = \{ \mathbf{v} \in \mathbf{H}(div; \Omega) \; : \; \mathbf{v} \cdot \mathbf{n} = 0 \; \text{on } \Gamma\_N \} \; \; \; \tag{18}$$

together with the Lagrangian *L* defined on **V***<sup>N</sup>* × *L*2(Ω) as

$$L(\mathbf{v}, q) \equiv f(\mathbf{v}) + \langle q\_\prime \nabla \cdot \mathbf{v} \rangle = \frac{1}{2} \int\_{\Omega} S(\mathbf{v} - \mathbf{u}^\mathrm{I}) \cdot (\mathbf{v} - \mathbf{u}^\mathrm{I}) \, d\mathbf{x} + \int\_{\Omega} q \nabla \cdot \mathbf{v} \, d\mathbf{x}.$$

A stationary point (**u**, *λ*) of *L* solves the following saddle–point problem

$$
\int\_{\Omega} \mathbf{S} \mathbf{u} \cdot \mathbf{v} \, d\mathbf{x} + \int\_{\Omega} \lambda \nabla \cdot \mathbf{v} \, d\mathbf{x} = \int\_{\Omega} \mathbf{S} \mathbf{u}^{\mathrm{I}} \cdot \mathbf{v} \, d\mathbf{x}, \quad \forall \, \mathbf{v} \in \mathbf{V}\_{N\prime} \tag{19}
$$

$$\int\_{\Omega} q \nabla \cdot \mathbf{u} \, d\mathbf{x} = 0, \quad \forall q \in L\_2(\Omega), \tag{20}$$

where *λ* **need not satisfy boundary conditions**. The solution **u** is the minimizer of *J*, and now it is obtained from the enlarged space **V***<sup>N</sup>* where free divergence is not required. Instead, the condition ∇ · **u** = 0 is relaxed by the introduction of the Lagrange multiplier *λ* so that **u** must satisfy the weaker condition (20). To solve (19)–(20) we introduce a method which has shown to be very effective for solving Stokes problems in computational fluid dynamics (Glowinski, 2003). The idea is as follows: assuming that (**u**, *λ*) is a solution of the problem (19)–(20), the vector field **u** is decomposed as **u** = **u**<sup>I</sup> + **u***λ*, where **u**<sup>I</sup> is the given initial vector field, and **u***<sup>λ</sup>* ∈ **V***<sup>N</sup>* solves

$$
\int\_{\Omega} \mathbf{S} \mathbf{u}\_{\lambda} \cdot \mathbf{v} \, d\mathbf{x} = -\int\_{\Omega} \lambda \nabla \cdot \mathbf{v} \, d\mathbf{x}, \quad \forall \mathbf{v} \in \mathbf{V}\_{N}. \tag{21}
$$

Furthermore, **u***<sup>λ</sup>* must satisfy (20) which has the following equivalent strong version

$$-\nabla \cdot \mathbf{u}\_{\lambda} = \nabla \cdot \mathbf{u}^{\mathsf{I}}.\tag{22}$$

A key point is that problem (21)–(22) can be formulated as a functional equation. For this we introduce the linear operator *A* from *L*2(Ω) into *L*2(Ω) defined by

$$A\,\boldsymbol{q} = -\nabla \cdot \mathbf{u}\_{\boldsymbol{q}\,'} \tag{23}$$

where **u***<sup>q</sup>* ∈ **V***<sup>N</sup>* is the solution of

6 Will-be-set-by-IN-TECH

and equal to one. Clearly the best results were obtained with *α*<sup>3</sup> = 0.001. We will explain this behavior later on. But, for the moment, we want to emphasize that this algorithm produces satisfactory results almost everywhere, except on the boundary Γ*D*, where homogeneous

α<sup>3</sup> 0.001 0.01 0.1 1 100 1000 *er* 1.9 <sup>×</sup> <sup>10</sup>−<sup>2</sup> 9.6 <sup>×</sup> <sup>10</sup>−<sup>2</sup> 1.4 <sup>×</sup> <sup>10</sup>−<sup>1</sup> 5.2 <sup>×</sup> <sup>10</sup>−<sup>1</sup> 6.4 <sup>×</sup> <sup>10</sup>−<sup>1</sup> 9.8 <sup>×</sup> <sup>10</sup>−<sup>1</sup> *mdiv* 4.1 <sup>×</sup> <sup>10</sup>−<sup>2</sup> <sup>−</sup>6.1 <sup>×</sup> <sup>10</sup>−<sup>2</sup> 2.9 <sup>×</sup> <sup>10</sup>−<sup>1</sup> 5.4 <sup>×</sup> <sup>10</sup>−<sup>1</sup> 7.8 <sup>×</sup> <sup>10</sup>−<sup>1</sup> 9.8 <sup>×</sup> <sup>10</sup>−<sup>1</sup>

We can say that the main advantage of this traditional way to solve the problem is its simplicity, since it only involves the solutions of an elliptic partial differential equation (PDE). On the other hand, one of its major drawbacks is that inconsistent or incorrect boundary conditions, on truncated artificial boundaries, degrade the accuracy of the solution. In the

The second approach to solve the problem (5), or equivalently problem (6), is based on the usual methodology to solve constrained optimization problems. That is, we introduce the

*<sup>S</sup>*(**<sup>v</sup>** <sup>−</sup> **<sup>u</sup>**<sup>I</sup>

 Ω

where *λ* **need not satisfy boundary conditions**. The solution **u** is the minimizer of *J*, and now it is obtained from the enlarged space **V***<sup>N</sup>* where free divergence is not required. Instead, the condition ∇ · **u** = 0 is relaxed by the introduction of the Lagrange multiplier *λ* so that **u** must satisfy the weaker condition (20). To solve (19)–(20) we introduce a method which has shown to be very effective for solving Stokes problems in computational fluid dynamics (Glowinski, 2003). The idea is as follows: assuming that (**u**, *λ*) is a solution of the problem (19)–(20), the vector field **u** is decomposed as **u** = **u**<sup>I</sup> + **u***λ*, where **u**<sup>I</sup> is the given initial vector field, and

> Ω

**V***<sup>N</sup>* = { **v** ∈ **H**(*div*; Ω) : **v** · **n** = 0 on Γ*<sup>N</sup>* } , (18)

) · (**<sup>v</sup>** <sup>−</sup> **<sup>u</sup>**<sup>I</sup>

) *d***x** + Ω

*q*∇ · **u** *d***x** = 0, ∀ *q* ∈ *L*2(Ω), (20)

*q*∇ · **v** *d***x**.

*<sup>S</sup>***u**<sup>I</sup> · **<sup>v</sup>** *<sup>d</sup>***x**, <sup>∀</sup> **<sup>v</sup>** <sup>∈</sup> **<sup>V</sup>***N*, (19)

*λ*∇ · **v** *d***x**, ∀ **v** ∈ **V***N*. (21)

Table 1. Numerical performance of *E1–algorithm* for different values of *α*3.

rest of the chapter, we introduce some alternatives to overcome these problems.

2 Ω

*λ*∇ · **v** *d***x** =

A stationary point (**u**, *λ*) of *L* solves the following saddle–point problem

 Ω

 Ω

*S***u***<sup>λ</sup>* · **v** *d***x** = −

**4. A saddle–point formulation and the conjugate gradient algorithm**

together with the Lagrangian *L* defined on **V***<sup>N</sup>* × *L*2(Ω) as

*<sup>L</sup>*(**v**, *<sup>q</sup>*) <sup>≡</sup> *<sup>J</sup>*(**v**) + �*q*, ∇ · **<sup>v</sup>**� <sup>=</sup> <sup>1</sup>

*S***u** · **v** *d***x** +

 Ω

 Ω

Dirchlet boundary conditions were imposed.

**4.1 Derivation of the formulation**

space of vector functions

**u***<sup>λ</sup>* ∈ **V***<sup>N</sup>* solves

$$
\int\_{\Omega} \mathbf{S} \mathbf{u}\_q \cdot \mathbf{v} \, d\mathbf{x} = - \int\_{\Omega} q \nabla \cdot \mathbf{v} \, d\mathbf{x}, \quad \forall \mathbf{v} \in \mathbf{V}\_N. \tag{24}
$$

With this definition, it is clear, from (21)–(22), that the multiplier *λ* satisfies the functional equation

$$A\lambda = \nabla \cdot \mathbf{u}^{\mathsf{I}}.\tag{25}$$

#### **4.2 Conjugate gradient algorithm**

Operator *A* is selfadjoint, and strongly elliptic, since from (23) and (24) we have

$$\begin{aligned} \int\_{\Omega} q' A \, q \, \, d\mathbf{x} &= - \int\_{\Omega} q' \, \nabla \cdot \mathbf{u}\_{\boldsymbol{q}} = \int\_{\Omega} S \, \mathbf{u}\_{\boldsymbol{q}'} \cdot \mathbf{u}\_{\boldsymbol{q}} \, d\mathbf{x} & \quad \forall \, \boldsymbol{q}, \, \boldsymbol{q}' \in L\_2(\Omega), \\\int\_{\Omega} q \, A \, \boldsymbol{q} \, \, d\mathbf{x} &= \int\_{\Omega} S \, \mathbf{u}\_{\boldsymbol{q}} \cdot \mathbf{u}\_{\boldsymbol{q}} > c \, \|\mathbf{u}\_{\boldsymbol{q}}\|\_{L\_2(\Omega)}^2 \qquad \forall \, \boldsymbol{q} \neq 0 \quad (\boldsymbol{0} < c < \min\{a\_i\}). \end{aligned}$$

Therefore, the following iterative conjugate gradient algorithm may be used to solve the infinite dimensional problem (25):

1. Given *<sup>λ</sup>*<sup>0</sup> <sup>∈</sup> *<sup>L</sup>*2(Ω), solve for **<sup>u</sup>**<sup>0</sup> <sup>∈</sup> **<sup>V</sup>***<sup>N</sup>*

$$\int\_{\Omega} \mathbf{S} \cdot \mathbf{u}^{0} \cdot \mathbf{v} \, d\mathbf{x} = \int\_{\Omega} \mathbf{S} \cdot \mathbf{u}^{I} \cdot \mathbf{v} \, d\mathbf{x} - \int\_{\Omega} \boldsymbol{\lambda}^{0} \nabla \cdot \mathbf{v} \, d\mathbf{x}, \quad \forall \, \mathbf{v} \in \mathbf{V}\_{N}.$$
  $\text{Set } \mathbf{g}^{0} = -\nabla \cdot \mathbf{u}^{0} \quad \text{and} \quad d^{0} = -\mathbf{g}^{0}.$ 

2. For *<sup>k</sup>* <sup>≥</sup> 0, assuming we know *<sup>λ</sup>k*, *<sup>g</sup>k*, *<sup>d</sup>k*, **<sup>u</sup>***k*, find *<sup>λ</sup>k*<sup>+</sup>1, *<sup>g</sup>k*<sup>+</sup>1, *<sup>d</sup>k*<sup>+</sup>1, **<sup>u</sup>***k*<sup>+</sup>1, doing the following: Solve for **<sup>u</sup>***k*<sup>∈</sup> **<sup>V</sup>***<sup>N</sup>*

$$\int\_{\Omega} \mathbf{S} \underbrace{\mathbf{u}}\_{}^{k} \cdot \mathbf{v} \, d\mathbf{x} = - \int\_{\Omega} \, d^{k} \nabla \cdot \mathbf{v} \, d\mathbf{x}, \quad \forall \, \mathbf{v} \in \mathbf{V}\_{N}.$$

Set *<sup>w</sup><sup>k</sup>* <sup>=</sup> −∇ · **<sup>u</sup>***<sup>k</sup>* and compute *<sup>α</sup><sup>k</sup>* <sup>=</sup> �*gk*, *<sup>g</sup>k*� �*dk*, *<sup>w</sup>k*� . Compute *λk*+<sup>1</sup> = *λ<sup>k</sup>* + *α<sup>k</sup> dk*, **u***k*+<sup>1</sup> = **u***<sup>k</sup>* + *α<sup>k</sup>* **u***k*, *gk*+<sup>1</sup> = *g<sup>k</sup>* + *α<sup>k</sup> wk*.

3. If �*gk*, *<sup>g</sup>k*� ≤ *<sup>ε</sup>*�*g*0, *<sup>g</sup>*0�, take *<sup>λ</sup>* <sup>=</sup> *<sup>λ</sup>k*+<sup>1</sup> and **<sup>u</sup>** <sup>=</sup> **<sup>u</sup>***k*+<sup>1</sup> and stop. Otherwise, compute

$$d^{k+1} = -g^{k+1} + \beta\_k d^k \quad \text{where} \quad \beta\_k = \frac{\langle g^{k+1}, g^{k+1} \rangle}{\langle g^k, g^k \rangle}.$$

<sup>1</sup> 1.2 1.4 1.6 1.8 <sup>2</sup> 2.2 2.4 -0.2

<sup>0</sup> <sup>2</sup> <sup>4</sup> <sup>6</sup> <sup>8</sup> <sup>10</sup> <sup>12</sup> -2

Ω =  1.85 1.9 1.95 2 2.05

0.8 0.82 0.84 0.86 0.88 0.9 0.92 0.94 0.96 0.98 1

<sup>2</sup> cos

8 8.2 8.4 8.6 8.8 9 9.2 9.4 9.6 9.8 10 3*πx* <sup>10</sup> <sup>+</sup> 1

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6

,

<sup>2</sup> <sup>&</sup>lt; *<sup>y</sup>* <sup>&</sup>lt; <sup>10</sup>

.

Fig. 3. Exact field **u** = (*x*, −*z*) in red, adjusted field obtained by the *CG–algorithm* in blue.

Mass–Consistent Wind Field Models: Numerical Techniques by L2–Projection Methods 31

The "exact" wind field satisfies ∇ · **<sup>u</sup>** <sup>=</sup> 1.2 <sup>×</sup> <sup>10</sup>−16. Figure 4 shows the adjusted and "exact"

Fig. 4. "Exact" field for cosine topography in red, adjusted field obtained by the *CG–algorithm*

where *h*(*x*) is a function constructed via cubic splines, which interpolate discrete data over 10 Km of real topography of a certain region in Mexico, contained in a database (GTOPO, 1997). The "exact" wind field satisfies ∇ · **<sup>u</sup>** <sup>=</sup> 6.1 <sup>×</sup> <sup>10</sup>−16. Figure 5 shows the adjusted and "exact" wind fields. We have an excellent agreement in all cases, even on truncated artificial boundaries. The relative error and the computed mean divergence are about the same order as in example 2. Table 2 shows a summary of the numerical results obtained with

(*x*, *<sup>y</sup>*) <sup>∈</sup> **<sup>R</sup>**<sup>2</sup> : 0 <sup>&</sup>lt; *<sup>x</sup>* <sup>&</sup>lt; 10, *<sup>h</sup>*(*x*) <sup>&</sup>lt; *<sup>y</sup>* <sup>&</sup>lt; <sup>10</sup>

**Example 4.** Terrain elevation from real data. In this case, the domain is defined as

**Example 3.** Cosine–shape topography. In this case, we define the domain as follows

(*x*, *<sup>y</sup>*) <sup>∈</sup> **<sup>R</sup>**<sup>2</sup> : 0 <sup>&</sup>lt; *<sup>x</sup>* <sup>&</sup>lt; 10, <sup>1</sup>

0 0.2 0.4 0.6 0.8 1 1.2

> Ω =

wind fields.

in blue.

Do *k* = *k* + 1 and return to 2.

Above, �·, ·� indicates the usual scalar product in *L*2(Ω). Observe that the adjusted field **u** is also computed as an intermediate step in the algorithm. In this algorithm, no boundary conditions are imposed on *λ*, contrary to what it was done in the first approach. This fact has a very important effect in the numerical calculation.

#### **4.3 A mixed finite element method**

To approximate the functions in **V***<sup>N</sup>* and *L*2(Ω), considered in the previous algorithm, we use the *Bercovier–Pironneau* finite element approximation (Bercovier & Pironneau, 1979). Functions in *L*2(Ω) are approximated by continuous piecewise linear polynomials over a triangulation T*<sup>h</sup>* of Ω, while the elements in **V***<sup>N</sup>* are also approximated by linear polynomials but now over a twice finer triangulation T*h*/2 of Ω. The fine triangulation T*h*/2 is obtained from a regular subdivision of each triangle *T* ∈ T*h*. Then, the functional spaces **V***<sup>N</sup>* and *L*2(Ω) will be approximated, respectively, by the finite dimensional spaces

$$\begin{aligned} \mathbf{V}\_{N\hbar} &= \left\{ \mathbf{v}\_{\hbar} \in \mathbb{C}^{0}(\bar{\Omega})^{2} \;:\; \mathbf{v}\_{\hbar}|\_{T} \in P\_{1} \times P\_{1\prime} \quad \forall \; T \in \mathcal{T}\_{\hbar/2\prime} \quad \mathbf{v}\_{\hbar} \cdot \mathbf{n} = 0 \text{ on } \Gamma\_{N} \right\}, \\ L\_{\hbar} &= \left\{ q\_{\hbar} \in \mathbb{C}^{0}(\bar{\Omega}) \;:\; q\_{\hbar}|\_{T} \in P\_{1\prime} \quad \forall \; T \in \mathcal{T}\_{\hbar} \right\}, \end{aligned}$$

We apply this mixed method, particularly to solve the integral equations in steps 1 and 2, as well as for the calculation of the weak divergence to obtain *g*<sup>0</sup> in step 1 and *w<sup>k</sup>* in step 2. Those calculations require this mixed method, or any other stable finite element pair, to avoid instabilities in the numerical solution. Actually, the main cost of this algorithm is the solution at each iteration of the integral equation to get **u***<sup>k</sup>* and the calculation of *wk*. However, if the trapezoidal rule is applied to approximate the left hand side of the integral equations, we obtain a system of algebraic equations with diagonal matrix, and the cost to solve them is just a vector multiplication. We call this new algorithm the *CG–algorithm*.

**Example 2**. We consider again the initial horizontal field **u**<sup>I</sup> = (*x*, 0), as in Example 1 to test the performance of the *CG–algorithm*. In order to compare the numerical results with those obtained with the *E1–algorithm*, we chose *h* = 1/40 and *h*/2 = 1/80 in this case. To stop the iterations we choose *ε* = 10−<sup>8</sup> at step 3. Figure 3 shows the exact and the adjusted wind fields. The agreement is excellent this time, even at the vertical boundaries *<sup>x</sup>* <sup>=</sup> 1 and *<sup>x</sup>* <sup>=</sup> 2. The relative error and the average divergence are *er* <sup>=</sup> 5.9 <sup>×</sup> <sup>10</sup>−<sup>4</sup> and *mdiv* <sup>=</sup> <sup>−</sup>5.3 <sup>×</sup> <sup>10</sup>−12, respectively. Note that we got a significant improvement: nearly two orders of magnitude better on the relative error, and about ten orders of magnitude better on the average divergence. The improvement of the relative error is mainly due to the reduction of the error on truncated boundaries, while the enhancing of average divergence is mainly due to the iterative method, because it stops when it reaches the tolerance (i.e. when the norm of the divergence is small enough).

To test further the *CG–algorithm* we consider two, more "realistic", additional examples. The first one includes a domain with a topography of a cosine–shape, and the second one includes a domain with a real topography. In both cases, the "exact" wind field was obtained with a Stokes solver using the methodology described in (Glowinski, 2003). The initial wind field **u***<sup>I</sup>* was obtained dropping the vertical component of the "exact" one in both cases. Then, the vector wind field is recovered using the same discretization parameters as in example 2.

8 Will-be-set-by-IN-TECH

Above, �·, ·� indicates the usual scalar product in *L*2(Ω). Observe that the adjusted field **u** is also computed as an intermediate step in the algorithm. In this algorithm, no boundary conditions are imposed on *λ*, contrary to what it was done in the first approach. This fact has

To approximate the functions in **V***<sup>N</sup>* and *L*2(Ω), considered in the previous algorithm, we use the *Bercovier–Pironneau* finite element approximation (Bercovier & Pironneau, 1979). Functions in *L*2(Ω) are approximated by continuous piecewise linear polynomials over a triangulation T*<sup>h</sup>* of Ω, while the elements in **V***<sup>N</sup>* are also approximated by linear polynomials but now over a twice finer triangulation T*h*/2 of Ω. The fine triangulation T*h*/2 is obtained from a regular subdivision of each triangle *T* ∈ T*h*. Then, the functional spaces **V***<sup>N</sup>* and *L*2(Ω)

We apply this mixed method, particularly to solve the integral equations in steps 1 and 2, as well as for the calculation of the weak divergence to obtain *g*<sup>0</sup> in step 1 and *w<sup>k</sup>* in step 2. Those calculations require this mixed method, or any other stable finite element pair, to avoid instabilities in the numerical solution. Actually, the main cost of this algorithm is the solution at each iteration of the integral equation to get **u***<sup>k</sup>* and the calculation of *wk*. However, if the trapezoidal rule is applied to approximate the left hand side of the integral equations, we obtain a system of algebraic equations with diagonal matrix, and the cost to solve them is just

**Example 2**. We consider again the initial horizontal field **u**<sup>I</sup> = (*x*, 0), as in Example 1 to test the performance of the *CG–algorithm*. In order to compare the numerical results with those obtained with the *E1–algorithm*, we chose *h* = 1/40 and *h*/2 = 1/80 in this case. To stop the iterations we choose *ε* = 10−<sup>8</sup> at step 3. Figure 3 shows the exact and the adjusted wind fields. The agreement is excellent this time, even at the vertical boundaries *<sup>x</sup>* <sup>=</sup> 1 and *<sup>x</sup>* <sup>=</sup> 2. The relative error and the average divergence are *er* <sup>=</sup> 5.9 <sup>×</sup> <sup>10</sup>−<sup>4</sup> and *mdiv* <sup>=</sup> <sup>−</sup>5.3 <sup>×</sup> <sup>10</sup>−12, respectively. Note that we got a significant improvement: nearly two orders of magnitude better on the relative error, and about ten orders of magnitude better on the average divergence. The improvement of the relative error is mainly due to the reduction of the error on truncated boundaries, while the enhancing of average divergence is mainly due to the iterative method, because it stops when it reaches the tolerance (i.e. when the norm

To test further the *CG–algorithm* we consider two, more "realistic", additional examples. The first one includes a domain with a topography of a cosine–shape, and the second one includes a domain with a real topography. In both cases, the "exact" wind field was obtained with a Stokes solver using the methodology described in (Glowinski, 2003). The initial wind field **u***<sup>I</sup>* was obtained dropping the vertical component of the "exact" one in both cases. Then, the vector wind field is recovered using the same discretization parameters as in example 2.

<sup>2</sup> : **<sup>v</sup>***h*|*<sup>T</sup>* <sup>∈</sup> *<sup>P</sup>*<sup>1</sup> <sup>×</sup> *<sup>P</sup>*1, <sup>∀</sup> *<sup>T</sup>* ∈ T*h*/2, **<sup>v</sup>***<sup>h</sup>* · **<sup>n</sup>** <sup>=</sup> 0 on <sup>Γ</sup>*<sup>N</sup>*

 ,  ,

Do *k* = *k* + 1 and return to 2.

**4.3 A mixed finite element method**

**V***Nh* =

*Lh* = 

of the divergence is small enough).

**<sup>v</sup>***<sup>h</sup>* <sup>∈</sup> *<sup>C</sup>*0(Ω¯ )

a very important effect in the numerical calculation.

will be approximated, respectively, by the finite dimensional spaces

*qh* <sup>∈</sup> *<sup>C</sup>*0(Ω¯ ) : *qh*|*<sup>T</sup>* <sup>∈</sup> *<sup>P</sup>*1, <sup>∀</sup> *<sup>T</sup>* ∈ T*<sup>h</sup>*

a vector multiplication. We call this new algorithm the *CG–algorithm*.

Fig. 3. Exact field **u** = (*x*, −*z*) in red, adjusted field obtained by the *CG–algorithm* in blue.

**Example 3.** Cosine–shape topography. In this case, we define the domain as follows

$$\Omega = \left\{ (\mathbf{x}, y) \in \mathbb{R}^2 \; : \; 0 < \mathbf{x} < 10, \; \frac{1}{2} \cos \frac{3\pi \mathbf{x}}{10} + \frac{1}{2} < y < 10 \right\}.$$

The "exact" wind field satisfies ∇ · **<sup>u</sup>** <sup>=</sup> 1.2 <sup>×</sup> <sup>10</sup>−16. Figure 4 shows the adjusted and "exact" wind fields.

Fig. 4. "Exact" field for cosine topography in red, adjusted field obtained by the *CG–algorithm* in blue.

**Example 4.** Terrain elevation from real data. In this case, the domain is defined as

$$\Omega = \left\{ (\mathbf{x}, \mathbf{y}) \in \mathbb{R}^2 \, : \, 0 < \mathbf{x} < 10, \, \, h(\mathbf{x}) < \mathbf{y} < 10 \right\},$$

where *h*(*x*) is a function constructed via cubic splines, which interpolate discrete data over 10 Km of real topography of a certain region in Mexico, contained in a database (GTOPO, 1997). The "exact" wind field satisfies ∇ · **<sup>u</sup>** <sup>=</sup> 6.1 <sup>×</sup> <sup>10</sup>−16. Figure 5 shows the adjusted and "exact" wind fields. We have an excellent agreement in all cases, even on truncated artificial boundaries. The relative error and the computed mean divergence are about the same order as in example 2. Table 2 shows a summary of the numerical results obtained with

Then, from (27)–(28) we obtain

**5. Some extensions and future research**

From equations (19)–(20), we obtain Ω 

and its corresponding elliptic problem (right):

*<sup>S</sup>***<sup>u</sup>** − ∇*<sup>λ</sup>* <sup>=</sup> *<sup>S</sup>***u**<sup>I</sup>

**5.1 Alternative boundary conditions for the elliptic problem**

*<sup>S</sup>***<sup>u</sup>** − ∇*<sup>λ</sup>* <sup>−</sup> *<sup>S</sup>***u**<sup>I</sup>

 Ω , and ∇ · **<sup>u</sup>** <sup>=</sup> 0 in <sup>Ω</sup>, −∇ ·

· **v** *d***x** =

 ΓΓ*<sup>N</sup>*

The boundary integral in (30) vanishes in two cases, namely: when **v** · **n** = 0 or when *λ* = 0 on Γ Γ*N*. The first case is not possible since it holds only on Γ*N*, and the second case is not a good choice on vertical boundaries as we have already seen in Section 3. However, there is a possibility: decompose Γ*<sup>D</sup>* as the union of the vertical boundaries, Γ*V*, and the top boundary, Γ*T*. Now, at Γ*<sup>T</sup>* we still impose *λ* = 0, and on Γ*<sup>V</sup>* we impose the new boundary condition **<sup>u</sup>** · **<sup>n</sup>** <sup>=</sup> **<sup>u</sup>**<sup>I</sup> · **<sup>n</sup>**. This new boundary condition is reasonable, since we assume that **<sup>u</sup>**<sup>I</sup> is the horizontal part of **u**. Therefore, with this choice, we obtain the saddle–point problem (left)

*λ* **v** · **n** *d***Γ**, ∀ **v** ∈ **V***N*, (30)

<sup>=</sup> ∇ · **<sup>u</sup>**<sup>I</sup> in <sup>Ω</sup>, (32)

*q*∇ · **u** *d***x** = 0, ∀ *q* ∈ *L*2(Ω). (31)

*<sup>S</sup>*−1∇*<sup>λ</sup>* 

*λ* = 0 on Γ*T*, *λ* = 0 on Γ*T*, (33) **<sup>u</sup>** · **<sup>n</sup>** <sup>=</sup> **<sup>u</sup>**<sup>I</sup> · **<sup>n</sup>** on <sup>Γ</sup>*V*, <sup>−</sup>*S*−1∇*<sup>λ</sup>* · **<sup>n</sup>** <sup>=</sup> 0 on <sup>Γ</sup>*V*, (34) **<sup>u</sup>** · **<sup>n</sup>** <sup>=</sup> 0 on <sup>Γ</sup>*N*. <sup>−</sup>*S*−1∇*<sup>λ</sup>* · **<sup>n</sup>** <sup>=</sup> **<sup>u</sup>**<sup>I</sup> · **<sup>n</sup>** on <sup>Γ</sup>*N*. (35)

*<sup>A</sup>* (*B q*) = *<sup>A</sup> <sup>φ</sup><sup>q</sup>* <sup>=</sup> −∇ · (*S*−1∇*φq*) = *<sup>q</sup>*. (29)

This shows that *B* can be used as an optimal preconditioner. Therefore, the additional cost of the preconditioned conjugate gradient algorithm is the solution of an elliptic problem at each iteration. However, this additional cost is offsetted by two nice properties: a) the preconditioning must reduce drastically the number of iterations (from about 1000 to less than 20, based on previous experience in CFD); b) there is a significant reduction of degrees of freedom in the discrete version of the elliptic problem associated to operator *B*. This elliptic problem is solved in a coarse mesh, and it is four times smaller than the elliptic problem for

Mass–Consistent Wind Field Models: Numerical Techniques by L2–Projection Methods 33

In this section, we present some additional alternatives to look at the problem. We first consider a different set of boundary conditions on vertical truncated boundaries for the multiplier *λ*, and we show that it produces better results than the traditional approach. We also show that if we introduce ghost nodes on the truncated artificial boundaries, we get even a better improvement. Finally, we introduce radial basis functions to solve the elliptic problems for the multiplier, and show that this is a promising alternative for 3-D wind fields.

the multiplier *λ* in the 2–D case, and about eight times smaller in 3–D problems.

the *CG-algorithm*. All numerical calculations were performed in a DELL Latitude D610 2.13 GHz laptop with an Intel Pentium M processor and 2 GB of RAM.

Fig. 5. "Exact" field for real topography in red, adjusted field obtained by the *CG–algorithm* in blue.


Table 2. Performance of the *CG–algorithm* for three different cases.

#### **4.4 Preconditioned conjugate gradient method**

The CPU time to solve the problem with the *CG–algorithm*, at the level of accuracy shown in Table 2, is about twice the CPU time needed to solve the problem with the *E1–algorithm*. In order to make this algorithm more reliable we need to speed up the iterative algorithm to get at least a comparable computational efficiency. Fortunately, we have found a good preconditioner for the iterative algorithm. This preconditioner is an optimal one, and we are presently working in its computer implementation, so we only describe here the main ideas without presenting numerical results yet.

Let *B* : *L*2(Ω) → *L*2(Ω) be an operator defined by

$$B\,\boldsymbol{q} = \phi\_{\boldsymbol{q}}, \quad \text{where} \quad \phi\_{\boldsymbol{q}} \text{ solves}: \int\_{\Omega} \left(\boldsymbol{\mathcal{S}}^{-1} \nabla \phi\_{\boldsymbol{q}}\right) \cdot \nabla \psi \,d\mathbf{x} = \int\_{\Omega} \boldsymbol{q} \,\psi \,d\mathbf{x} \qquad \forall \,\psi \in H^{1}(\Omega) \tag{26}$$

Operator *<sup>B</sup>* is self-adjoint and elliptic, and satisfies *<sup>A</sup>* (*B q*) = *<sup>q</sup>*, inside <sup>Ω</sup>, for every *<sup>q</sup>* <sup>∈</sup> *<sup>L</sup>*2(Ω). An easy way to see these properties is considering the differential form of operators *A* and *B*:

$$A\,\boldsymbol{q} = -\nabla \cdot \mathbf{u}\_{\boldsymbol{q}} = -\nabla \cdot (\mathbb{S}^{-1} \nabla \boldsymbol{q})\_{\prime} \quad \text{since} \quad \mathbb{S} \,\mathbf{u}\_{\boldsymbol{q}} = \nabla \boldsymbol{q} \quad \text{in} \quad \Omega\_{\prime} \tag{27}$$

$$B\,q = \phi\_q = -[\nabla \cdot (\mathbb{S}^{-1} \nabla)]^{-1} q\_\prime \qquad \text{since} \quad -\nabla \cdot (\mathbb{S}^{-1} \nabla \phi\_q) = q \quad \text{in} \quad \Omega. \tag{28}$$

Then, from (27)–(28) we obtain

10 Will-be-set-by-IN-TECH

the *CG-algorithm*. All numerical calculations were performed in a DELL Latitude D610 2.13

Fig. 5. "Exact" field for real topography in red, adjusted field obtained by the *CG–algorithm* in

Example Case with *er mdiv* No. iters. CPU–time (sec.) 2 u(x,z) = (x,-z) 5.9 <sup>×</sup> <sup>10</sup>−<sup>4</sup> <sup>−</sup>5.3 <sup>×</sup> <sup>10</sup>−<sup>12</sup> <sup>1214</sup> 3.9 3 cosine topography 3.6 <sup>×</sup> <sup>10</sup>−<sup>6</sup> 9.8 <sup>×</sup> <sup>10</sup>−<sup>9</sup> <sup>955</sup> 3.3 4 real topography 4.1 <sup>×</sup> <sup>10</sup>−<sup>6</sup> 5.7 <sup>×</sup> <sup>10</sup>−<sup>11</sup> <sup>1000</sup> 3.7

The CPU time to solve the problem with the *CG–algorithm*, at the level of accuracy shown in Table 2, is about twice the CPU time needed to solve the problem with the *E1–algorithm*. In order to make this algorithm more reliable we need to speed up the iterative algorithm to get at least a comparable computational efficiency. Fortunately, we have found a good preconditioner for the iterative algorithm. This preconditioner is an optimal one, and we are presently working in its computer implementation, so we only describe here the main ideas

8 8.5 9 9.5 10

8 8.5 9 9.5 10

GHz laptop with an Intel Pentium M processor and 2 GB of RAM.

<sup>0</sup> <sup>2</sup> <sup>4</sup> <sup>6</sup> <sup>8</sup> <sup>10</sup> <sup>12</sup> <sup>1</sup>

Table 2. Performance of the *CG–algorithm* for three different cases.

 Ω 

*<sup>S</sup>*−1∇*φ<sup>q</sup>*

Operator *<sup>B</sup>* is self-adjoint and elliptic, and satisfies *<sup>A</sup>* (*B q*) = *<sup>q</sup>*, inside <sup>Ω</sup>, for every *<sup>q</sup>* <sup>∈</sup> *<sup>L</sup>*2(Ω). An easy way to see these properties is considering the differential form of operators *A* and *B*:

· ∇*ψ d***x** =

*A q* <sup>=</sup> −∇·**u***<sup>q</sup>* <sup>=</sup> −∇ · (*S*−1<sup>∇</sup> *<sup>q</sup>*), since *<sup>S</sup>* **<sup>u</sup>***<sup>q</sup>* <sup>=</sup> <sup>∇</sup>*<sup>q</sup>* in <sup>Ω</sup>, (27) *B q* <sup>=</sup> *<sup>φ</sup><sup>q</sup>* <sup>=</sup> <sup>−</sup>[∇ · (*S*−1∇)]−<sup>1</sup> *<sup>q</sup>*, since −∇· (*S*−1∇*φq*) = *<sup>q</sup>* in <sup>Ω</sup>. (28)

 Ω

*<sup>q</sup> <sup>ψ</sup> <sup>d</sup>***<sup>x</sup>** <sup>∀</sup> *<sup>ψ</sup>* <sup>∈</sup> *<sup>H</sup>*1(Ω) (26)

**4.4 Preconditioned conjugate gradient method**

without presenting numerical results yet.

*B q* = *φq*, where *φ<sup>q</sup>* solves :

Let *B* : *L*2(Ω) → *L*2(Ω) be an operator defined by

blue.

$$A\left(B\,\boldsymbol{q}\right) = A\,\phi\_{\boldsymbol{q}} = -\nabla \cdot \left(\boldsymbol{S}^{-1}\nabla\phi\_{\boldsymbol{q}}\right) = \boldsymbol{q}.\tag{29}$$

This shows that *B* can be used as an optimal preconditioner. Therefore, the additional cost of the preconditioned conjugate gradient algorithm is the solution of an elliptic problem at each iteration. However, this additional cost is offsetted by two nice properties: a) the preconditioning must reduce drastically the number of iterations (from about 1000 to less than 20, based on previous experience in CFD); b) there is a significant reduction of degrees of freedom in the discrete version of the elliptic problem associated to operator *B*. This elliptic problem is solved in a coarse mesh, and it is four times smaller than the elliptic problem for the multiplier *λ* in the 2–D case, and about eight times smaller in 3–D problems.

#### **5. Some extensions and future research**

In this section, we present some additional alternatives to look at the problem. We first consider a different set of boundary conditions on vertical truncated boundaries for the multiplier *λ*, and we show that it produces better results than the traditional approach. We also show that if we introduce ghost nodes on the truncated artificial boundaries, we get even a better improvement. Finally, we introduce radial basis functions to solve the elliptic problems for the multiplier, and show that this is a promising alternative for 3-D wind fields.

#### **5.1 Alternative boundary conditions for the elliptic problem**

From equations (19)–(20), we obtain

$$\int\_{\Omega} \left( \mathbf{S} \mathbf{u} - \nabla \lambda - \mathbf{S} \mathbf{u}^{\mathrm{I}} \right) \cdot \mathbf{v} \, d\mathbf{x} = \int\_{\Gamma \times \Gamma\_{N}} \lambda \, \mathbf{v} \cdot \mathbf{n} \, d\Gamma\_{\prime} \quad \forall \, \mathbf{v} \in \mathbf{V}\_{N\prime} \tag{30}$$

$$\int\_{\Omega} q \nabla \cdot \mathbf{u} \, d\mathbf{x} = 0, \quad \forall q \in L\_2(\Omega). \tag{31}$$

The boundary integral in (30) vanishes in two cases, namely: when **v** · **n** = 0 or when *λ* = 0 on Γ Γ*N*. The first case is not possible since it holds only on Γ*N*, and the second case is not a good choice on vertical boundaries as we have already seen in Section 3. However, there is a possibility: decompose Γ*<sup>D</sup>* as the union of the vertical boundaries, Γ*V*, and the top boundary, Γ*T*. Now, at Γ*<sup>T</sup>* we still impose *λ* = 0, and on Γ*<sup>V</sup>* we impose the new boundary condition **<sup>u</sup>** · **<sup>n</sup>** <sup>=</sup> **<sup>u</sup>**<sup>I</sup> · **<sup>n</sup>**. This new boundary condition is reasonable, since we assume that **<sup>u</sup>**<sup>I</sup> is the horizontal part of **u**. Therefore, with this choice, we obtain the saddle–point problem (left) and its corresponding elliptic problem (right):

$$\mathbf{S} \cdot \mathbf{S} \mathbf{u} - \nabla \lambda = \mathbf{S} \mathbf{u}^{\mathrm{I}}, \quad \text{and} \quad \nabla \cdot \mathbf{u} = 0 \quad \text{in } \Omega, \qquad -\nabla \cdot \left( \mathbf{S}^{-1} \nabla \lambda \right) = \nabla \cdot \mathbf{u}^{\mathrm{I}} \quad \text{in } \Omega, \tag{32}$$

$$
\lambda = 0 \quad \text{on } \Gamma\_{\overline{T}\prime} \tag{33}
$$

$$\mathbf{u} \cdot \mathbf{n} = \mathbf{u}^{\mathrm{I}} \cdot \mathbf{n} \quad \text{on } \Gamma\_{V\_{\prime}} \qquad \qquad \qquad \qquad -\mathbb{S}^{-1} \nabla \lambda \cdot \mathbf{n} = 0 \quad \text{on} \quad \Gamma\_{V\_{\prime}} \qquad \tag{34}$$

$$\mathbf{u} \cdot \mathbf{n} = 0 \quad \text{on } \Gamma\_N. \qquad \qquad \qquad -\mathbf{S}^{-1} \nabla \lambda \cdot \mathbf{n} = \mathbf{u}^\mathrm{I} \cdot \mathbf{n} \quad \text{on} \quad \Gamma\_N. \tag{35}$$

Instead of looking for boundary conditions on Γ*D*, we may enforce mass conservation by

Mass–Consistent Wind Field Models: Numerical Techniques by L2–Projection Methods 35

Ω

 Γ*<sup>D</sup>*

the compatibility condition associated to the the above Poisson–Neumann–like problem. Therefore, this problem has a unique solution *<sup>λ</sup>* <sup>∈</sup> *<sup>H</sup>*1(Ω)/**R**, and its variational formulation

Observe that when *q* = 1 in *H*1(Ω)/**R**, we recover (39). However, the computational solution of this problem is not trivial, since the matrix associated to the discrete version is semidefinite. On the other hand, the symmetry of the matrix is lost because the boundary integral in the right–hand side has the unknown *λ*. A way to overcome this computational problem is to introduce "ghost nodes", around and beyond of the nonphysical truncated boundary Γ*D*. Then, we may impose *λ<sup>h</sup>* = 0 and/or *∂λh*/*∂***n** = 0 on the outer layer of those ghost nodes. At the end, we discard the solution on the ghost nodes, and we only keep the solution values on the actual nodes. Actually, this is a well–known way to deal with differential equations in

**Example 6.** We consider one more time the problem from example 1 with the same discretization parameters. We incorporate two layers of **ghost nodes** and impose *λ* = 0 on the outer layer. The recovered wind field obtained is such that the relative error and average weak divergence are *er* <sup>=</sup> 2.1 <sup>×</sup> <sup>10</sup>−<sup>5</sup> and *mdiv* <sup>=</sup> 1.6 <sup>×</sup> <sup>10</sup>−6, respectively. The figure with the comparison of the adjusted wind field and the exact wind field is not shown, because it is very similar to Figure 5. Instead, we summarize in Table 4 the results for this example with

> Case *E1–algortihm E2–algortithm Ghost-Nodes CG–algorithm er* 1.9 <sup>×</sup> <sup>10</sup>−<sup>2</sup> 4.0 <sup>×</sup> <sup>10</sup>−<sup>4</sup> 2.1 <sup>×</sup> <sup>10</sup>−<sup>5</sup> 5.4 <sup>×</sup> <sup>10</sup>−<sup>4</sup> *mdiv* 4.1 <sup>×</sup> <sup>10</sup>−<sup>2</sup> 1.8 <sup>×</sup> <sup>10</sup>−<sup>2</sup> 1.6 <sup>×</sup> <sup>10</sup>−<sup>6</sup> <sup>−</sup>5.2 <sup>×</sup> <sup>10</sup>−<sup>12</sup>

Table 4 shows how boundary conditions degrade numerical calculations. It is observed that the solution improves each time the Dirichlet boundary condition *λ* = 0 is applied to a smaller section of the non-physical boundary. This is not surprising, since this boundary condition introduces a large artificial gradient, mainly on vertical truncated boundaries, when calculating the term <sup>∇</sup>*<sup>λ</sup>* in order to get **<sup>u</sup>** <sup>=</sup> **<sup>u</sup>**<sup>I</sup> <sup>+</sup> *<sup>S</sup>*−1∇*<sup>λ</sup>* at the corresponding boundary nodes. At this time, and taking in account the performance of every algorithm, we may recommend to use either the classical approach with ghost nodes or the saddle point problem approach with the conjugate gradient algorithm, specially if we do not have enough information at

Table 4. Comparison of numerical solutions obtained with different algorithms.

∇ · **<sup>u</sup>**<sup>I</sup> *<sup>d</sup>***<sup>x</sup>** <sup>=</sup>

(**u**<sup>I</sup> <sup>+</sup> *<sup>S</sup>*−1∇*λ*) · **<sup>n</sup>** *<sup>d</sup>*<sup>Γ</sup> <sup>=</sup> 0. (39)

*<sup>q</sup>* (**u**<sup>I</sup> <sup>+</sup> *<sup>S</sup>*−1∇*λ*) · **<sup>n</sup>** *<sup>d</sup>*<sup>Γ</sup> , <sup>∀</sup> *<sup>q</sup>* <sup>∈</sup> *<sup>H</sup>*1(Ω)/**R**. (40)

**<sup>u</sup>**<sup>I</sup> · **<sup>n</sup>** *<sup>d</sup>*Γ. Actually, this is

 Γ

asking any solution of (37)–(38) to satisfy

Equations (37)–(39) imply the identity

(*S*−1∇*λ*) · ∇*q d***<sup>x</sup>** <sup>=</sup> <sup>−</sup>

domains with truncated boundaries.

the different algorithms.

truncated boundaries.

 Ω  Γ

is: Given **<sup>u</sup>**<sup>I</sup> <sup>∈</sup> **<sup>L</sup>**2(Ω), find *<sup>λ</sup>* <sup>∈</sup> *<sup>H</sup>*1(Ω)/**<sup>R</sup>** such that

 Ω

**u** · **n** *d*Γ =

**<sup>u</sup>**<sup>I</sup> · ∇*q d***<sup>x</sup>** <sup>+</sup>

 Γ*<sup>D</sup>*

The finite element algorithm for the elliptic problem is: Given **u**<sup>I</sup> *<sup>h</sup>* ∈ **L***h*, find *λh*∈*Hh* such that

$$\int\_{\Omega} \mathbf{S}^{-1} \nabla \lambda\_{\hbar} \cdot \nabla q \, d\mathbf{x} = -\int\_{\Omega} \mathbf{u}\_{\hbar}^{\mathrm{I}} \cdot \nabla q \, d\mathbf{x} + \int\_{\Gamma\_{V}} q \, \mathbf{u}\_{\hbar}^{\mathrm{I}} \cdot \mathbf{n} \, d\Gamma\_{\prime} \quad \forall q \in H\_{h\nu} \tag{36}$$

where *Hh* is defined as in (13), but with *q* = 0 on Γ*<sup>T</sup>* instead of Γ*D*. Equation (36) differs from equation (14) only by the boundary integral on Γ*V*. We call (36), together with (15), the *E2–algorithm*.

**Example 5.** Let us consider one more time the problem introduced in example 1 and in example 3, with the same discretization parameter, *h* = 1/80. The recovered wind field for both cases is better than the one obtained with the *E1–algorithm*, since the vertical component is recovered fairly well, not only in the interior of the domain but also at the vertical boundaries. Figure 6 shows the exact and recovered wind fields. Table 3 shows a summary of the results obtained in examples 1, 2, 3 and 5. For the case with exact wind field **u** = (*x*, −*z*) the immediate effect of this improvement is the reduction of the relative error by two orders of magnitude. However, we do not obtain a comparable reduction of the mean divergence. For the problems with cosine topography occurs the opposite. Actually, the numerical results show that the most effective algorithm to reduce both, the relative error and the average divergence is the *CG–algorithm*.

Fig. 6. Exact field (red) and adjusted field obtained by the *E2–algorithm* (blue). Left: case with exact field **u** = (*x*, −*z*). Right: case with cosine topography.


Table 3. Summary of the numerical results obtained in examples 1, 2, 3 and 5.

#### **5.2 Ghost nodes**

Given that **u** belongs to **V** and satisfies (2), then *λ* satisfies the equations

$$-\nabla \cdot (\mathbf{S}^{-1} \nabla \lambda) = \nabla \cdot \mathbf{u}^{\mathrm{I}} , \quad \text{in} \quad \Omega \,, \tag{37}$$

$$-(\mathbb{S}^{-1}\nabla\lambda)\cdot\mathbf{n} = \mathbf{u}^{\mathrm{I}}\cdot\mathbf{n}, \quad \text{on} \quad \Gamma\_N. \tag{38}$$

12 Will-be-set-by-IN-TECH

*<sup>h</sup>* · ∇*q d***x** +

where *Hh* is defined as in (13), but with *q* = 0 on Γ*<sup>T</sup>* instead of Γ*D*. Equation (36) differs from equation (14) only by the boundary integral on Γ*V*. We call (36), together with (15), the

**Example 5.** Let us consider one more time the problem introduced in example 1 and in example 3, with the same discretization parameter, *h* = 1/80. The recovered wind field for both cases is better than the one obtained with the *E1–algorithm*, since the vertical component is recovered fairly well, not only in the interior of the domain but also at the vertical boundaries. Figure 6 shows the exact and recovered wind fields. Table 3 shows a summary of the results obtained in examples 1, 2, 3 and 5. For the case with exact wind field **u** = (*x*, −*z*) the immediate effect of this improvement is the reduction of the relative error by two orders of magnitude. However, we do not obtain a comparable reduction of the mean divergence. For the problems with cosine topography occurs the opposite. Actually, the numerical results show that the most effective algorithm to reduce both, the relative error and the average

 Γ*<sup>V</sup> q* **u**<sup>I</sup>

<sup>0</sup> <sup>2</sup> <sup>4</sup> <sup>6</sup> <sup>8</sup> <sup>10</sup> <sup>12</sup> -2

Fig. 6. Exact field (red) and adjusted field obtained by the *E2–algorithm* (blue). Left: case with

**Ex. Case with Algorithm** *er mdiv* **No. iters. CPU time(s)** <sup>1</sup> **<sup>u</sup>** = (*x*, <sup>−</sup>*z*) *E1–algorithm* 1.9 <sup>×</sup> <sup>10</sup>−<sup>2</sup> 4.1 <sup>×</sup> <sup>10</sup>−<sup>2</sup> — 1.78 <sup>2</sup> **<sup>u</sup>** = (*x*, <sup>−</sup>*z*) *CG–algorithm* 5.9 <sup>×</sup> <sup>10</sup>−<sup>4</sup> <sup>−</sup>5.3 <sup>×</sup> <sup>10</sup>−<sup>12</sup> 1214 3.90 <sup>5</sup> **<sup>u</sup>** = (*x*, <sup>−</sup>*z*) *E2–algorithm* 4.0 <sup>×</sup> <sup>10</sup>−<sup>4</sup> 1.8 <sup>×</sup> <sup>10</sup>−<sup>2</sup> — 1.78 3 cosine topography *CG–algorithm* 3.6 <sup>×</sup> <sup>10</sup>−<sup>6</sup> 9.8 <sup>×</sup> <sup>10</sup>−<sup>9</sup> 955 3.30 5 cosine topography *E2–algorithm* 9.2 <sup>×</sup> <sup>10</sup>−<sup>2</sup> 3.6 <sup>×</sup> <sup>10</sup>−<sup>4</sup> — 2.08

*<sup>h</sup>* ∈ **L***h*, find *λh*∈*Hh* such that

*<sup>h</sup>* · **n** *d*Γ, ∀ *q* ∈ *Hh*, (36)

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8

8 8.2 8.4 8.6 8.8 9 9.2 9.4 9.6 9.8 10

, in Ω, (37)

<sup>−</sup>(*S*−1∇*λ*) · **<sup>n</sup>** <sup>=</sup> **<sup>u</sup>**<sup>I</sup> · **<sup>n</sup>**, on <sup>Γ</sup>*N*. (38)

The finite element algorithm for the elliptic problem is: Given **u**<sup>I</sup>

 Ω **u**I

1.9 1.95 <sup>2</sup> 2.05 0.8

Table 3. Summary of the numerical results obtained in examples 1, 2, 3 and 5.

−∇ · (*S*−1∇*λ*) = ∇ · **<sup>u</sup>**<sup>I</sup>

Given that **u** belongs to **V** and satisfies (2), then *λ* satisfies the equations

0.85 0.9 0.95 1

exact field **u** = (*x*, −*z*). Right: case with cosine topography.

*<sup>S</sup>*−1∇*λ<sup>h</sup>* · ∇*q d***<sup>x</sup>** <sup>=</sup> <sup>−</sup>

 Ω

divergence is the *CG–algorithm*.

<sup>1</sup> 1.2 1.4 1.6 1.8 <sup>2</sup> 2.2 2.4 -0.2

*E2–algorithm*.

0 0.2 0.4 0.6 0.8 1 1.2

**5.2 Ghost nodes**

Instead of looking for boundary conditions on Γ*D*, we may enforce mass conservation by asking any solution of (37)–(38) to satisfy

$$\int\_{\Gamma} \mathbf{u} \cdot \mathbf{n} \, d\Gamma = \int\_{\Gamma\_D} (\mathbf{u}^\mathrm{I} + \mathcal{S}^{-1} \nabla \lambda) \cdot \mathbf{n} \, d\Gamma = 0. \tag{39}$$

Equations (37)–(39) imply the identity Ω ∇ · **<sup>u</sup>**<sup>I</sup> *<sup>d</sup>***<sup>x</sup>** <sup>=</sup> Γ **<sup>u</sup>**<sup>I</sup> · **<sup>n</sup>** *<sup>d</sup>*Γ. Actually, this is the compatibility condition associated to the the above Poisson–Neumann–like problem. Therefore, this problem has a unique solution *<sup>λ</sup>* <sup>∈</sup> *<sup>H</sup>*1(Ω)/**R**, and its variational formulation is: Given **<sup>u</sup>**<sup>I</sup> <sup>∈</sup> **<sup>L</sup>**2(Ω), find *<sup>λ</sup>* <sup>∈</sup> *<sup>H</sup>*1(Ω)/**<sup>R</sup>** such that

$$\int\_{\Omega} (\mathbb{S}^{-1} \nabla \lambda) \cdot \nabla q \, d\mathbf{x} = -\int\_{\Omega} \mathbf{u}^{\mathrm{I}} \cdot \nabla q \, d\mathbf{x} + \int\_{\Gamma\_D} q \left( \mathbf{u}^{\mathrm{I}} + \mathbb{S}^{-1} \nabla \lambda \right) \cdot \mathbf{n} \, d\Gamma, \quad \forall q \in H^{1}(\Omega) / \mathbb{R}. \tag{40}$$

Observe that when *q* = 1 in *H*1(Ω)/**R**, we recover (39). However, the computational solution of this problem is not trivial, since the matrix associated to the discrete version is semidefinite. On the other hand, the symmetry of the matrix is lost because the boundary integral in the right–hand side has the unknown *λ*. A way to overcome this computational problem is to introduce "ghost nodes", around and beyond of the nonphysical truncated boundary Γ*D*. Then, we may impose *λ<sup>h</sup>* = 0 and/or *∂λh*/*∂***n** = 0 on the outer layer of those ghost nodes. At the end, we discard the solution on the ghost nodes, and we only keep the solution values on the actual nodes. Actually, this is a well–known way to deal with differential equations in domains with truncated boundaries.

**Example 6.** We consider one more time the problem from example 1 with the same discretization parameters. We incorporate two layers of **ghost nodes** and impose *λ* = 0 on the outer layer. The recovered wind field obtained is such that the relative error and average weak divergence are *er* <sup>=</sup> 2.1 <sup>×</sup> <sup>10</sup>−<sup>5</sup> and *mdiv* <sup>=</sup> 1.6 <sup>×</sup> <sup>10</sup>−6, respectively. The figure with the comparison of the adjusted wind field and the exact wind field is not shown, because it is very similar to Figure 5. Instead, we summarize in Table 4 the results for this example with the different algorithms.


Table 4. Comparison of numerical solutions obtained with different algorithms.

Table 4 shows how boundary conditions degrade numerical calculations. It is observed that the solution improves each time the Dirichlet boundary condition *λ* = 0 is applied to a smaller section of the non-physical boundary. This is not surprising, since this boundary condition introduces a large artificial gradient, mainly on vertical truncated boundaries, when calculating the term <sup>∇</sup>*<sup>λ</sup>* in order to get **<sup>u</sup>** <sup>=</sup> **<sup>u</sup>**<sup>I</sup> <sup>+</sup> *<sup>S</sup>*−1∇*<sup>λ</sup>* at the corresponding boundary nodes. At this time, and taking in account the performance of every algorithm, we may recommend to use either the classical approach with ghost nodes or the saddle point problem approach with the conjugate gradient algorithm, specially if we do not have enough information at truncated boundaries.

∑*n*

0.2

**6. Conclusions**

0.4

z

0.6

0.8

1

*<sup>j</sup>*=<sup>1</sup> *<sup>λ</sup><sup>j</sup> <sup>φ</sup>*(�**<sup>x</sup>** <sup>−</sup> **<sup>x</sup>***j*�), where the unknown vector {*λi*}*<sup>n</sup>*

*n* ∑ *j*=1

*n* ∑ *j*=1

**u***h*(**x**) = **u**<sup>I</sup>

(**x**) +

*<sup>λ</sup><sup>j</sup>* <sup>L</sup> *<sup>φ</sup>*(�**x***<sup>i</sup>* <sup>−</sup> **<sup>x</sup>***j*�) = ∇·**u**<sup>I</sup>

Mass–Consistent Wind Field Models: Numerical Techniques by L2–Projection Methods 37

*n* ∑ *j*=1

obtained a relative error *er* <sup>=</sup> 1.98 <sup>×</sup> <sup>10</sup>−<sup>4</sup> and mean divergence *mdiv* <sup>=</sup> <sup>−</sup>5.59 <sup>×</sup> <sup>10</sup>−6.

1 1.5

Fig. 7. Collocation points (left), and comparison of exact and recovered wind fields (right).

We studied the problem of generating an adjusted wind field from horizontal wind data by different numerical techniques. We have shown that boundary conditions can significantly affect numerical solutions depending of how we treat artificial truncated boundaries. The usual methodology (*E1–algorithm*) does not produce satisfactory results close to the vertical boundaries due to the high gradients introduced by the term *<sup>S</sup>*−1∇*<sup>λ</sup>* in (2), when homogeneous Dirichlet boundary conditions are imposed there. The formulation of the problem as a saddle–point one, with a functional equation that has a self–adjoint and strongly elliptic operator, allows the use of an iterative conjugate gradient algorithm (the *CG–algorithm*). This new methodology, in the context of mass consistent models, produces

0

0.2

0.4

0.6

0.8

1

**Example 8**. As a last example we include a 3–D numerical calculation using radial basis functions. The exact syntectic wind field for this example is **u** = (*x*, *y*, −2*z*). We dropped the vertical component so that **u**<sup>I</sup> = (*x*, *y*, 0). A multiquadric kernel with *c* = 11.33 was used, and we chose *α*<sup>1</sup> = *α*<sup>2</sup> = *α*<sup>3</sup> = 1. The collocation points were obtained by a 5 × 5 × 5 regular subdivision of Ω = (1, 2) × (0, 1) × (0, 1). Figure 7 shows the collocation points and the comparison between the exact and recovered wind field. The agreement is excellent, we

L *λh*(**x***i*) =

B *λh*(**x***i*) =

Therefore the recovered wind field **u***<sup>h</sup>* is given by

<sup>1</sup> 1.2 1.4 1.6 1.8 <sup>2</sup>

<sup>2</sup> <sup>0</sup>

*<sup>i</sup>*=<sup>1</sup> satisfies the system of equations

*<sup>λ</sup><sup>j</sup> <sup>S</sup>*−1∇*φ*(�**<sup>x</sup>** <sup>−</sup> **<sup>x</sup>***j*�). (47)

<sup>1</sup> 1.2 1.4 1.6 1.8 <sup>2</sup>

1 1.5 2

*λ<sup>j</sup>* B *φ*(�**x***<sup>i</sup>* − **x***j*�) = *g*(**x***i*), *i* = *ni* + 1, . . . , *n*. (46)

(**x***i*), *i* = 1, . . . , *ni*, (45)

#### **5.3 Approximation with radial basis functions**

A function <sup>Φ</sup> : **<sup>R</sup>***<sup>d</sup>* <sup>→</sup> **<sup>R</sup>** is called a radial basis function (RBF) if <sup>Φ</sup>(**x**) = *<sup>φ</sup>*(�**x**�) where the kernel *<sup>φ</sup>* is a scalar function *<sup>φ</sup>* : **<sup>R</sup>**<sup>+</sup> <sup>→</sup> **<sup>R</sup>**, and *<sup>d</sup>* the spatial dimension. Usually �**x**� is denoted by *r*, and typical functions used in applications are, among others:


.

The constant value *c* es called the *shape parameter*. The radial basis function method was first introduced in the 1970s for multivariate scattered data approximation, (Hardy, 1971). This interpolation problem is defined as follows:

Given a set of points {**x***j*}*<sup>n</sup> <sup>j</sup>*=<sup>1</sup> <sup>⊂</sup> <sup>Ω</sup> <sup>⊂</sup> **<sup>R</sup>***d*, approximate the function *<sup>f</sup>*(**x**) from the set of values *fj* = *f*(**x***j*). A simple form is to define

$$s(\mathbf{x}) = \sum\_{j=1}^{k} \lambda\_j \phi(||\mathbf{x} - \mathbf{x}\_j||) + p(\mathbf{x}).\tag{41}$$

where *p* is a polynomial which depends on the specific RBF. Then the interpolation condition

$$s(\mathbf{x}\_i) = \sum\_{j=1}^k \lambda\_j \phi(||\mathbf{x}\_i - \mathbf{x}\_j||) + p(\mathbf{x}\_i) = f\_{i\prime} \quad \text{i} = 1, \ldots, n,\tag{42}$$

gives an algebraic system of equations for *<sup>λ</sup>* <sup>=</sup> {*λi*}*<sup>k</sup> <sup>j</sup>*=1. However, the corresponding matrix could be ill conditioned and, in some cases, even rank deficient, and special techniques are needed, like preconditioning and least squares (Buhmann, 2003), (Wendland, 2005).

In the last two decades, the main focus of the applications seems to have slowly shifted from scattered data approximation to the numerical solution of PDE. Radial basis function collocation methods for solving PDE are truly meshfree algorithms, in the sense that collocation points can be chosen freely and no connectivity between the points is needed or used (Kansa, 1990), (Narcowich & Ward, 1994). The main attraction of RBF collocation method to solve PDE is that it can be extended directly to solve 3–D problems. Moreover, due to the absence of a grid, these techniques are better suited than classical methods to cope with problems having complex boundaries. So, we think that the RBF collocation method is a good choice to study our problem, and we want to explore this alternative.

Let us denote by L the linear elliptic differential operator for the multiplier *λ*, and B the boundary operator. Suppose that we want to solve the problem

$$\mathcal{L}\,\lambda = -\nabla \cdot (\mathbb{S}^{-1}\nabla \lambda) = \nabla \cdot \mathbf{u}^{\mathrm{I}} \quad \text{in} \quad \Omega,\tag{43}$$

$$
\mathcal{B}\,\lambda = \emptyset \quad \text{on} \quad \Gamma\_N. \tag{44}
$$

We consider a set of collocation points {**x***i*}*<sup>n</sup> <sup>i</sup>*=1, with *ni* points in the interior and *nb* points on the boundary, so that *n* = *ni* + *n <sup>f</sup>* . We look for an approximate solution *λh*(**x**) = ∑*n <sup>j</sup>*=<sup>1</sup> *<sup>λ</sup><sup>j</sup> <sup>φ</sup>*(�**<sup>x</sup>** <sup>−</sup> **<sup>x</sup>***j*�), where the unknown vector {*λi*}*<sup>n</sup> <sup>i</sup>*=<sup>1</sup> satisfies the system of equations

$$\mathcal{L}\,\lambda\_{\rm li}(\mathbf{x}\_{i}) = \sum\_{j=1}^{n} \lambda\_{j}\,\mathcal{L}\,\phi(\|\mathbf{x}\_{i} - \mathbf{x}\_{j}\|) = \nabla \cdot \mathbf{u}^{\rm I}(\mathbf{x}\_{i}), \quad i = 1, \ldots, n\_{i} \tag{45}$$

$$\mathcal{B}\,\lambda\_{\mathbb{H}}(\mathbf{x}\_{i}) = \sum\_{j=1}^{n} \lambda\_{j}\mathcal{B}\,\phi(\|\mathbf{x}\_{i} - \mathbf{x}\_{j}\|) = \mathcal{g}(\mathbf{x}\_{i}), \quad i = n\_{i} + 1, \ldots, n. \tag{46}$$

Therefore the recovered wind field **u***<sup>h</sup>* is given by

14 Will-be-set-by-IN-TECH

A function <sup>Φ</sup> : **<sup>R</sup>***<sup>d</sup>* <sup>→</sup> **<sup>R</sup>** is called a radial basis function (RBF) if <sup>Φ</sup>(**x**) = *<sup>φ</sup>*(�**x**�) where the kernel *<sup>φ</sup>* is a scalar function *<sup>φ</sup>* : **<sup>R</sup>**<sup>+</sup> <sup>→</sup> **<sup>R</sup>**, and *<sup>d</sup>* the spatial dimension. Usually �**x**� is denoted

The constant value *c* es called the *shape parameter*. The radial basis function method was first introduced in the 1970s for multivariate scattered data approximation, (Hardy, 1971). This

where *p* is a polynomial which depends on the specific RBF. Then the interpolation condition

could be ill conditioned and, in some cases, even rank deficient, and special techniques are

In the last two decades, the main focus of the applications seems to have slowly shifted from scattered data approximation to the numerical solution of PDE. Radial basis function collocation methods for solving PDE are truly meshfree algorithms, in the sense that collocation points can be chosen freely and no connectivity between the points is needed or used (Kansa, 1990), (Narcowich & Ward, 1994). The main attraction of RBF collocation method to solve PDE is that it can be extended directly to solve 3–D problems. Moreover, due to the absence of a grid, these techniques are better suited than classical methods to cope with problems having complex boundaries. So, we think that the RBF collocation method is a good

Let us denote by L the linear elliptic differential operator for the multiplier *λ*, and B the

on the boundary, so that *n* = *ni* + *n <sup>f</sup>* . We look for an approximate solution *λh*(**x**) =

needed, like preconditioning and least squares (Buhmann, 2003), (Wendland, 2005).

*<sup>j</sup>*=<sup>1</sup> <sup>⊂</sup> <sup>Ω</sup> <sup>⊂</sup> **<sup>R</sup>***d*, approximate the function *<sup>f</sup>*(**x**) from the set of values

*λ<sup>j</sup> φ*(�**x** − **x***j*�) + *p*(**x**). (41)

*<sup>j</sup>*=1. However, the corresponding matrix

*λ<sup>j</sup> φ*(�**x***<sup>i</sup>* − **x***j*�) + *p*(**x***i*) = *fi*, *i* = 1, . . . , *n*, (42)

<sup>L</sup> *<sup>λ</sup>* <sup>=</sup> −∇ · (*S*−1∇*λ*) = ∇ · **<sup>u</sup>**<sup>I</sup> in <sup>Ω</sup>, (43) B *λ* = *g* on Γ*N*. (44)

*<sup>i</sup>*=1, with *ni* points in the interior and *nb* points

**5.3 Approximation with radial basis functions**

1. Multiquadrics, *<sup>φ</sup>*(*r*) = <sup>√</sup>

2. Gaussians, *φ*(*r*) = *e*−*c r*<sup>2</sup>

Given a set of points {**x***j*}*<sup>n</sup>*

3. Thin plate splines, *φ*(*r*) = *r*<sup>2</sup> ln(*r*).

4. Inverse multiquadrics, *<sup>φ</sup>*(*r*) = 1/√*c*<sup>2</sup> <sup>+</sup> *<sup>r</sup>*2.

interpolation problem is defined as follows:

*s*(**x***i*) =

*k* ∑ *j*=1

gives an algebraic system of equations for *<sup>λ</sup>* <sup>=</sup> {*λi*}*<sup>k</sup>*

*fj* = *f*(**x***j*). A simple form is to define

by *r*, and typical functions used in applications are, among others:

*s*(**x**) =

choice to study our problem, and we want to explore this alternative.

boundary operator. Suppose that we want to solve the problem

We consider a set of collocation points {**x***i*}*<sup>n</sup>*

*k* ∑ *j*=1

*c*<sup>2</sup> + *r*2.

.

$$\mathbf{u}\_h(\mathbf{x}) = \mathbf{u}^\mathrm{I}(\mathbf{x}) + \sum\_{j=1}^n \lambda\_j \, \mathbf{S}^{-1} \nabla \phi(\|\mathbf{x} - \mathbf{x}\_j\|). \tag{47}$$

**Example 8**. As a last example we include a 3–D numerical calculation using radial basis functions. The exact syntectic wind field for this example is **u** = (*x*, *y*, −2*z*). We dropped the vertical component so that **u**<sup>I</sup> = (*x*, *y*, 0). A multiquadric kernel with *c* = 11.33 was used, and we chose *α*<sup>1</sup> = *α*<sup>2</sup> = *α*<sup>3</sup> = 1. The collocation points were obtained by a 5 × 5 × 5 regular subdivision of Ω = (1, 2) × (0, 1) × (0, 1). Figure 7 shows the collocation points and the comparison between the exact and recovered wind field. The agreement is excellent, we obtained a relative error *er* <sup>=</sup> 1.98 <sup>×</sup> <sup>10</sup>−<sup>4</sup> and mean divergence *mdiv* <sup>=</sup> <sup>−</sup>5.59 <sup>×</sup> <sup>10</sup>−6.

Fig. 7. Collocation points (left), and comparison of exact and recovered wind fields (right).

#### **6. Conclusions**

We studied the problem of generating an adjusted wind field from horizontal wind data by different numerical techniques. We have shown that boundary conditions can significantly affect numerical solutions depending of how we treat artificial truncated boundaries. The usual methodology (*E1–algorithm*) does not produce satisfactory results close to the vertical boundaries due to the high gradients introduced by the term *<sup>S</sup>*−1∇*<sup>λ</sup>* in (2), when homogeneous Dirichlet boundary conditions are imposed there. The formulation of the problem as a saddle–point one, with a functional equation that has a self–adjoint and strongly elliptic operator, allows the use of an iterative conjugate gradient algorithm (the *CG–algorithm*). This new methodology, in the context of mass consistent models, produces

Finardi, S., Tinarelli, G., Nanni, A., Brusasca, G. & Carboni, G. (2010). Evaluation of a 3–d flow

Mass–Consistent Wind Field Models: Numerical Techniques by L2–Projection Methods 39

Glowinski, R. (2003). *Numerical Methods for Fluids (Part 3), Handbook of Numerical Analysis,*

Hardy, R. L. (1971). Multiquadric equations of topography and other irregular surfaces, *J.*

Kansa, E. J. (1990). Multiquadrics a scattered data approximation scheme with applications

Kitada, T., Igarashi, K. & Owada, M. (1986). Numerical analysis of air pollution in a

Kitada, T., Kaki, A., H., U. & K., P. L. (1983). Estimation of the vertical air motion from limited horizontal wind data–a numerical experiment, *Atmos. Environ* 17(11): 181–2192. Montero, G., Rodríguez, E., Montenegro, R., Escobar, J. M. & González-Yuste, J. M. (2005).

Narcowich, F. J. & Ward, J. D. (1994). Generalized hermite interpolation via matrix-valued

Núñez, M. A., Flores, C. & Juárez, H. (2007). Interpolation of hydrodynamic velocity data

Núñez, M. A., Flores, C. & Juárez, L. H. (2006). A study of hydrodynamic mass–consistent

Pennel, W. T. (1983). *An evaluation of the role of numerical wind field models in wind turbine*

Potter, B. & Butler, B. (2009). Using wind models to more effectively manage wildfire, *Fire*

Ratto, C. F. (1996). An overview of mass-consistent models, *in* D. P. Lalas & C. F.Ratto (eds),

Ratto, C. F., Festa, R., Romeo, C., Frumento, O. A. & Galluzzi, M. (1994). Mass–consistent

Ross, D. G., Smith, I. N., Manins, P. C. & Fox, D. G. (1988). Diagnostic wind field modeling for complex terrain: Model development and testing, *J. Appl. Meteor.* 27: 785–796. Ruhnau, P. & Schnorr, C. (2007). Optical stokes flow estimation: an imaging–based control

*siting. Technical Report PNL–SA–11129*, Battelle Memorial Institute, Pacific Northwest

*Modeling of Atmosphere Flow Fields*, World Scientific Publications, Place of publication,

models for wind fields over complex terrain: The state of the art, *Environ. Software*

to computational fluid dynamics ii: Solutions to parabolic, hyperbolic and elliptic

combined field of land/sea breeze and mountain/valley wind, *J. Clim. Appl. Meteorol.*

Genetic algorithms for an improved parameter estimation with local refinement of

GTOPO, . (1997). Gtopo30 documentation, section 7, *U. S. Geological Survey* pp. 211–224.

partial differential equations, *Comput. Math. Appl.* 19(8/9): 147–161.

tetrahedral meshes in a wind model, *Adv. Eng. Softw.* 36: 3–10.

conditionally positive definite functions, *Math. Comp.* 63: 661–687.

with the continuity equation, *J. Comput. Meth. Sci. Eng.* 7(1): 21–42.

models, *J. Comput. Meth. Sci. Eng.* 6: 1078–1089.

Laboratory, Richland, Washington.

approach, *Exp. Fluids* 42(1): 61–78.

*management today* 69(2): 40–46.

pp. 379–400.

9(4): 247–268.

*and Algorithms*, Springer–Verlag, Berlin.

*volume IX*, North–Holland, Amsterdam.

*Geophys. Res.* 76(8): 1905–1915.

25(6): 767–784.

URL: *www.scd.ucar.edu/ dss/datasests/ds758.0hmtl*

and pollutant dispersion modelling system to estimate climatological ground level concentrations in complex coastal sites, *Int. J. Environ. Pollut.* 16(1–6): 472–482. Flores, C. F., Juárez, L. H., Nuñez, M. A. & Sandoval, M. L. (2010). Algorithms for vector field generation in mass consistent models, *J. Numer. Methods for PDE* 26(4): 826–842. Girault, V. & Raviart, P. A. (1986). *Finite Element Methods for the Navier–Stokes Equations: Theory*

much better results, it does not involve the solution of differential equations, and it does not require boundary conditions for the multiplier. An optimal preconditioner for the conjugate gradient algorithm was introduced, and we hope, based on previous experience with the solution of the Stokes problem, to reduce the number of iterations of nearly a thousand to a few tens.

In an attempt to improve the numerical results obtained with the traditional approach, we introduced new boundary conditions for the multiplier on vertical boundaries. These boundary conditions are demonstrated to be physically and mathematically consistent. The numerical reconstruction of the wind field was improved by two orders of magnitude, but a comparable reduction on the weak divergence is not always obtained. However, the introduction of "ghost nodes" produces more satisfactory results, reducing both the relative error and the mean weak divergence by two or more orders of magnitude. On the other hand, we have just started to explore meshfree methods. In particular, radial basis collocation methods seem to be a very simple reliably alternative to the reconstruction of three–dimensional wind fields, according to the preliminary numerical results shown in this work.

The application of the different alternatives and methodologies presented here to the more realistic three–dimensional case is a continuation of the present work. Another interesting issue is the potential extension and application of these methodologies to other experimental fields, such as fluid dynamics and computer vision. In particular, the reconstruction of solenoidal velocity fields from experimental data, obtained through the *particle image velocimetry* technique, is an important issue, (Adrian, 2005). Its relation with computer vision is established by *optical flow* estimation, (Ruhnau & Schnorr, 2007).

#### **7. Acknowledgements**

Authors wish to express our deep appreciation to InTech for the kind invitation to contribute with this chapter. No doubt this was a great motivation. Special thanks to Ms. Tajana Jevtic and Jana Sertic, InTech process managers, for their guidance, patience and kindness throughout the editorial process.

#### **8. References**

Adrian, R. J. (2005). Twenty years of particle image velocimetry, *Exp. Fluids* 39(2): 159–169.


16 Will-be-set-by-IN-TECH

much better results, it does not involve the solution of differential equations, and it does not require boundary conditions for the multiplier. An optimal preconditioner for the conjugate gradient algorithm was introduced, and we hope, based on previous experience with the solution of the Stokes problem, to reduce the number of iterations of nearly a thousand to

In an attempt to improve the numerical results obtained with the traditional approach, we introduced new boundary conditions for the multiplier on vertical boundaries. These boundary conditions are demonstrated to be physically and mathematically consistent. The numerical reconstruction of the wind field was improved by two orders of magnitude, but a comparable reduction on the weak divergence is not always obtained. However, the introduction of "ghost nodes" produces more satisfactory results, reducing both the relative error and the mean weak divergence by two or more orders of magnitude. On the other hand, we have just started to explore meshfree methods. In particular, radial basis collocation methods seem to be a very simple reliably alternative to the reconstruction of three–dimensional wind fields, according to the preliminary numerical results shown in this

The application of the different alternatives and methodologies presented here to the more realistic three–dimensional case is a continuation of the present work. Another interesting issue is the potential extension and application of these methodologies to other experimental fields, such as fluid dynamics and computer vision. In particular, the reconstruction of solenoidal velocity fields from experimental data, obtained through the *particle image velocimetry* technique, is an important issue, (Adrian, 2005). Its relation with computer vision

Authors wish to express our deep appreciation to InTech for the kind invitation to contribute with this chapter. No doubt this was a great motivation. Special thanks to Ms. Tajana Jevtic and Jana Sertic, InTech process managers, for their guidance, patience and kindness

Adrian, R. J. (2005). Twenty years of particle image velocimetry, *Exp. Fluids* 39(2): 159–169. Bercovier, M. & Pironneau, O. (1979). Error estimates for the finite element method solution of the stokes problem in the primitive variables, *Numer. Math.* 33: 211–224. Buhmann, M. D. (2003). *Radial basis functions: theory and implementations*, Cambridge

Castino, F., Rusca, L. & Solari, G. (2003). Wind climate micro–zoning: a pilot application to liguria region (north–western italy), *J. Wind. Eng. Ind. Aerodyn.* 91(11): 1353–1375. Ciarlet, P. G. (2002). *The Finite Element Method for Elliptic Problems. Re–edited as Vol. 40, SIAM*,

Ferragut, L., Montenegro, R., Montero, G., Rodríguez, E., Asensio, M. L. & Escobar, J. M.

(2010). Comparison between 2.5-d and 3-d realistic models for wind field adjustmen,

is established by *optical flow* estimation, (Ruhnau & Schnorr, 2007).

a few tens.

work.

**7. Acknowledgements**

**8. References**

throughout the editorial process.

University Press, United Kingdom.

North–Holland Amsterdam (1970), Philadelphia, PA.

*J. Wind. Eng. Ind. Aerodyn.* 98(10–11): 548–558.


**3** 

Hwataik Han *Kookmin University* 

*Korea* 

**Ventilation Effectiveness Measurements** 

Ventilation effectiveness has been defined in various ways by many investigators. The term *ventilation efficiency* was first used by Yaglou and Witheridge (1937). They defined it as the ratio of the carbon dioxide concentration in a room to that in the extract duct. The ventilation was considered to be effective if the air contaminants in high concentration level are captured by the exhaust before it spreads out into the room. This definition has been the

The mathematical concepts of age and residence time were introduced in investigations of mixing characteristics in reactors by chemical engineers such as Danckwerts (1958) and Spalding (1958). They mentioned the similarity between the mixing of gases in reactors and the mixing of air in ventilated rooms. Sandberg (1981) first applied the concept of *age of air* to ventilation studies. He summarized various definitions of ventilation efficiency including relative efficiency, absolute efficiency, steady state efficiency, and transient efficiency. The sooner the supply air reaches a particular point in the room, the greater the air change efficiency at that point. This concept has been widely accepted by many researchers and

The ASHRAE Handbook (2009) states that *ventilation effectiveness* is a description of an air distribution system's ability to remove internally generated pollutants from a building, zone, or space, whereas *air change effectiveness* is defined as a description of a system's ability to deliver ventilation air to a building, zone, or space. Thus, ventilation effectiveness indicates the effectiveness of exhaust, whereas the air change effectiveness indicates the effectiveness of supply. However, the terminology *ventilation effectiveness* commonly

In this chapter, we provide a one-to-one analogy between exhaust effectiveness and supply effectiveness using the concept of the age of air. The meanings of local and overall values of supply and exhaust effectiveness need be understood appropriately in conjunction with the

We also extend the theory of the local mean age of air and local mean residual lifetime of contaminant to a space with multiple inlets and outlets. Theoretical considerations are given to derive the relations between the LMAs from individual inlets and the combined LMA of total supply air. In addition, the relations between the LMRs toward individual outlets and the combined LMR of total exhaust air are considered. These relations can be used to investigate the effect of each supply inlet and/or the contribution of each exhaust outlet in a

cornerstone of various definitions of ventilation efficiency ever since.

organizations throughout the world including ASHRAE and AIVC.

includes both supply and exhaust characteristics.

aforementioned definitions of ventilation effectiveness.

**1. Introduction** 

**Using Tracer Gas Technique** 


Hwataik Han *Kookmin University Korea* 

#### **1. Introduction**

18 Will-be-set-by-IN-TECH

40 Fluid Dynamics, Computational Modeling and Applications

Sasaki, Y. (1958). An objective analysis based on the variational method, *Journal Met. Soc. Japan*

Sherman, C. A. (1978). A mass–consistent model for wind fields over complex terrain, *J. Appl.*

Wang, Y., Williamson, C., Garvey, D., Chang, S. & Cogan, J. (2005). Application of a multigrid

Wendland, H. (2005). *Scattered data approximation*, Cambridge University Press, United

method to a mass–consistent diagnostic wind model, *J. Appl. Meteor.* 44(7): 1078–1089.

36: 77–88.

Kingdom.

*Meteor.* Vol. 17(3): 312–319.

Ventilation effectiveness has been defined in various ways by many investigators. The term *ventilation efficiency* was first used by Yaglou and Witheridge (1937). They defined it as the ratio of the carbon dioxide concentration in a room to that in the extract duct. The ventilation was considered to be effective if the air contaminants in high concentration level are captured by the exhaust before it spreads out into the room. This definition has been the cornerstone of various definitions of ventilation efficiency ever since.

The mathematical concepts of age and residence time were introduced in investigations of mixing characteristics in reactors by chemical engineers such as Danckwerts (1958) and Spalding (1958). They mentioned the similarity between the mixing of gases in reactors and the mixing of air in ventilated rooms. Sandberg (1981) first applied the concept of *age of air* to ventilation studies. He summarized various definitions of ventilation efficiency including relative efficiency, absolute efficiency, steady state efficiency, and transient efficiency. The sooner the supply air reaches a particular point in the room, the greater the air change efficiency at that point. This concept has been widely accepted by many researchers and organizations throughout the world including ASHRAE and AIVC.

The ASHRAE Handbook (2009) states that *ventilation effectiveness* is a description of an air distribution system's ability to remove internally generated pollutants from a building, zone, or space, whereas *air change effectiveness* is defined as a description of a system's ability to deliver ventilation air to a building, zone, or space. Thus, ventilation effectiveness indicates the effectiveness of exhaust, whereas the air change effectiveness indicates the effectiveness of supply. However, the terminology *ventilation effectiveness* commonly includes both supply and exhaust characteristics.

In this chapter, we provide a one-to-one analogy between exhaust effectiveness and supply effectiveness using the concept of the age of air. The meanings of local and overall values of supply and exhaust effectiveness need be understood appropriately in conjunction with the aforementioned definitions of ventilation effectiveness.

We also extend the theory of the local mean age of air and local mean residual lifetime of contaminant to a space with multiple inlets and outlets. Theoretical considerations are given to derive the relations between the LMAs from individual inlets and the combined LMA of total supply air. In addition, the relations between the LMRs toward individual outlets and the combined LMR of total exhaust air are considered. These relations can be used to investigate the effect of each supply inlet and/or the contribution of each exhaust outlet in a

Therefore, the *local supply index* can be understood as the ratio of the LMA at the exhaust to that at the point, and the *local exhaust index* as the ratio of the LMR at the supply to that at the point. The overall room effectiveness can be defined similarly. The definitions of supply and exhaust effectiveness are shown in Table 1. The subscript P is the location of interest,

It will be proved later in this chapter that the room averages of LMA and LMR are identical. Therefore, the overall room supply effectiveness and exhaust effectiveness should be the same. We do not need to distinguish the overall values, but we call this the *room ventilation effectiveness*. Note that supply effectiveness and exhaust effectiveness are meaningful only

SUPPLY EFFECTIVENESS EXHAUST EFFECTIVENESS

φ

< > φ

*p*

ε

Table 1. Definitions of supply and exhaust effectiveness using LMA and LMR




< >=

Tracer gas techniques have been widely used to measure air change rates and the air change effectiveness in a ventilated zone. Any measurable gas can be used as a tracer gas. It is desirable to follow air movements faithfully and for the gas to be nonreactive with other materials. Etheridge & Sandberg (1996) suggested that an ideal tracer gas should have the

A wide variety of gases have been employed as tracers. The characteristics of the most commonly used tracer gasses are given in Table 2. Carbon dioxide is a good tracer gas since it has a molecular weight similar to air and is mixed well with air. However, it has a background concentration of approximately 350 ppm, and it is produced by people and the

ε

τ

φ

Residual Life Time of Air

Local Exhaust Index *n* sup

φ

 φ= =

= LMR at supply/LMR at P

*P P*

*n* τ

φ

< >

*<sup>P</sup>* = Local Mean Residual Lifetime at P

= Room Average of LMR

Overall Room Exhaust Effectiveness

(1)

θ *ex* φsup *<sup>n</sup>* = =τ

and < > indicates the spatial average over the entire space.

for local values.

Age of Air

*<sup>P</sup>* = Local Mean Age at P

= LMA at exhaust/LMA at P

Overall Room Supply Effectiveness

Local Supply Index

*n* τ

θ

< >

**3. Tracer gas technique** 

following characteristics:


*n ex <sup>p</sup> P P*

τ θ

θ θ= =

= Room Average of LMA

θ

α

α

< >=

**3.1 Tracer gases** 


< > θ

space with multiple inlets and outlets. Three examples of tracer gas applications are included in this chapter.

#### **2. Definitions of ventilation effectiveness**

#### **2.1 Age and residual-life-time**

Consider a point P in a room with one supply air inlet and one return air exhaust. The *age of air* is the length of time required for the supply air to reach the point. As air can reach the point through various paths, the mean value of the ages at the point is called the *local mean age* (LMA) of the air at P. Likewise, the length of time required for the contaminant located at P to reach an exhaust is called the *residual lifetime of the contaminant* at P. The mean value through various paths is the *local mean residual lifetime* (LMR).

Local mean age represents the un-freshness of supply air so that it can be used as a local supply index at the point. Local mean residual lifetime represents the slowness of removal of the contaminant generated at the point, and can be used to represent a local exhaust index. The LMA and LMR represent the local supply and exhaust effectiveness, respectively, at the point in the room. We note that they depend on the room airflow pattern only, and should not be dependent on the source distribution of a contaminant in the space unless the contaminant concentration alters the airflow characteristics of the room.

Fig. 1. Concept of age and residual lifetime of indoor air.

#### **2.2 Supply and exhaust effectiveness**

A complete mixing condition is considered to be a reference condition we can use to define ventilation effectiveness. The *air change rate* is the number of room volumes of air supplied in one hour, and is defined as *Q/V* where *Q* is the volumetric flow rate of air into the room and *V* is the room volume. The room *nominal time constant* τ*<sup>n</sup>* is the inverse of the air change rate.

The local supply and exhaust indices are defined as the ratios of the local mean age and the local mean residual lifetime compared to the nominal time constant, respectively. These local indices can exceed 100% and can be as large as infinity.

Notice that the LMA at the exhaust means the total residence time of supply air in the space, which is the same as the LMR at the supply. We note that these values are equal to the nominal time constant.

$$
\theta\_{\rm ex} = \phi\_{\rm sup} = \pi\_n \tag{1}
$$

Therefore, the *local supply index* can be understood as the ratio of the LMA at the exhaust to that at the point, and the *local exhaust index* as the ratio of the LMR at the supply to that at the point. The overall room effectiveness can be defined similarly. The definitions of supply and exhaust effectiveness are shown in Table 1. The subscript P is the location of interest, and < > indicates the spatial average over the entire space.

It will be proved later in this chapter that the room averages of LMA and LMR are identical. Therefore, the overall room supply effectiveness and exhaust effectiveness should be the same. We do not need to distinguish the overall values, but we call this the *room ventilation effectiveness*. Note that supply effectiveness and exhaust effectiveness are meaningful only for local values.


Table 1. Definitions of supply and exhaust effectiveness using LMA and LMR

#### **3. Tracer gas technique**

#### **3.1 Tracer gases**

42 Fluid Dynamics, Computational Modeling and Applications

space with multiple inlets and outlets. Three examples of tracer gas applications are

Consider a point P in a room with one supply air inlet and one return air exhaust. The *age of air* is the length of time required for the supply air to reach the point. As air can reach the point through various paths, the mean value of the ages at the point is called the *local mean age* (LMA) of the air at P. Likewise, the length of time required for the contaminant located at P to reach an exhaust is called the *residual lifetime of the contaminant* at P. The mean value

Local mean age represents the un-freshness of supply air so that it can be used as a local supply index at the point. Local mean residual lifetime represents the slowness of removal of the contaminant generated at the point, and can be used to represent a local exhaust index. The LMA and LMR represent the local supply and exhaust effectiveness, respectively, at the point in the room. We note that they depend on the room airflow pattern only, and should not be dependent on the source distribution of a contaminant in the space unless the

**P**

Residence time

A complete mixing condition is considered to be a reference condition we can use to define ventilation effectiveness. The *air change rate* is the number of room volumes of air supplied in one hour, and is defined as *Q/V* where *Q* is the volumetric flow rate of air into the room and *V*

The local supply and exhaust indices are defined as the ratios of the local mean age and the local mean residual lifetime compared to the nominal time constant, respectively. These

Notice that the LMA at the exhaust means the total residence time of supply air in the space, which is the same as the LMR at the supply. We note that these values are equal to the

τ

*<sup>n</sup>* is the inverse of the air change rate.

Residual life time

φp=LMRp

included in this chapter.

**2.1 Age and residual-life-time** 

**2. Definitions of ventilation effectiveness** 

through various paths is the *local mean residual lifetime* (LMR).

contaminant concentration alters the airflow characteristics of the room.

Age

Fig. 1. Concept of age and residual lifetime of indoor air.

**2.2 Supply and exhaust effectiveness** 

nominal time constant.

is the room volume. The room *nominal time constant*

local indices can exceed 100% and can be as large as infinity.

θp=LMAp

Tracer gas techniques have been widely used to measure air change rates and the air change effectiveness in a ventilated zone. Any measurable gas can be used as a tracer gas. It is desirable to follow air movements faithfully and for the gas to be nonreactive with other materials. Etheridge & Sandberg (1996) suggested that an ideal tracer gas should have the following characteristics:


A wide variety of gases have been employed as tracers. The characteristics of the most commonly used tracer gasses are given in Table 2. Carbon dioxide is a good tracer gas since it has a molecular weight similar to air and is mixed well with air. However, it has a background concentration of approximately 350 ppm, and it is produced by people and the

The injection and distribution system releases an appropriate amount of tracer gas and distributes it into the zones. There are several means of releasing tracer gas, either manually or automatically. A graduated syringe or other containers of known volume may be used for simple manual injections. For automated injection systems, a compressed tracer gas supply is connected to a gas line with an electronic mass flow controller, or other tracer gas flow

An automatic distributing system includes a tubing network that dispenses a tracer gas via manifolds and automated valves, and pressure-operated valves that stop the flow from entering the tubing network when the tubing is not pressurized. There should be no leaks in the tubing. A mixing fan is frequently used for good mixing of tracer gases within a zone.

Air sampling can be achieved either manually or automatically. Manual samplers may include syringes, flexible bottles, or sampling bags with a capacity of at least three times the minimum sampler size of the gas analyzer used. Automatic samplers may utilize either a sampling network or automated samplers. Sampling networks consist of tubing, a manifold or selection switch that is typically solenoid-driven, and a pump that draws air samples through the network. Tracer gas molecules should not adhere to the tubing or manifold surfaces. Materials that absorb tracer gas may cause major inaccuracies in the measurement. There are various types of gas analyzers based on principles such as infrared spectroscopy, gas chromatography, or mass spectroscopy. A gas analyzer should be suited to the tracer

t

t

t

There are three commonly used methods of injecting a tracer gas: step-up, step-down, and pulse methods. The *step-up method* introduces a tracer gas at a given time and onward until

Table 3. Tracer injection methods and the corresponding concentration responses.

*<sup>M</sup>* C(t)

*<sup>M</sup>* C(t)

*<sup>M</sup>* C(t)

INJECTION MONITORING

time elapsed t

time elapsed t

t time elapsed

C(∞)

C(0)

**3.2.1 Injection and distribution** 

rate measurement and control devices.

**3.2.2 Sampling and monitoring** 

Step-up method

Step-down method

Pulse method

**3.3 Tracer injection methods** 

gas used, and the concentration range studied.

combustion of fuels in occupied spaces. The effect of the production should be compensated. Hydrogen gas and water vapor have also been also used as tracer gases by Dufton and Marley (1935). They pointed out problems related to phase change and adhesion on surfaces when using water vapor. Sulfur hexafluoride is also a common tracer gas used in various ventilated spaces. It is not present in normal ambient air and can be used at very low concentrations. This minimizes the amount of tracer gas needed for a test. However, it has a molecular weight approximately five times that of air, and it should be diluted and/or well mixed with the surrounding air during injection.


Table 2. Characteristics of commonly used tracer gases (Sandberg, 1981)

#### **3.2 Tracer gas system**

A general tracer gas system is composed of an injection and distribution system, a sampling and monitoring system, and a data acquisition and control system. An example of a typical experimental setup is shown in Fig. 2 (ASTM, 1993).

Fig. 2. Typical tracer gas experimental system.

#### **3.2.1 Injection and distribution**

44 Fluid Dynamics, Computational Modeling and Applications

combustion of fuels in occupied spaces. The effect of the production should be compensated. Hydrogen gas and water vapor have also been also used as tracer gases by Dufton and Marley (1935). They pointed out problems related to phase change and adhesion on surfaces when using water vapor. Sulfur hexafluoride is also a common tracer gas used in various ventilated spaces. It is not present in normal ambient air and can be used at very low concentrations. This minimizes the amount of tracer gas needed for a test. However, it has a molecular weight approximately five times that of air, and it should be diluted and/or

> Analytical method

GC-ECD

GC-ECD

A general tracer gas system is composed of an injection and distribution system, a sampling and monitoring system, and a data acquisition and control system. An example of a typical

> Data Acquisition and Control System

Purge Gas

Ventilated Space

dioxide 44 -56.6 1.98 IR 0.05-2000 Variable Slight

Detection range (ppm)

0.05-2000

0.05-2000 0.00002-0.5

> Distribution System

> > Injection System

0.001-0.05

Background

concentration Toxicity

well mixed with the surrounding air during injection.

Freon12 121 -29.8 5.13 IR

hexafluoride 146 -50.8 6.18 IR

experimental setup is shown in Fig. 2 (ASTM, 1993).

Sampling System

Monitoring System

Fig. 2. Typical tracer gas experimental system.

hexane 338 57.0 GC-ECD 10-8

Table 2. Characteristics of commonly used tracer gases (Sandberg, 1981)

Boiling point (C)

Density (15) (kg/*m*3)

Helium 4 -268.9 0.17 MS 5.24

oxide 44 -88.5 1.85 IR 0.05-2000 0.03

Gas Molecular

Carbon

Nitrous

Sulphur

Perfluoron

**3.2 Tracer gas system** 

weight

The injection and distribution system releases an appropriate amount of tracer gas and distributes it into the zones. There are several means of releasing tracer gas, either manually or automatically. A graduated syringe or other containers of known volume may be used for simple manual injections. For automated injection systems, a compressed tracer gas supply is connected to a gas line with an electronic mass flow controller, or other tracer gas flow rate measurement and control devices.

An automatic distributing system includes a tubing network that dispenses a tracer gas via manifolds and automated valves, and pressure-operated valves that stop the flow from entering the tubing network when the tubing is not pressurized. There should be no leaks in the tubing. A mixing fan is frequently used for good mixing of tracer gases within a zone.

#### **3.2.2 Sampling and monitoring**

Air sampling can be achieved either manually or automatically. Manual samplers may include syringes, flexible bottles, or sampling bags with a capacity of at least three times the minimum sampler size of the gas analyzer used. Automatic samplers may utilize either a sampling network or automated samplers. Sampling networks consist of tubing, a manifold or selection switch that is typically solenoid-driven, and a pump that draws air samples through the network. Tracer gas molecules should not adhere to the tubing or manifold surfaces. Materials that absorb tracer gas may cause major inaccuracies in the measurement. There are various types of gas analyzers based on principles such as infrared spectroscopy, gas chromatography, or mass spectroscopy. A gas analyzer should be suited to the tracer gas used, and the concentration range studied.

Table 3. Tracer injection methods and the corresponding concentration responses.

#### **3.3 Tracer injection methods**

There are three commonly used methods of injecting a tracer gas: step-up, step-down, and pulse methods. The *step-up method* introduces a tracer gas at a given time and onward until

concentration at the exhaust can be obtained from the total tracer generation rate in the space, which is the product of the generation rate per volume times the space volume. Thus,

*ex ex*

<sup>∞</sup> = = <sup>∞</sup>

*P p C C*

<sup>∞</sup> <sup>0</sup>

In order to measure the local mean residual lifetime at P, the injection point should be at P and the monitoring point should be at the exhaust. The LMR can be obtained using the

In a step-up method, the exhaust concentration reaches a steady state value ( ) *<sup>P</sup> Cex* ∞ as time goes to infinity. The mass balance should be satisfied; thus, the steady concentration at the exhaust should be equal to the total mass generation divided by the airflow rate, / *M Q* .

> ( ) <sup>1</sup> / <sup>1</sup> ( )

*P*

*C t dt M Q*

*<sup>M</sup> Q C t dt <sup>M</sup>*

( )

*C V M*

*P*

< ∞> <sup>=</sup>

*P*

< ∞> <sup>=</sup> <sup>∞</sup>

*C C* ( ) ( )

*P ex*

*n*

τ

where *M* is the contaminant generation rate at P. The first term in the integral is the total generation rate, and the second term is the rate of contaminant leaving the room through the extract duct. The integration of the difference up to the steady state results in the amount of contaminant left inside the room, which is called the *internal hold-up*. This is the product of the average room concentration times the room volume (Sandberg, 1981). Then, Eq. (4) can

*P ex*

<sup>=</sup> <sup>0</sup>

*p*

sup

*P*

*P*

sup 0 sup

*C t dt*

sup

*t C t dt*

*C t dt*

0

*p*

φ

= −

*ex <sup>p</sup>*

0

∞

= −⋅

sup

*P*

*P*

Table 4. Equations to calculate LMA and LMR for three injection methods

( )

*C t dt*

( ) (0)

( )

( )

0 sup ( ) <sup>1</sup>

<sup>∞</sup> = −

*C*

*P*

*P*

*C*

0

Therefore, the LMR using a step-up method can be written as

φ<sup>∞</sup>

∞

0

∞ <sup>⋅</sup> <sup>=</sup> 

*p*

*p*

θ<sup>∞</sup>

*p*

θ

equations in Table 4 similar to LMA equations.

θ

α

Step-up

Step-down

Pulse

**4.2 LMR measurements** 

be written as

θ

θ

( ) ( )

LOCAL MEAN AGE LOCAL MEAN RESIDUAL-LIFE-

φ

φ <sup>∞</sup> =

φ

(3)

TIME

*ex <sup>p</sup> <sup>P</sup>*

*ex <sup>p</sup> <sup>P</sup> ex C t dt C*

0

∞

0

*<sup>p</sup> <sup>P</sup> ex*

∞ <sup>⋅</sup> <sup>=</sup> 

 <sup>∞</sup> = − <sup>∞</sup>

( ) <sup>1</sup>

*ex*

*C*

( ) *P*

( ) (0) *P*

( )

( )

(4)

(5)

*P ex*

*t C t dt*

*C t dt*

*C t dt*

it reaches a steady state. The concentration response at a monitoring point is observed continuously. As a steady state is reached, the concentration is maintained at the steady state value. The *step-down method* is the opposite of the step-up method. Tracer injection is stopped abruptly and the concentration decay is monitored at a monitoring point. The concentration decays exponentially and approaches a background concentration. The concentration decay method is frequently used to measure air change rate starting from a uniform mixing of room air. Finally, the *pulse method* introduces a certain amount of tracer gas in a short period of time. A peak concentration response is detected at a monitoring point with a time delay. The concentration decays down to an initial concentration after the peak. Table 3 shows concentration responses according to the three injection methods.

#### **4. Measurements of ventilation effectiveness**

#### **4.1 LMA measurements**

In order to measure the local mean age at point P, the tracer injection point should be at a supply diffuser and the monitoring point is at point P as shown in Fig. 3. LMA can be obtained by integrating the area above the concentration curve (shaded area) divided by the steady state concentration after a step-up tracer injection. Similarly, it is the area under the concentration curve for a step-down method. In the case of a pulse method, it can be calculated using the first moment of the area under the concentration curve. The equations used to calculate LMAs are shown in Table 4 for three injection procedures. The equations are different from one injection method to another, but the result should be the same. The superscripts and the subscripts of the concentrations indicate the injection and the monitoring points, respectively.

Fig. 3. Injection and monitoring points for LMA and transient step-up response.

It is known that the LMA distribution in a space is equivalent to the steady concentration distribution with uniformly-distributed sources in the space (Han, 1992). The proof is given in the appendix. Thus,

$$\theta\_p = \frac{\overline{C}\_p(\approx)}{\dot{m}} \tag{2}$$

where *m* is the tracer generation rate per unit volume. In Eq. (2), *C* has an over-bar rather than a superscript, which represents a uniform tracer injection throughout the entire space. The local supply index, which is the ratio of the LMAs at the exhaust and at P, is calculated using the ratio of the steady concentrations with over-bars at those points. The steady concentration at the exhaust can be obtained from the total tracer generation rate in the space, which is the product of the generation rate per volume times the space volume. Thus,

$$\alpha\_p = \frac{\theta\_{\varepsilon\chi}}{\theta\_p} = \frac{\overline{\mathbb{C}}\_{\varepsilon\chi}(\Leftrightarrow)}{\overline{\mathbb{C}}\_p(\Leftrightarrow)}\tag{3}$$


Table 4. Equations to calculate LMA and LMR for three injection methods

#### **4.2 LMR measurements**

46 Fluid Dynamics, Computational Modeling and Applications

it reaches a steady state. The concentration response at a monitoring point is observed continuously. As a steady state is reached, the concentration is maintained at the steady state value. The *step-down method* is the opposite of the step-up method. Tracer injection is stopped abruptly and the concentration decay is monitored at a monitoring point. The concentration decays exponentially and approaches a background concentration. The concentration decay method is frequently used to measure air change rate starting from a uniform mixing of room air. Finally, the *pulse method* introduces a certain amount of tracer gas in a short period of time. A peak concentration response is detected at a monitoring point with a time delay. The concentration decays down to an initial concentration after the peak. Table 3 shows concentration responses according to the three injection methods.

In order to measure the local mean age at point P, the tracer injection point should be at a supply diffuser and the monitoring point is at point P as shown in Fig. 3. LMA can be obtained by integrating the area above the concentration curve (shaded area) divided by the steady state concentration after a step-up tracer injection. Similarly, it is the area under the concentration curve for a step-down method. In the case of a pulse method, it can be calculated using the first moment of the area under the concentration curve. The equations used to calculate LMAs are shown in Table 4 for three injection procedures. The equations are different from one injection method to another, but the result should be the same. The superscripts and the subscripts of the concentrations indicate the injection and the

EXHAUST

Fig. 3. Injection and monitoring points for LMA and transient step-up response.

*p C m*

θ

It is known that the LMA distribution in a space is equivalent to the steady concentration distribution with uniformly-distributed sources in the space (Han, 1992). The proof is given

where *m* is the tracer generation rate per unit volume. In Eq. (2), *C* has an over-bar rather than a superscript, which represents a uniform tracer injection throughout the entire space. The local supply index, which is the ratio of the LMAs at the exhaust and at P, is calculated using the ratio of the steady concentrations with over-bars at those points. The steady

*<sup>p</sup>*( )

Cp sup (t)

<sup>t</sup> LMAP

<sup>∞</sup> <sup>=</sup> (2)

Cp sup (∞)

**4. Measurements of ventilation effectiveness** 

**4.1 LMA measurements** 

monitoring points, respectively.

Monitoring

Age of Air

**P**

Injection

in the appendix. Thus,

SUPPLY

In order to measure the local mean residual lifetime at P, the injection point should be at P and the monitoring point should be at the exhaust. The LMR can be obtained using the equations in Table 4 similar to LMA equations.

In a step-up method, the exhaust concentration reaches a steady state value ( ) *<sup>P</sup> Cex* ∞ as time goes to infinity. The mass balance should be satisfied; thus, the steady concentration at the exhaust should be equal to the total mass generation divided by the airflow rate, / *M Q* . Therefore, the LMR using a step-up method can be written as

$$\begin{split} \phi\_p &= \int\_0^\infty \mathbf{1} - \frac{\mathbf{C}\_{ex}}{\dot{M} / Q} dt \\ &= \frac{1}{\dot{M}} \int\_0^\infty \dot{M} - Q \cdot \mathbf{C}\_{ex}{}^P(t) dt \end{split} \tag{4}$$

where *M* is the contaminant generation rate at P. The first term in the integral is the total generation rate, and the second term is the rate of contaminant leaving the room through the extract duct. The integration of the difference up to the steady state results in the amount of contaminant left inside the room, which is called the *internal hold-up*. This is the product of the average room concentration times the room volume (Sandberg, 1981). Then, Eq. (4) can be written as

$$\begin{split} \phi\_p &= \frac{<\mathcal{C}^P(\Leftrightarrow) > V}{\dot{M}} \\ &= \frac{<\mathcal{C}^P(\Leftrightarrow) > \mathfrak{r}\_n}{\mathcal{C}\_{\text{cc}}^{-P}(\Leftrightarrow)} \end{split} \tag{5}$$

used for LMA measurements, a monitoring point should be fixed at the exhaust, and a tracer should be injected at every point in the space repeatedly. The concentration response by simultaneous tracer injections can be obtained by superimposing every injection source present over the entire space, since the concentration equation is linear. Therefore, the room average exhaust effectiveness is the ratio of the exhaust concentration to the room average concentration with a uniformly distributed source superimposed in the space, which is identical to the steady method used to determine overall supply effectiveness. This concludes the proof that supply effectiveness equals the overall exhaust effectiveness of a given space,

> < >=< > ε

The room mean age or the room mean residual lifetime can also be obtained from the transient concentration responses at exhaust according to Table 5 for different tracer

> θ φ

θ

θ

Table 5. Equations used to calculate RMA and RMR for three injection methods

 φ <sup>∞</sup> < >=< >=

<sup>⋅</sup> < >=< >=

 φ

When there are multiple supply inlets, the LMA from one inlet is different than those from the other inlets. Consider a ventilated space configuration with two supply inlets as shown in Fig. 5. The airflow rates through the inlets are *Qa* and *Qb*, respectively, and room air is

Suppose we inject a tracer gas only at inlet *a* by a step-up method. The supply concentration is assumed to be 1.0 at inlet *a* and 0.0 at inlet *b*. The concentration response at P, *tC* )( *<sup>a</sup>*

shown in Fig. 5. *<sup>a</sup> LMAP* is the area above the curve (left-hatched area). Subscript P represents a monitoring point, and superscript *a* represents an injection location. The steady concentration ∞)(*<sup>a</sup> CP* has a value ranging between zero and unity because the supply

The response after a step-up injection at inlet *b* can be characterized similarly by *<sup>b</sup> LMAP* and ∞)(*<sup>b</sup> CP* . In this case, the non-dimensional supply concentration is 0.0 at inlet *a* and 1.0 at

 α

ROOM MEAN AGE = ROOM MEAN RESIDUAL-LIFE-TIME

0

<sup>∞</sup> < >=< >= ⋅ −

0

0

∞

0

∞

2 ( )

*ex*

*Q t C t dt V C t dt*

( ) <sup>1</sup> ( ) *<sup>Q</sup> C t ex <sup>t</sup> dt V C*

sup

2 sup

*ex*

sup

( ) (0) *<sup>Q</sup> C t ex t dt V C*

( )

*<sup>P</sup>* , is

<sup>∞</sup>

(9)

and that the room mean age of air is identical to the room mean residual lifetime:

Step-up method sup

injection methods (Kuehn *et al*., 1998).

Step-down method

Pulse method

exhausted through an outlet on the other side of the space.

concentration at inlet *a* is non-dimensionalized.

**5. Multiple inlets and outlets 5.1 LMA from multiple inlets** 

The local exhaust index can be obtained either from the definition of the LMR ratio, or by the ratio of the room average concentration to the exhaust concentration when a source is located at P. Thus,

$$\varepsilon\_p = \frac{\mathbb{C}\_{\text{ex}}{}^p(\Leftrightarrow)}{<\mathbb{C}^p(\Leftrightarrow)>} \tag{6}$$

Equation (6) looks quite similar to the classical definition by Yaglou and Witheridge (1937). They also defined ventilation efficiency as the ratio of the room average concentration to exhaust concentration for a given contaminant source. They understood this quantity as the overall efficiency of the room, not as a local efficiency at the given source location, though. The ratio shown in Eq. (6) is not the room exhaust index, but the local exhaust index for a given source located at P. Various definitions have been proposed for removal effectiveness by several authors (Sandberg and Sjoberg, 1983; Skaaret ,1986). Although there have been many studies on the measurement of LMA (Shaw *et al*., 1992; Han *et al*., 1999; Xing *et al*., 2001), the distributions of LMR have rarely been measured experimentally (Han *et al*., 2002).

Fig. 4. Injection and monitoring points for LMR and transient step-up response.

#### **4.3 Overall ventilation effectiveness**

The overall room effectiveness is the spatial average of local values over the entire space. As previously discussed, LMA can be obtained using transient and steady approaches. The steady method indicates that the spatial average of LMA is the spatial average of the steady concentration distribution with uniformly distributed tracer sources of unit strength, as follows:

$$<\theta> = \frac{<\overline{\mathbb{C}}(\simeq)>}{\dot{m}}\tag{7}$$

Therefore, the overall supply effectiveness; i.e., the ratio of LMA at exhaust to the room average LMA, equals the ratio of the concentration at exhaust to the spatial average of the steady concentration as follows:

$$<\alpha> = \frac{\overline{C}\_{\text{ex}}(\circ \circ)}{<\overline{C}(\circ \circ)>} \tag{8}$$

On the other hand, to obtain the overall exhaust effectiveness, LMR should be obtained at every internal point to calculate its spatial average over the entire space. Unlike the method

The local exhaust index can be obtained either from the definition of the LMR ratio, or by the ratio of the room average concentration to the exhaust concentration when a source is

> *ex <sup>p</sup> <sup>p</sup> C C*

 <sup>∞</sup> <sup>=</sup> < ∞>

Equation (6) looks quite similar to the classical definition by Yaglou and Witheridge (1937). They also defined ventilation efficiency as the ratio of the room average concentration to exhaust concentration for a given contaminant source. They understood this quantity as the overall efficiency of the room, not as a local efficiency at the given source location, though. The ratio shown in Eq. (6) is not the room exhaust index, but the local exhaust index for a given source located at P. Various definitions have been proposed for removal effectiveness by several authors (Sandberg and Sjoberg, 1983; Skaaret ,1986). Although there have been many studies on the measurement of LMA (Shaw *et al*., 1992; Han *et al*., 1999; Xing *et al*., 2001), the distributions of LMR have rarely been measured experimentally (Han *et al*., 2002).

> CP ex(t)

ε

Monitoring

Fig. 4. Injection and monitoring points for LMR and transient step-up response.

θ

α<sup>∞</sup> < >=

The overall room effectiveness is the spatial average of local values over the entire space. As previously discussed, LMA can be obtained using transient and steady approaches. The steady method indicates that the spatial average of LMA is the spatial average of the steady concentration distribution with uniformly distributed tracer sources of unit

> *C*( ) *m*

( ) ( ) *Cex C*

< ∞>

< >= (7)

(8)

< ∞>

Therefore, the overall supply effectiveness; i.e., the ratio of LMA at exhaust to the room average LMA, equals the ratio of the concentration at exhaust to the spatial average of the

On the other hand, to obtain the overall exhaust effectiveness, LMR should be obtained at every internal point to calculate its spatial average over the entire space. Unlike the method

EXHAUST Injection

**P**

**4.3 Overall ventilation effectiveness** 

strength, as follows:

steady concentration as follows:

Residual life time

( ) ( ) *p*

(6)

<sup>t</sup> LMRP

CP ex(∞)

located at P. Thus,

SUPPLY

used for LMA measurements, a monitoring point should be fixed at the exhaust, and a tracer should be injected at every point in the space repeatedly. The concentration response by simultaneous tracer injections can be obtained by superimposing every injection source present over the entire space, since the concentration equation is linear. Therefore, the room average exhaust effectiveness is the ratio of the exhaust concentration to the room average concentration with a uniformly distributed source superimposed in the space, which is identical to the steady method used to determine overall supply effectiveness. This concludes the proof that supply effectiveness equals the overall exhaust effectiveness of a given space, and that the room mean age of air is identical to the room mean residual lifetime:

$$<\xi \not\simeq \mathcal{Q} \tag{9}$$

The room mean age or the room mean residual lifetime can also be obtained from the transient concentration responses at exhaust according to Table 5 for different tracer injection methods (Kuehn *et al*., 1998).


Table 5. Equations used to calculate RMA and RMR for three injection methods

#### **5. Multiple inlets and outlets**

#### **5.1 LMA from multiple inlets**

When there are multiple supply inlets, the LMA from one inlet is different than those from the other inlets. Consider a ventilated space configuration with two supply inlets as shown in Fig. 5. The airflow rates through the inlets are *Qa* and *Qb*, respectively, and room air is exhausted through an outlet on the other side of the space.

Suppose we inject a tracer gas only at inlet *a* by a step-up method. The supply concentration is assumed to be 1.0 at inlet *a* and 0.0 at inlet *b*. The concentration response at P, *tC* )( *<sup>a</sup> <sup>P</sup>* , is shown in Fig. 5. *<sup>a</sup> LMAP* is the area above the curve (left-hatched area). Subscript P represents a monitoring point, and superscript *a* represents an injection location. The steady concentration ∞)(*<sup>a</sup> CP* has a value ranging between zero and unity because the supply concentration at inlet *a* is non-dimensionalized.

The response after a step-up injection at inlet *b* can be characterized similarly by *<sup>b</sup> LMAP* and ∞)(*<sup>b</sup> CP* . In this case, the non-dimensional supply concentration is 0.0 at inlet *a* and 1.0 at

concentration is the average of individual exhaust concentrations weighted by the airflow

1

**t**

LMAp

LMRp

b

= + (13)

LMRp

a+b

a

(14)

(15)

**Cb**

**Ca+bP (t)**

**P (∞)**

**Ca**

**P (∞)**

**Cex(t)**

**Qa**

**Qb**

Fig. 6. Local mean residual lifetime at P and concentration responses at the exhausts.

() () () *a b P P ex a b Q Q Ct C t C t Q Q*

> ( ) <sup>1</sup> ( )

*a P*

( ) <sup>1</sup> ( )

*b*

*C dt*

( ) ( ) () ()

( )

*a b ab P P PP a b ab*

*Q Q QQ C C Ct Ct Q Q QQ dt*

*ex*

*C*

Therefore, the combined LMR is the weighted average of the individual LMRs. The weighting factors are the percentages of the contaminant removal rates through the corresponding exhaust outlets. They can be understood as the contribution factors of the

∞+ ∞ − + <sup>=</sup> <sup>∞</sup>

*P P*

The individual LMRs to the outlets can be obtained by integrating the areas above the

0

*<sup>P</sup> <sup>a</sup> <sup>a</sup> <sup>P</sup> <sup>P</sup>*

∞

*C t LMR dt C*

= − <sup>∞</sup>

= − <sup>∞</sup>

*b b P P*

Similarly, the combined LMR can be obtained by the area above the average exhaust concentration curve. The combined LMR can be rearranged using Eq. (13), and can be

*C t dt LMR*

∞

0

rates through the outlets, as follows:

**P**

corresponding concentration curves:

expressed with the individual LMRs as follows:

0

∞

*ex <sup>P</sup>*

*C t LMR dt C*

= − <sup>∞</sup>

( ) <sup>1</sup> ( )

*ex*

( ) ( ) ( ) ( )

∞ ∞ = ⋅+ ⋅ ∞ ∞

*aa ba a b*

*CQ CQ LMR LMR CQ CQ*

*P P*

*ex ex a b a b P P*

*M M LMR LMR M M*

=⋅ +⋅

 

0

individual outlets for a given tracer source at P.

∞

LMRp

a

LMRp

b

**Q**

inlet *b*. We note that the steady state concentration ∞)(*<sup>b</sup> CP* is complementary to ∞)(*<sup>a</sup> CP* , since the inlet concentration boundary conditions are switched.

In the case of simultaneous tracer injections at both inlets, the concentration response at point P is given by the addition of the concentration responses from individual injections, as shown in Fig. 5.

Fig. 5. Local mean age from individual supply inlets and concentration responses at P.

$$\mathbf{C}\_{P}^{\
u+b}(t) = \mathbf{C}\_{P}^{\
u}(t) + \mathbf{C}\_{P}^{\
u}(t) \tag{10}$$

This is because the indoor airflow pattern remains unchanged and the governing equation is linear with respect to concentration. Concentrations reach 1.0 at all internal points as a steady state is reached. Thus,

$$1 = C\_p \, ^a(\ast \ast) + C\_p \, ^b(\ast \ast) \tag{11}$$

The combined LMA is the area above the combined concentration curve, which is the shaded area in Fig. 5. The relations between the LMAs can be derived as follows (Han *et al*., 2010):

$$LMA\_P = \mathbb{C}\_P{}^a \text{(\rightsquigarrow)} \cdot LMA\_P{}^a + \mathbb{C}\_P{}^b \text{(\rightsquigarrow)} \cdot LMA\_P{}^b \tag{12}$$

Therefore, the combined LMA is the weighted average of the LMAs from each individual inlet, and the weighting factors for calculating the average are the corresponding steady state concentrations at the point. The steady state concentrations can be considered to be the contribution factors of the corresponding inlets for characterizing the supply air conditions at the point.

#### **5.2 LMR to multiple outlets**

If there are multiple outlets, the contribution of each outlet is different with respect to eliminating contaminants generated in a space, depending on the relative source locations. Consider a case with two outlets with exhaust flow rates of *Qa* and *Qb* as shown in Fig. 6. The time for the contaminant generated at P to reach one exhaust, *<sup>a</sup> LMRP* , is different from the time to reach the other, *<sup>b</sup> LMRP* . The total amount of contaminants exhausted by one outlet is different from that exhausted by the other. Figure 6 shows concentration responses at the exhausts according to a step-up injection at point P. The combined exhaust

inlet *b*. We note that the steady state concentration ∞)(*<sup>b</sup> CP* is complementary to ∞)(*<sup>a</sup> CP* ,

In the case of simultaneous tracer injections at both inlets, the concentration response at point P is given by the addition of the concentration responses from individual injections, as

1

**t**

**Cp**

**Cp**

**a (∞)**

**Cp**

**b(∞)**

**a+b (t)**

LMAp

LMAp

( ) ( ) *a ab b LMA C LMA C LMA PP P P P* = ∞⋅ + ∞⋅ (12)

b

<sup>+</sup> = + (10)

LMAp

a+b

a

**CP(t)**

**Q**

Fig. 5. Local mean age from individual supply inlets and concentration responses at P.

This is because the indoor airflow pattern remains unchanged and the governing equation is linear with respect to concentration. Concentrations reach 1.0 at all internal points as a

1 () () *a b* = ∞+ ∞ *C C P P* (11)

The combined LMA is the area above the combined concentration curve, which is the shaded area in Fig. 5. The relations between the LMAs can be derived as follows (Han *et* 

Therefore, the combined LMA is the weighted average of the LMAs from each individual inlet, and the weighting factors for calculating the average are the corresponding steady state concentrations at the point. The steady state concentrations can be considered to be the contribution factors of the corresponding inlets for characterizing the supply air conditions

If there are multiple outlets, the contribution of each outlet is different with respect to eliminating contaminants generated in a space, depending on the relative source locations. Consider a case with two outlets with exhaust flow rates of *Qa* and *Qb* as shown in Fig. 6. The time for the contaminant generated at P to reach one exhaust, *<sup>a</sup> LMRP* , is different from the time to reach the other, *<sup>b</sup> LMRP* . The total amount of contaminants exhausted by one outlet is different from that exhausted by the other. Figure 6 shows concentration responses at the exhausts according to a step-up injection at point P. The combined exhaust

since the inlet concentration boundary conditions are switched.

**P**

() () () *ab a b C tCtCt P PP*

LMAp

steady state is reached. Thus,

**5.2 LMR to multiple outlets** 

a

LMAp

b

shown in Fig. 5.

**Qa**

**Qb**

*al*., 2010):

at the point.

concentration is the average of individual exhaust concentrations weighted by the airflow rates through the outlets, as follows:

Fig. 6. Local mean residual lifetime at P and concentration responses at the exhausts.

$$\mathbf{C}\_{ex}(t) = \frac{Q\_a}{Q} \mathbf{C}\_a^{\;P}(t) + \frac{Q\_b}{Q} \mathbf{C}\_b^{\;P}(t) \tag{13}$$

The individual LMRs to the outlets can be obtained by integrating the areas above the corresponding concentration curves:

$$\begin{aligned} LMR\_{P}{}^{a} &= \int\_{0}^{\infty} \mathbf{1} - \frac{\mathbf{C}\_{a}{}^{P} \{t\}}{\mathbf{C}\_{a}{}^{P} \{\r\mathbf{s}\}} dt \\ LMR\_{P}{}^{b} &= \int\_{0}^{\infty} \mathbf{1} - \frac{\mathbf{C}\_{b}{}^{P} \{t\} dt}{\mathbf{C}\_{b}{}^{P} \{\r\mathbf{s}\} dt} \end{aligned} \tag{14}$$

Similarly, the combined LMR can be obtained by the area above the average exhaust concentration curve. The combined LMR can be rearranged using Eq. (13), and can be expressed with the individual LMRs as follows:

$$\begin{split} \text{LMR}\_{P} &= \int\_{0}^{\infty} \mathbf{1} - \frac{\mathbf{C}\_{ex}(t)}{\mathbf{C}\_{ex}(\circ \circ \circ)} dt \\ &= \int\_{0}^{\infty} \frac{\left(\frac{\mathbf{Q}\_{a}}{Q} \mathbf{C}\_{a}{}^{P}(\circ \circ) + \frac{\mathbf{Q}\_{b}}{Q} \mathbf{C}\_{b}{}^{P}(\circ \circ)\right) - \left(\frac{\mathbf{Q}\_{a}}{Q} \mathbf{C}\_{a}{}^{P}(t) + \frac{\mathbf{Q}\_{b}}{Q} \mathbf{C}\_{b}{}^{P}(t)\right)}{\mathbf{C}\_{ex}(\circ \circ)} dt \\ &= \frac{\mathbf{C}\_{a}{}^{P}(\circ \circ) \mathbf{Q}\_{a}}{\mathbf{C}\_{ex}(\circ \circ) \mathbf{Q}} \cdot LMR\_{P}{}^{a} + \frac{\mathbf{C}\_{b}{}^{P}(\circ \circ) \mathbf{Q}\_{a}}{\mathbf{C}\_{ex}(\circ \circ) \mathbf{Q}} \cdot LMR\_{P}{}^{b} \\ &= \frac{\dot{M}\_{a}}{M} \cdot LMR\_{P}{}^{a} + \frac{\dot{M}\_{b}}{M} \cdot LMR\_{P}{}^{b} \end{split} \tag{15}$$

Therefore, the combined LMR is the weighted average of the individual LMRs. The weighting factors are the percentages of the contaminant removal rates through the corresponding exhaust outlets. They can be understood as the contribution factors of the individual outlets for a given tracer source at P.

The air change per hour (ACH) is the ratio of volumetric flow rate to the volume of the room. By a simple mathematical manipulation, the relation of ACH between the model and

2

The airflow and temperature conditions of the chamber were adjusted and checked until the steady state was reached. The sampling tube was positioned at a monitoring point using the three-dimensional (3-D) traversing system. Using a syringe, 3 mL of SF6 gas was injected into the supply duct. The gas monitor started to take data at the same time as the gas injection, which works on the principle of electron capture gas chromatography. Concentration data were recorded every 70 s until the concentration fell within 1% of the maximum concentration. The same measurement was repeated with a delayed injection by 35 s to double the number of data points. The sampling port was then moved to the next position, and the aforementioned procedure was repeated to cover the entire cross-section at

In order to investigate the effect of thermal buoyancy, three different temperature conditions were tested: isothermal, cooling, and heating. The experimental conditions and measurements are summarized in Table 6. The values of the corresponding full scale

> 17.5 1.032 345 62.6 (15.6) 24.4 (24.4) 0.1 24.4 (24.4) 24.4 (24.4) 0.1 24.4 (24.4) 24.4 0 (0) 3374 0

<sup>1</sup> () () *ACH ACH <sup>m</sup> <sup>f</sup> <sup>N</sup>*= (20)

Isothermal Cooling Heating

11.2 0.813 272 49.4 (12.4) -1.6 (19.2) 1.2 27.6 (22.9) 46.0 (25.2) 4.4 24.3 (22.5) 22.2 47.6 (6.0) 3128 0.1215

22.5 1.209 404 73.3 (18.3) 57.1 (30.5) 1.9 41.4 (28.5) -3.7 (22.9) 2.2 23.4 (26.3) 26.7 -60.8 (-7.6) 3293 -0.0690

the prototype becomes

**6.1.4 Experimental procedure** 

the center of the chamber.

situation are shown in parentheses.

Pressure drop across nozzle [mmH2O] Supply velocity at diffuser [m/s] Supply airflow rate [m3/h] Air change per hour [ACH] Supply air temperature [°C]

STD of supply air temperature [°C] Exhaust air temperature [°C] Wall temperature [°C] STD of wall temperature [°C] Ceiling temperature [°C] Mean temperature [°C] Twall – Tsupply [°C] Reynolds number Archimedes number

\*Numbers in ( ) indicate the corresponding values in full-scale situations.

Table 6. Experimental conditions and measurements (Han, 1999)

#### **6. Examples of tracer gas applications**

#### **6.1 Effect of supply air temperature on LMA distributions 6.1.1 Problem description**

It is often observed that fresh air supplied to a space is bypassed directly to an exhaust without contributing to effective room ventilation. Bypass affects the ventilation effectiveness of the room significantly. It is quite common in office buildings, especially when warm air is supplied from ceiling diffusers in the winter season. The following example considers the effect of supply air temperature on LMA distribution in a rectangular space with a diffuser and a return grill on the ceiling.

#### **6.1.2 Experimental setup**

A schematic of the experimental chamber is shown in Fig. 7. The chamber measures 1.95 m × 1.95 m × 1.45 m. The height of 1.45 m is about one-half of a full-scale office room. The interior surfaces (walls and floors) are made of aluminum panels. By circulating temperature-controlled fluid through the passages embedded in each panel, the temperatures of the walls and floors are precisely controlled. The ceiling is insulated with polystyrene insulation boards of 50 mm thickness. Air is supplied to the chamber through three linear sections that measures 0.635 m in length and 0.0508 m in width each. The three sections are aligned to form a 0.0508 m × 1.905 m straight slot inlet in the ceiling. The return slot is identical to the inlet and is also placed in the ceiling. This configuration produces a two-dimensional (2-D) flow in the chamber. A detailed description of the physical structure and the control system of the chamber is given by Corpron (1992).

#### **6.1.3 Similitude**

In this study, The Reynolds number and Archimedes number are considered important in simulating the full-scale conditions. These dimensionless numbers are defined as follows:

$$\text{Re} = \frac{\rho uL}{\mu} \quad \text{or} \quad \frac{inertia\,force}{viscous\,force} \tag{16}$$

$$\text{Ar} = \frac{\beta \text{g} \, \text{L} \, \text{LT}}{\text{u}^2} \quad \text{or} \quad \frac{\text{buoyancy force}}{\text{inertia force}} \tag{17}$$

For a half-scale model, the characteristic length is related as *L NL <sup>m</sup>* = *<sup>f</sup>* , where *N* equals 0.5. Subscript *m* stands for the model and *f* represents the full scale. As the thermodynamic properties ρ, μ, and β are assumed to be constant for both, the characteristic velocity of the model needs to be increased by a factor of 1/*N*. Also, the temperature difference needs to be increased by a factor of <sup>3</sup> 1 /*N* for similarity. Thus,

$$
\mu\_m = \frac{1}{N} \mu\_f \tag{18}
$$

$$\left(\Delta T\right)\_m = \frac{1}{N^3} (\Delta T)\_f \tag{19}$$

The air change per hour (ACH) is the ratio of volumetric flow rate to the volume of the room. By a simple mathematical manipulation, the relation of ACH between the model and the prototype becomes

$$\text{(ACH)}\_{m} = \frac{1}{N^2} \text{(ACH)}\_{f} \tag{20}$$

#### **6.1.4 Experimental procedure**

52 Fluid Dynamics, Computational Modeling and Applications

It is often observed that fresh air supplied to a space is bypassed directly to an exhaust without contributing to effective room ventilation. Bypass affects the ventilation effectiveness of the room significantly. It is quite common in office buildings, especially when warm air is supplied from ceiling diffusers in the winter season. The following example considers the effect of supply air temperature on LMA distribution in a rectangular

A schematic of the experimental chamber is shown in Fig. 7. The chamber measures 1.95 m × 1.95 m × 1.45 m. The height of 1.45 m is about one-half of a full-scale office room. The interior surfaces (walls and floors) are made of aluminum panels. By circulating temperature-controlled fluid through the passages embedded in each panel, the temperatures of the walls and floors are precisely controlled. The ceiling is insulated with polystyrene insulation boards of 50 mm thickness. Air is supplied to the chamber through three linear sections that measures 0.635 m in length and 0.0508 m in width each. The three sections are aligned to form a 0.0508 m × 1.905 m straight slot inlet in the ceiling. The return slot is identical to the inlet and is also placed in the ceiling. This configuration produces a two-dimensional (2-D) flow in the chamber. A detailed description of the physical structure

In this study, The Reynolds number and Archimedes number are considered important in simulating the full-scale conditions. These dimensionless numbers are defined as follows:

Re *uL inertia force*

<sup>2</sup> Ar *gL T buoyancy force u inertia force*

For a half-scale model, the characteristic length is related as *L NL <sup>m</sup>* = *<sup>f</sup>* , where *N* equals 0.5. Subscript *m* stands for the model and *f* represents the full scale. As the thermodynamic

model needs to be increased by a factor of 1/*N*. Also, the temperature difference needs to be

1

3 <sup>1</sup> () () *T T <sup>m</sup> <sup>f</sup> <sup>N</sup>*

ρ

β

increased by a factor of <sup>3</sup> 1 /*N* for similarity. Thus,

μ

*viscous force*

are assumed to be constant for both, the characteristic velocity of the

= ∝ (16)

<sup>Δ</sup> = ∝ (17)

*u u <sup>m</sup> <sup>f</sup> <sup>N</sup>*= (18)

Δ= Δ (19)

**6. Examples of tracer gas applications** 

**6.1.1 Problem description** 

**6.1.2 Experimental setup** 

**6.1.3 Similitude** 

properties

ρ, μ, and β

**6.1 Effect of supply air temperature on LMA distributions** 

space with a diffuser and a return grill on the ceiling.

and the control system of the chamber is given by Corpron (1992).

The airflow and temperature conditions of the chamber were adjusted and checked until the steady state was reached. The sampling tube was positioned at a monitoring point using the three-dimensional (3-D) traversing system. Using a syringe, 3 mL of SF6 gas was injected into the supply duct. The gas monitor started to take data at the same time as the gas injection, which works on the principle of electron capture gas chromatography. Concentration data were recorded every 70 s until the concentration fell within 1% of the maximum concentration. The same measurement was repeated with a delayed injection by 35 s to double the number of data points. The sampling port was then moved to the next position, and the aforementioned procedure was repeated to cover the entire cross-section at the center of the chamber.

In order to investigate the effect of thermal buoyancy, three different temperature conditions were tested: isothermal, cooling, and heating. The experimental conditions and measurements are summarized in Table 6. The values of the corresponding full scale situation are shown in parentheses.


\*Numbers in ( ) indicate the corresponding values in full-scale situations.

Table 6. Experimental conditions and measurements (Han, 1999)

125.9 138.0 130.9 139.3 126.5 125.9

119.4

145.0

122.4

127.0

118.0

10

119.4

124.2 124.2 126.0 128.7 124.5 125.1

**6.1.6 Concluding remarks** 

121.7

127.3

125.9

120.7

temperature.

117.3

124.1

20

30

40

50

120.6

10

20

30

40

50

116.0

124.6

113.0

125.1 132.8 148.2 145.9 143.7 126.4

10 20 30 40 50 60 70

10 20 30 40 50 60 70

the velocity patterns obtained by Liang (1994).

enhanced compared to the isothermal condition.

120.7 127.7 126.1 126.8 124.7 126.4 108.8 116.9 94.4 122.0 86.4 68.1

Fig. 8. LMA distribution and velocity vector fields for isothermal condition.

124.3 127.7 132.8 126.6 125.5 102.8

93.9

Fig. 9. Local mean age distributions for cooling and heating conditions (Han, 1999).

86.0

311.1

10

(a) Cooling condition. (b) Heating condition.

Using a pulsed injection method using SF6 tracer gas, LMA distributions were measured in a half-scale thermal chamber. Boundary conditions were applied that simulated isothermal, heating, and cooling conditions by controlling the supply air temperature and the wall

1. The LMA distribution was found to be closely related to the velocity distribution in the chamber. The results for LMA distributions are in good agreement qualitatively with

2. For an isothermal condition, the largest LMA occurred at the center and not at the corners, which indicates that there was a large recirculating zone at the center. During a cooling operation, supply air penetrated deeply into the chamber, and mixing was

3. For a heating condition, there was a large variation of local mean ages due to thermal stratification in the chamber. It can be concluded that local ventilation effectiveness in

the lower part of a room can be very poor under heating operations.

20

30

40

50

505.3 468.3 270.9 158.7 149.4 145.8

434.3

161.2

138.4

332.0

215.7

134.5

389.4 542.1 308.4 137.9 66.2 144.2

10 20 30 40 50 60 70

487.0 387.6 230.9 111.3 58.5 56.1

431.9

292.0

105.6

82.0

428.5

125.1

66.8

126.2 60.0

124.0

112.1

123.1

129.2

129.2

126.2 84.9

82.8

104.0

(a) LMA distribution (Han, 1999). (b) Velocity vectors (Liang, 1994).

135.7

85.2

76.1

130.7

147.0

118.9

115.2 95.3

Fig. 7. Cross-section of a thermal chamber (Han, 1999).

#### **6.1.5 Results and discussion**

Figure 8 shows LMA contours superimposed with the LMA data measured at 36 locations for an isothermal condition. The LMA distribution is nearly uniform over the entire crosssection except at corners. The maximum LMA was observed at the center, which indicates there is a large recirculation in the middle of the chamber. A velocity vector drawing by Liang (1994) is also shown in Fig. 8. The air jet from the supply inlet moves downward and leaves the chamber through the exhaust after making a large clockwise circulation in the space. The supply jet is attached to the right wall, and then separates before it hits the floor. This tendency for flows to attach to walls is known as the Coanda effect. It is interesting to note that the distribution of local mean age in the space shows a good overall picture of the airflow pattern in the space.

For a cooling condition, the supply air temperature is lower than the room temperature and the buoyant force acts downward, which is the same as the direction of the supply air. Figure 9(a) shows the LMA distribution in the chamber. Assisted by buoyancy, the mixing of the flow is enhanced and the local mean age values are more uniformly distributed compared to the isothermal condition. The location of maximum LMA is shifted downward and to the right in comparison with the isothermal condition. The maximum LMA value is less than that in the isothermal case.

In a heating condition, the thermal buoyancy opposes the inertial effect. The local mean age distribution is shown in Fig. 9(b). A large variation in LMA can be observed in the chamber because of thermal stratification. The variation is small at the upper part and large at the lower part of the chamber. The air jet from the supply port does not seem to penetrate into the space effectively; rather, it short circuits to the exhaust duct. Liang (1994) observed that the flow field under the heating condition was unstable and the supply jet oscillated slowly within the chamber. Because of the oscillatory behavior, velocity vectors could not be measured in the experiment, and only the frequency of the oscillatory motion was reported. The room mean ages obtained by integrating the local mean age values over the entire space give 118 s, 120 s, and 234 s for the isothermal, cooling, and heating conditions, respectively.

Figure 8 shows LMA contours superimposed with the LMA data measured at 36 locations for an isothermal condition. The LMA distribution is nearly uniform over the entire crosssection except at corners. The maximum LMA was observed at the center, which indicates there is a large recirculation in the middle of the chamber. A velocity vector drawing by Liang (1994) is also shown in Fig. 8. The air jet from the supply inlet moves downward and leaves the chamber through the exhaust after making a large clockwise circulation in the space. The supply jet is attached to the right wall, and then separates before it hits the floor. This tendency for flows to attach to walls is known as the Coanda effect. It is interesting to note that the distribution of local mean age in the space shows a good overall picture of the

For a cooling condition, the supply air temperature is lower than the room temperature and the buoyant force acts downward, which is the same as the direction of the supply air. Figure 9(a) shows the LMA distribution in the chamber. Assisted by buoyancy, the mixing of the flow is enhanced and the local mean age values are more uniformly distributed compared to the isothermal condition. The location of maximum LMA is shifted downward and to the right in comparison with the isothermal condition. The maximum LMA value is

In a heating condition, the thermal buoyancy opposes the inertial effect. The local mean age distribution is shown in Fig. 9(b). A large variation in LMA can be observed in the chamber because of thermal stratification. The variation is small at the upper part and large at the lower part of the chamber. The air jet from the supply port does not seem to penetrate into the space effectively; rather, it short circuits to the exhaust duct. Liang (1994) observed that the flow field under the heating condition was unstable and the supply jet oscillated slowly within the chamber. Because of the oscillatory behavior, velocity vectors could not be measured in the experiment, and only the frequency of the oscillatory motion was reported. The room mean ages obtained by integrating the local mean age values over the entire space give 118 s, 120 s, and 234 s for the isothermal, cooling, and heating conditions, respectively.

Fig. 7. Cross-section of a thermal chamber (Han, 1999).

**6.1.5 Results and discussion** 

airflow pattern in the space.

less than that in the isothermal case.

Fig. 8. LMA distribution and velocity vector fields for isothermal condition.

Fig. 9. Local mean age distributions for cooling and heating conditions (Han, 1999).

#### **6.1.6 Concluding remarks**

Using a pulsed injection method using SF6 tracer gas, LMA distributions were measured in a half-scale thermal chamber. Boundary conditions were applied that simulated isothermal, heating, and cooling conditions by controlling the supply air temperature and the wall temperature.


For the LMR measurements, tracer gas was injected through a porous sphere at a point in the chamber and the transient tracer gas concentration variation was measured at the exhaust. After the tracer gas was exhausted completely from the chamber, the injection sphere was moved to another position. The procedure was repeated for other internal points. There are 15 injection points equally spaced in the center plane of the chamber. For LMA measurements, all the experimental conditions were identical to the LMR case, but tracer gas was injected at a supply duct. Then, transient tracer gas concentration variation was measured at the internal points. Experiments were conducted for three different

Flow visualization results are shown in Fig. 11 for three different exhaust locations. The air change rates are 12ACHs. For Case 1, the air supplied in the horizontal direction moved toward the lower-left exhaust in the diagonal direction. The room air formed two large recirculating flows at the upper-left and lower-right corners. For Case 2, the air supplied along the ceiling changed its direction by the opposite wall and made a large counterclockwise circulation in the chamber. We note that the airflow pattern was quite similar to a complete mixing condition. For Case 3, supplied air faced directly toward the exhaust. The room air was mostly stagnant, but with a slow recirculation due to the viscous action of the

(a) Case 1 (b) Case 2 (c) Case 3

Contours of LMA and LMR are plotted in Fig. 12 for Case 1. It can be seen that both of the distributions are closely related to the airflow pattern shown in Fig. 11(a). LMA and LMR values are large within recirculating zones. LMA is small adjacent to the supply inlet and large adjacent to the exhaust, whereas LMR is small adjacent to the exhaust and large adjacent to the supply inlet. Figure 13 shows LMA and LMR distributions for Case 2. The LMA near ceiling is relatively small, whereas the LMR near floor is small. Both have large values within a large recirculation zone at the center. Figure 14 shows the results for Case 3. The LMA and LMR are small near ceiling, and large adjacent to the floor. The airflow pattern in Fig. 11(c) indicates there was large stagnant recirculation in the lower part of the space. The tracer gas diffused out into the lower part could not be exhausted

Fig. 11. Flow visualization results for three inlet-outlet configurations.

exhaust locations under isothermal room temperature conditions.

**6.2.3 Experimental procedure** 

**6.2.4 Results and discussion** 

bypassing flow along the ceiling.

effectively.

Further research needs to be done to improve the tracer gas technique and to apply the technique to various applications.

#### **6.2 Effect of inlet/outlet configurations on LMA and LMR distributions 6.2.1 Problem description**

The airflow pattern in a ventilated space varies according to the locations of supply inlets and exhaust outlets. In this example, LMA and LMR distributions are measured and compared in a rectangular enclosure with three different inlet/outlet configurations. A supply slot is fixed at the top of a right wall, and an exhaust slot is varied at the bottom-left (Case 1), bottom-right (Case 2), and top-left (Case 3) locations.

#### **6.2.2 Experimental setup**

The experimental chamber has dimensions of 1.8 m x 1.2 m x 0.9 m. There is a supply slot on the top of the right wall, and an exhaust slot at one of the three locations. The supply and exhaust slots are 0.025 m in width, and supply air was discharged horizontally. The airflow rate ranged from 4 to 76 ACH. The pressure inside the chamber was maintained neutral by an exhaust fan in order to minimize infiltration through the envelope. Sulfur hexafluoride at 30% concentration was used as a tracer gas. Using a syringe, 10 mL of SF6 was injected into a polystyrene tube, and the gas was mixed with a continuous stream of nitrogen. The diluted tracer gas was discharged at a point in the chamber through a porous sphere 40 mm in diameter connected at the end of the injection tube. A tracer gas detector is a multi-gas monitor based on the non-disperse infrared (NDIR) absorption principle. To visualize airflow patterns in the chamber, helium bubbles were discharged into a supply air duct, and a sheet of light was illuminated through a glass window along the center of the chamber. A schematic diagram of the experimental setup is shown in Fig. 10.

Fig. 10. Schematic diagram of experimental setup (Han *et al*., 2002).

#### **6.2.3 Experimental procedure**

56 Fluid Dynamics, Computational Modeling and Applications

Further research needs to be done to improve the tracer gas technique and to apply the

The airflow pattern in a ventilated space varies according to the locations of supply inlets and exhaust outlets. In this example, LMA and LMR distributions are measured and compared in a rectangular enclosure with three different inlet/outlet configurations. A supply slot is fixed at the top of a right wall, and an exhaust slot is varied at the bottom-left

The experimental chamber has dimensions of 1.8 m x 1.2 m x 0.9 m. There is a supply slot on the top of the right wall, and an exhaust slot at one of the three locations. The supply and exhaust slots are 0.025 m in width, and supply air was discharged horizontally. The airflow rate ranged from 4 to 76 ACH. The pressure inside the chamber was maintained neutral by an exhaust fan in order to minimize infiltration through the envelope. Sulfur hexafluoride at 30% concentration was used as a tracer gas. Using a syringe, 10 mL of SF6 was injected into a polystyrene tube, and the gas was mixed with a continuous stream of nitrogen. The diluted tracer gas was discharged at a point in the chamber through a porous sphere 40 mm in diameter connected at the end of the injection tube. A tracer gas detector is a multi-gas monitor based on the non-disperse infrared (NDIR) absorption principle. To visualize airflow patterns in the chamber, helium bubbles were discharged into a supply air duct, and a sheet of light was illuminated through a glass window along the center of the chamber. A

Halogen Lamp for Flow Visualization

Computer Digital Camera

Fig. 10. Schematic diagram of experimental setup (Han *et al*., 2002).

y P

x

Mircro - Manometer

Injection Point

He-bubble Generator

> Supply Plenum

Connector

Syringe

SF6 gas

Screen Nozzle

Mircro -

Flow meter

Supply Fan

Nitrogen Gas

**6.2 Effect of inlet/outlet configurations on LMA and LMR distributions** 

(Case 1), bottom-right (Case 2), and top-left (Case 3) locations.

schematic diagram of the experimental setup is shown in Fig. 10.

technique to various applications.

**6.2.1 Problem description** 

**6.2.2 Experimental setup** 

SF gas <sup>6</sup> Monitor

Exhaust Fan Mixing

Box

Manometer

To Outdoor

Computer

For the LMR measurements, tracer gas was injected through a porous sphere at a point in the chamber and the transient tracer gas concentration variation was measured at the exhaust. After the tracer gas was exhausted completely from the chamber, the injection sphere was moved to another position. The procedure was repeated for other internal points. There are 15 injection points equally spaced in the center plane of the chamber.

For LMA measurements, all the experimental conditions were identical to the LMR case, but tracer gas was injected at a supply duct. Then, transient tracer gas concentration variation was measured at the internal points. Experiments were conducted for three different exhaust locations under isothermal room temperature conditions.

#### **6.2.4 Results and discussion**

Flow visualization results are shown in Fig. 11 for three different exhaust locations. The air change rates are 12ACHs. For Case 1, the air supplied in the horizontal direction moved toward the lower-left exhaust in the diagonal direction. The room air formed two large recirculating flows at the upper-left and lower-right corners. For Case 2, the air supplied along the ceiling changed its direction by the opposite wall and made a large counterclockwise circulation in the chamber. We note that the airflow pattern was quite similar to a complete mixing condition. For Case 3, supplied air faced directly toward the exhaust. The room air was mostly stagnant, but with a slow recirculation due to the viscous action of the bypassing flow along the ceiling.

Fig. 11. Flow visualization results for three inlet-outlet configurations.

Contours of LMA and LMR are plotted in Fig. 12 for Case 1. It can be seen that both of the distributions are closely related to the airflow pattern shown in Fig. 11(a). LMA and LMR values are large within recirculating zones. LMA is small adjacent to the supply inlet and large adjacent to the exhaust, whereas LMR is small adjacent to the exhaust and large adjacent to the supply inlet. Figure 13 shows LMA and LMR distributions for Case 2. The LMA near ceiling is relatively small, whereas the LMR near floor is small. Both have large values within a large recirculation zone at the center. Figure 14 shows the results for Case 3. The LMA and LMR are small near ceiling, and large adjacent to the floor. The airflow pattern in Fig. 11(c) indicates there was large stagnant recirculation in the lower part of the space. The tracer gas diffused out into the lower part could not be exhausted effectively.

Fig. 15. Effect of air change rate on room mean ventilation effectiveness for three cases.

The distributions of LMA and LMR were obtained in a rectangular space with three different inlet and outlet configurations, and the corresponding airflow patterns were

1. The distributions of LMA and LMR show different characteristics, but both are closely

2. LMA values are small adjacent to supply inlets, and large adjacent to return-air exhausts. LMR values are small adjacent to exhausts, and large adjacent to supply air

3. Compared to Cases 1 and 2, Case 3 shows poor overall room ventilation effectiveness, since the supply air jet is directed toward the exhaust outlet located at the opposite side. 4. The overall ventilation effectiveness depends not only on supply-exhaust

The concept of local mean residual lifetime of the contaminant can be used in designing the layouts of exhausts and contaminant sources in a building such as a smoking zone, whereas concept of local mean age can be used in designing a proper distribution of fresh supply air

A space with multiple inlets is considered. It has a pentagonal shape with two inlets and a single outlet, which models a simplified livestock building. The LMAs from individual inlets are obtained by injecting a tracer gas at each inlet separately, and the combined LMA is obtained by injecting a tracer gas at both inlets simultaneously. This example is intended

The experimental chamber is pentagonal in shape with a height of 1.4 m and a width of 3.0 m. It is roughly a one-third scale model of a livestock building. The chamber has a length of 0.15

to verify the relation previously derived theoretically between the LMAs.

**6.2.5 Concluding remarks** 

inlets, as expected.

into an occupied zone.

**6.3.1 Problem description** 

**6.3.2 Experimental setup** 

related to the airflow pattern in the space.

configurations, but also on the air change rate.

**6.3 LMA distributions in a space with multiple inlets** 

visualized.

Fig. 12. LMA and LMR distributions for Case 1 (Han *et al*., 2002).

Fig. 13. LMA and LMR distributions for Case 2 (Han *et al*., 2002).

Fig. 14. LMA and LMR distributions for Case 3 (Han *et al.*, 2002).

Figure 15 shows room mean ventilation effectiveness for various air change rates. For Case 1, ventilation effectiveness decreased as the air change rate increased, but remained nearly constant for large air change rates over 20. It varies between 0.8 and 1.0, which is similar to a complete mixing condition. Note that the room ventilation effectiveness is 1 for complete mixing conditions, and 2 for perfect piston flow conditions. For Case 2, as the air change rate increased, the effectiveness increased initially and decreased slowly afterward. The effectiveness remained nearly constant for ACH over 20, similar to Case 1. However, the ventilation effectiveness in Case 3 is significantly lower compared to Cases 1 and 2, especially when the air change rate was low. This is due to the fact that the supply jet was not mixed well with the air in the chamber.

Fig. 15. Effect of air change rate on room mean ventilation effectiveness for three cases.

#### **6.2.5 Concluding remarks**

58 Fluid Dynamics, Computational Modeling and Applications

(a) LMA contour (s). (b) LMR contour (s).

(a) LMA contour (s). (b) LMR contour (s).

(a) LMA contour (s). (b) LMR contour (s).

Figure 15 shows room mean ventilation effectiveness for various air change rates. For Case 1, ventilation effectiveness decreased as the air change rate increased, but remained nearly constant for large air change rates over 20. It varies between 0.8 and 1.0, which is similar to a complete mixing condition. Note that the room ventilation effectiveness is 1 for complete mixing conditions, and 2 for perfect piston flow conditions. For Case 2, as the air change rate increased, the effectiveness increased initially and decreased slowly afterward. The effectiveness remained nearly constant for ACH over 20, similar to Case 1. However, the ventilation effectiveness in Case 3 is significantly lower compared to Cases 1 and 2, especially when the air change rate was low. This is due to the fact that the supply jet was

Fig. 12. LMA and LMR distributions for Case 1 (Han *et al*., 2002).

Fig. 13. LMA and LMR distributions for Case 2 (Han *et al*., 2002).

Fig. 14. LMA and LMR distributions for Case 3 (Han *et al.*, 2002).

not mixed well with the air in the chamber.

The distributions of LMA and LMR were obtained in a rectangular space with three different inlet and outlet configurations, and the corresponding airflow patterns were visualized.


The concept of local mean residual lifetime of the contaminant can be used in designing the layouts of exhausts and contaminant sources in a building such as a smoking zone, whereas concept of local mean age can be used in designing a proper distribution of fresh supply air into an occupied zone.

#### **6.3 LMA distributions in a space with multiple inlets 6.3.1 Problem description**

A space with multiple inlets is considered. It has a pentagonal shape with two inlets and a single outlet, which models a simplified livestock building. The LMAs from individual inlets are obtained by injecting a tracer gas at each inlet separately, and the combined LMA is obtained by injecting a tracer gas at both inlets simultaneously. This example is intended to verify the relation previously derived theoretically between the LMAs.

#### **6.3.2 Experimental setup**

The experimental chamber is pentagonal in shape with a height of 1.4 m and a width of 3.0 m. It is roughly a one-third scale model of a livestock building. The chamber has a length of 0.15

rate was 27 CMH. Each figure shows three injection cases: injection at vent *a* (Case *a*), at

The concentrations increased rapidly initially, and reached a constant steady state. The steady concentrations indicate the effective supply airflow rates contributing to the ventilation at the point by each supply inlet in a relative sense. The steady concentration in Case *b* is greater than that in Case *a*, which means the ventilation performance at point P was influenced more by the supply air from vent *b* than by the supply air from vent *a*. At the exhaust, the steady concentrations in Cases *a* and *b* are nearly the same, since the airflow rates of the two inlets are the same. We note that the non-dimensional steady concentrations at the exhaust could be determined by the relative airflow rates from the two supply inlets.

> Vent a Vent b Both

CO2

(a) At point P. (b) At exhaust. Fig. 17. Concentration responses at P and at exhaust after step-up injections at the inlets

concentration(ppm)

0 200 400 600 800

Contour Graph 3

0.16 0.53 0.76

0.17

0.05

X Data 0.0 0.5 1.0 1.5 2.0 2.5 3.0

Length[m]

**0.7**

0.15 0.46 0.66

0.01

0.09 0.01 0.66

**0.5**

**0.9**

**0.3**

**0.1**

0.02 0.12 0.76

0.01

Time(s)

Vent a Vent b Both *b ex c a ex c ba ex c* +

*b Pc a Pc ba Pc* +

vent *b* (Case *b*), and at vents *a* and *b* (Case *c*).

0 200 400 600 800

Contour Graph 6

0.98

0.47

0.98

0.95

0.98 0.54

0.99

0.91

Y Data

Height[m]

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4

(a) Injection at vent *a*. (b) Injection at vent *b*.

The steady concentrations were obtained by taking the averages of the fluctuating concentrations for a certain period of time after reaching the steady state. Figure 18 shows the spatial distributions of the steady concentrations measured at internal points. The concentration values have been made dimensionless by dividing those by the steady concentrations obtained in Case *c*. Iso-concentration contours were drawn based on the numerical values measured at the grid points. In Fig. 18(a), non-dimensional steady concentrations by vent *a* are greater than 0.5 at an upper part of the space, and less than 0.5

Col 5

0.02 0.07 0.92

0.99 0.34

0.34

Fig. 18. Spatial distributions of steady concentrations (Han *et al.,* 2011).

X Data 0.0 0.5 1.0 1.5 2.0 2.5 3.0

0.24

**0.5**

**0.7**

**0.9**

0.98

0.99

0.88 0.24

Length[m]

**0.3**

Time(s)

(Han *et al.,* 2011).

0.97

0.71 0.23

CO2

> Y Data

Height[m]

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4

Col 5

**0.1**

concentration(ppm)

m, and thus can be considered 2-D. There are five openings in the model, and three are used for our experiment. Two openings on the left wall (vents *a* and *b*) are used as supply inlets, and one opening on the opposite wall (vent *d*) is used as an exhaust. Vents that are not used for the experiment have been carefully sealed. The sizes of all openings are 0.05 m x 0.15 m. A schematic diagram of the experimental setup is shown in Fig. 16.

Fig. 16. Schematic diagram of experimental setup (Han *et al.,* 2011).

Carbon dioxide was used as a tracer gas. Injection ports were installed in both supply ducts upstream of the flow nozzles to ensure that the tracer gas was well mixed with incoming air streams. The amount of tracer gas was controlled by mass flow controllers (MFCs). A MFC contains a thermal mass flow meter that measures the air temperature rise across an internal heater. The range is between 0 and 10 L/min, and the error is reported to be below 1% of the measured values. A step-up method was adopted for tracer injection using a MFC. The tracer gas injection rate was held constant until a steady state condition was reached. The gas detector was an infrared single gas analyzer, and the sampling interval was 1.6 s. The range of the monitor was 20,000 ppm maximum, and the accuracy is 1% of the range.

#### **6.3.3 Experimental procedure**

Three cases of tracer injections were applied: injection at vent *a*, injection at vent *b*, and injection at vents *a* and *b* simultaneously. In order to obtain local mean age distributions in the space, tracer concentration responses were measured at 19 internal points evenly distributed in the space, including point P. The airflow rate was varied from 16.2 to 54 CMH. The airflow rates of vents *a* and *b* are maintained to be the same.

#### **6.3.4 Results and discussion**

Figure 17 shows concentration responses measured at point P and at the exhaust. The concentrations have been obtained by subtracting the background concentration, which is the average concentration measured before a tracer injection is applied. The total airflow

m, and thus can be considered 2-D. There are five openings in the model, and three are used for our experiment. Two openings on the left wall (vents *a* and *b*) are used as supply inlets, and one opening on the opposite wall (vent *d*) is used as an exhaust. Vents that are not used for the experiment have been carefully sealed. The sizes of all openings are 0.05 m x 0.15 m. A

schematic diagram of the experimental setup is shown in Fig. 16.

Nozzle

ΔP

Flow straightener

Nozzle

ΔP

manometer

Fig. 16. Schematic diagram of experimental setup (Han *et al.,* 2011).

Supply duct

Transition

Vent *a*

**Supply**

1.5m

Carbon dioxide was used as a tracer gas. Injection ports were installed in both supply ducts upstream of the flow nozzles to ensure that the tracer gas was well mixed with incoming air streams. The amount of tracer gas was controlled by mass flow controllers (MFCs). A MFC contains a thermal mass flow meter that measures the air temperature rise across an internal heater. The range is between 0 and 10 L/min, and the error is reported to be below 1% of the measured values. A step-up method was adopted for tracer injection using a MFC. The tracer gas injection rate was held constant until a steady state condition was reached. The gas detector was an infrared single gas analyzer, and the sampling interval was 1.6 s. The

range of the monitor was 20,000 ppm maximum, and the accuracy is 1% of the range.

CMH. The airflow rates of vents *a* and *b* are maintained to be the same.

Three cases of tracer injections were applied: injection at vent *a*, injection at vent *b*, and injection at vents *a* and *b* simultaneously. In order to obtain local mean age distributions in the space, tracer concentration responses were measured at 19 internal points evenly distributed in the space, including point P. The airflow rate was varied from 16.2 to 54

Figure 17 shows concentration responses measured at point P and at the exhaust. The concentrations have been obtained by subtracting the background concentration, which is the average concentration measured before a tracer injection is applied. The total airflow

Straightener

Damper

Vent *e*

**P**

Vent *b* Vent *d*

3m

0.25m

Computer

1.4m

**Exhaust**

0.6m

0.05m 0.05m

0.05m

Vent *c*

CO2 gas monitor

Flange

Mass flow controller

Fan

Inverter Mico-

Injector

3 way valve

Tracer gas

**6.3.3 Experimental procedure** 

**6.3.4 Results and discussion** 

Trace gas tube rate was 27 CMH. Each figure shows three injection cases: injection at vent *a* (Case *a*), at vent *b* (Case *b*), and at vents *a* and *b* (Case *c*).

The concentrations increased rapidly initially, and reached a constant steady state. The steady concentrations indicate the effective supply airflow rates contributing to the ventilation at the point by each supply inlet in a relative sense. The steady concentration in Case *b* is greater than that in Case *a*, which means the ventilation performance at point P was influenced more by the supply air from vent *b* than by the supply air from vent *a*. At the exhaust, the steady concentrations in Cases *a* and *b* are nearly the same, since the airflow rates of the two inlets are the same. We note that the non-dimensional steady concentrations at the exhaust could be determined by the relative airflow rates from the two supply inlets.

Fig. 17. Concentration responses at P and at exhaust after step-up injections at the inlets (Han *et al.,* 2011).

Contour Graph 3

Fig. 18. Spatial distributions of steady concentrations (Han *et al.,* 2011).

The steady concentrations were obtained by taking the averages of the fluctuating concentrations for a certain period of time after reaching the steady state. Figure 18 shows the spatial distributions of the steady concentrations measured at internal points. The concentration values have been made dimensionless by dividing those by the steady concentrations obtained in Case *c*. Iso-concentration contours were drawn based on the numerical values measured at the grid points. In Fig. 18(a), non-dimensional steady concentrations by vent *a* are greater than 0.5 at an upper part of the space, and less than 0.5

through the corresponding inlets. The longer the individual LMA, the longer the corresponding supply air resides in the space. The combined LMAs are in the midst of individual LMAs, all of which vary linearly with respect to the nominal time constant. The weighted averages calculated from the individual LMAs are also shown in the figure. Notice that the weighting factors at exhaust are both 0.5 in this case, since the airflow rates are the same for both inlets. Theoretically, the combined LMA should be the same with the nominal time constant regardless of the airflow rates, which is shown by a solid line in the figure. We

note that the combined LMAs appeared within the 10% error band.

LMAa LMAb LMAa+b Eq. *<sup>c</sup> LMAP <sup>b</sup> LMAP <sup>a</sup> LMAP*

*Eqn* )12(.

LMA[s]

(a) LMA at P. (b) LMA at exhaust.

Fig. 20. Local mean ages of supply air at point P and at the exhaust as a function of nominal

In this example, a case of multiple inlets was considered. The relations between LMA values from individual inlets and the combined LMA were obtained experimentally in a simplified model space simulating livestock applications. The following conclusions are drawn from

1. Our experimental results confirmed the theoretical relation between the individual LMAs and the combined LMA of the total supply air. The weighting factors are the steady concentrations obtained with a continuous step-up tracer injection at the

2. At every point in the space, the non-dimensional steady concentrations are complimentary to each other. The non-dimensional steady concentration at a point can be considered as a relative contribution factor of an individual inlet to the supply

3. The spatial distribution of an individual LMA indicates how fast the supply air from the corresponding inlet can reach the space, and it is closely related to the airflow pattern in

4. These experimental procedures were verified by the fact that the overall local mean

The concepts and the relations developed in this study can be applied to various

ages at the exhaust are in good agreement with the nominal time constants.

applications to quantify supply characteristics of individual inlets.

0 30 60 90 120 150 180

Nominal time constant[s]

LMAa LMAb LMAa+b Eq. *<sup>c</sup> LMAex <sup>b</sup> LMAex <sup>a</sup> LMAex*

*Eqn* )12(.

±10%

0 40 80 120 160

Nominal time constant[s]

0

time constant.

these results.

the space.

**6.3.5 Concluding remarks** 

corresponding supply inlets.

characteristics at the point.

40

80

120

LMA[s]

160

200

at a lower part of the space. The concentration distributions by vent *b* are the opposite, as shown in Fig. 18(b). The non-dimensional steady concentrations are complimentary to each other; i.e., the sum of the concentrations is unity at any point in the space.

LMA contours from individual inlets are shown in Fig. 19. Figure 19(a) shows small LMAs in the vicinity of vent *a*, and large LMAs near vent *b* and at the upper-right corner. Figure 19(b) shows small values starting from vent *b* along the floor up to the exit on the right, and large values at three corners in the upper part of the space. By following the contour lines, we can visualize the approximate airflow pattern in the space directed toward the exit. Figure 19(c) shows the combined LMA by simultaneous injections at both supply inlets. The LMAs are small along the floor near vent *a*, and along the left part of the roof near vent *b*. The combined LMA can be calculated from the individual LMAs according to Eq. (12). The distribution is shown in Fig. 19(d) and can be compared to Fig. 19(c). The overall patterns are in good agreement.

Fig. 19. Spatial distributions of local mean ages for Cases *a*, *b*, and *c,* with weighted averages (Han *et al.,* 2011).

The local mean ages at P are shown in Fig. 20(a) for various airflow rates. The airflow rate is expressed with the nominal time constant, which is the inverse of the air change rate. As the nominal time increased, both *<sup>a</sup> LMAP* and *<sup>b</sup> LMAP* increased linearly. The slope of *<sup>a</sup> LMAP* is greater than that of *<sup>b</sup> LMAP* . The combined LMA by total supply air, *<sup>c</sup> LMAP* , falls between the two sets. The figure also shows the LMA data calculated from the individual LMAs using Eq. (12).

The LMA values at the exhaust are expressed with respect to the nominal time constant in Fig. 20(b). The individual LMAs at exhaust indicate the residence time of the air supplied through the corresponding inlets. The longer the individual LMA, the longer the corresponding supply air resides in the space. The combined LMAs are in the midst of individual LMAs, all of which vary linearly with respect to the nominal time constant. The weighted averages calculated from the individual LMAs are also shown in the figure. Notice that the weighting factors at exhaust are both 0.5 in this case, since the airflow rates are the same for both inlets. Theoretically, the combined LMA should be the same with the nominal time constant regardless of the airflow rates, which is shown by a solid line in the figure. We note that the combined LMAs appeared within the 10% error band.

Fig. 20. Local mean ages of supply air at point P and at the exhaust as a function of nominal time constant.

#### **6.3.5 Concluding remarks**

62 Fluid Dynamics, Computational Modeling and Applications

at a lower part of the space. The concentration distributions by vent *b* are the opposite, as shown in Fig. 18(b). The non-dimensional steady concentrations are complimentary to each

LMA contours from individual inlets are shown in Fig. 19. Figure 19(a) shows small LMAs in the vicinity of vent *a*, and large LMAs near vent *b* and at the upper-right corner. Figure 19(b) shows small values starting from vent *b* along the floor up to the exit on the right, and large values at three corners in the upper part of the space. By following the contour lines, we can visualize the approximate airflow pattern in the space directed toward the exit. Figure 19(c) shows the combined LMA by simultaneous injections at both supply inlets. The LMAs are small along the floor near vent *a*, and along the left part of the roof near vent *b*. The combined LMA can be calculated from the individual LMAs according to Eq. (12). The distribution is shown in Fig. 19(d) and can be compared to Fig. 19(c). The overall patterns

other; i.e., the sum of the concentrations is unity at any point in the space.

are in good agreement.

219 117 16

134 121 33

**60**

84 70 52

55

30

**30**

30

26 81 21

(Han *et al.,* 2011).

using Eq. (12).

**180**

**90**

129

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4

> 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4

Height[m]

Height[m]

0.0 0.5 1.0 1.5 2.0 2.5 3.0

0.0 0.5 1.0 1.5 2.0 2.5 3.0

85 77 55

51

51

**90**

Length[m]

**60**

**90**

129 62 39

159

**180**

265

163 121 43

**180**

150

228 212 38

**90**

93 72 69

71 93 56

85

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4

(c) Simultaneous injection at vents *a* and *b.* (d) Weighted average of (a) and (b).

Fig. 19. Spatial distributions of local mean ages for Cases *a*, *b*, and *c,* with weighted averages

The local mean ages at P are shown in Fig. 20(a) for various airflow rates. The airflow rate is expressed with the nominal time constant, which is the inverse of the air change rate. As the nominal time increased, both *<sup>a</sup> LMAP* and *<sup>b</sup> LMAP* increased linearly. The slope of *<sup>a</sup> LMAP* is greater than that of *<sup>b</sup> LMAP* . The combined LMA by total supply air, *<sup>c</sup> LMAP* , falls between the two sets. The figure also shows the LMA data calculated from the individual LMAs

The LMA values at the exhaust are expressed with respect to the nominal time constant in Fig. 20(b). The individual LMAs at exhaust indicate the residence time of the air supplied

Height[m]

(a) Injection at vent *a*. (b) Injection at vent *b*.

Height[m]

27 63 119 **30**

49 68 127

**120**

30

**60**

79 74 55

51

31 52 43

**30**

35

0.0 0.5 1.0 1.5 2.0 2.5 3.0

**90**

79 102 88

**90**

**60**

43

46

93 84 111

95 101 66

**120**

77

**90**

**120**

76

151 95 90

158 97 56

Length[m]

0.0 0.5 1.0 1.5 2.0 2.5 3.0

80 81 51

**90**

46

57

Length[m]

Length[m]

In this example, a case of multiple inlets was considered. The relations between LMA values from individual inlets and the combined LMA were obtained experimentally in a simplified model space simulating livestock applications. The following conclusions are drawn from these results.


The concepts and the relations developed in this study can be applied to various applications to quantify supply characteristics of individual inlets.

( )1 *C C v D m m* ⋅∇ =∇⋅ ∇ +

Therefore, the steady concentration divided by the source strength equals the local mean age

( )

0

Refrigerating, and Air-Conditioning Engineers, Atlanta, USA.

and Ventilation Centre, Coventry, United Kingdom.

(0) *C C dt C m*

AIVC, (1990). A Guide to Air Change Efficiency, Technical Note AIVC28, Air Infiltration

ASHRAE, (2009). *ASHRAE Handbook-Fundamentals*, American Society of Heating,

ASTM, (1993). Standard Test Methods for Determining Air Change in a Single Zone by

Corpron, M. H. (1992). *Design and Characterization of a Ventilation Chamber*, M.S. Thesis.

Danckwerts, P. V. (1958). Local Residence-Times in Continuous-Flow Systems, *Chemical* 

Dufton, A. F. & Marley, W. G. (1935). Measurement of Rate of Air Change, Institution of

Etheridge, D. & Sandberg, M. (1996). *Building Ventilation: Theory and Measurement*, John

Han, H. (1992). Calculation of Ventilation Effectiveness Using Steady-State Concentration

Han, H.; Kuehn, T. H. & Kim, Y. I. (1999). Local Mean Age Measurements for Heating,

Han, H.; Choi, S. H. & Lee, W. W. (2002). Distribution of Local Supply and Exhaust

Han, H.; Shin, C. Y.; Lee, I.B. & Kwon, K. S. (2010). Local Mean Ages of Air in a Room with Multiple Inlets, *Int. J. of Air-Conditioning and Refrigeration*, Vol. 18, No. 1, pp. 15-21. Han, H.; Shin, C. Y.; Lee, I.B. & Kwon, K. S. (2011). Tracer Gas Experiment for Local Mean

Kuehn, T. H.; Ramsey, J. W. & Threlkeld, J. L. (1998). *Thermal Environmental Engineering*, 3rd

Liang, H. (1994). Room Air Movement and Contaminant Transport, *Ph.D. Thesis*. Department of Mechanical Engineering, University of Minnesota, Minneapolis, USA. Sandberg, M. (1981). What is Ventilation Efficiency, *Building and Environment*, Vol. 16, No. 2,

Distributions and Turbulent Airflow Patterns in a Half Scale Office Building, Proc. of Int.l Symp. on Room Air Convection and Ventilation Effectiveness, pp. 187-191,

Cooling, and Isothermal Supply Air Conditions, *ASHRAE Trans*., Vol. 105, Pt. 2, pp.

Effectiveness according to Room Airflow Patterns, *International Journal of Air-*

Ages of Air from Individual Supply Inlets in a Space with Multiple Inlets, *Building* 

Means of a Tracer Gas Dilution, E741-93, American Society for Testing and

Department of Mechanical Engineering, University of Minnesota, Minneapolis, USA.

in the space:

**9. References** 

Materials, USA.

Tokyo, Japan.

275-282.

pp. 123-135.

*Engineering Science*, Vol. 9, pp. 78-79.

Wiley & Sons, New York, USA.

Heating and Ventilating Engineers, Vol. 1, p. 645.

*conditioning and Refrigeration,* Vol. 10, No. 4, pp. 177-183.

*and Environment*, Vol. 46, pp. 2462-2471.

ed., Prentice Hall, London, United Kingdom.

(A4)

<sup>∞</sup> <sup>∞</sup> ∴ = (A5)

#### **7. Conclusion**

The purpose of ventilation is to supply fresh air to an occupied space and to effectively remove contaminants generated within the space. Ventilation performance is determined not only by the air change rate, but also by the ventilation effectiveness.

This study dealt with ventilation effectiveness based on the concept of the age of air. Ventilation effectiveness was categorized into supply effectiveness and exhaust effectiveness. The local supply index was represented by the local mean age of supply air; similarly, the local exhaust index was represented by local mean residual-life-time of contaminant. Overall room ventilation effectiveness was expressed as one value, regardless of supply and exhaust, because the room average of the local supply index was found to be identical to that of the local exhaust index.

The age concept has been extended to a space with multiple inlet and outlet openings. Theoretical derivations were made to obtain the relations between the LMAs from individual inlets and the combined LMA of total supply air, as well as the relations between the LMRs toward individual outlets and the overall LMR of the total exhaust air. Those relations can be used to investigate the effect of each supply inlet among many inlets, and the contribution of each exhaust outlet among many outlets in a space with multiple inlets and outlets.

The tracer gas technique provided a powerful tool in our ventilation studies for measuring the ventilation effectiveness of a conditioned space as well as to evaluate the performance of diffusers and exhaust grills. The ventilation theories provided in this chapter can be applied to various applications to provide good indoor air quality and to save ventilation energy use in buildings.

#### **8. Appendix**

It can be easily proved that the local mean age distribution in a space is equivalent to the steady concentration distribution with uniformly distributed sources of unit strength in the space.

The general equation that governs the transient concentration distribution can be expressed as

$$\frac{\partial \mathbf{C}}{\partial t} + \vec{v} \cdot \nabla \mathbf{C} = \nabla \cdot (D \nabla \mathbf{C}) + \dot{m} \tag{A1}$$

where *D* is the diffusion coefficient of the contaminant in air. Consider the case of a stepdown procedure with no contaminant source in the space. By integrating Eq. (A1) from zero to infinity with *m* equal to zero, we obtain

$$\text{C}(\rightsquigarrow)-\text{C}(\text{O}) + \vec{v} \cdot \nabla \int\_{0}^{\rightsquigarrow} \text{C}dt = \nabla \cdot \left(\text{D}\nabla \int\_{0}^{\rightsquigarrow} \text{C}dt\right) \tag{A2}$$

The steady concentration is zero; thus, Eq. (A2) can be rewritten as

$$\left[\vec{v}\cdot\nabla\left[\int\_{0}^{\infty}\frac{\mathsf{C}}{\mathsf{C}(\mathsf{O})}dt\right] = \nabla\cdot\left(D\nabla\left[\int\_{0}^{\infty}\frac{\mathsf{C}}{\mathsf{C}(\mathsf{O})}dt\right]\right) + 1\tag{A3}$$

The expression in the bracket is the local mean age under a step-down procedure.

On the other hand, Eq. (A1) can be simplified for steady concentration with uniformly distributed sources. As *m* is constant through the space, the equation can be simplified as

$$
\vec{v} \cdot \nabla \frac{\mathbf{C}}{\dot{m}} = \nabla \cdot (D \nabla \frac{\mathbf{C}}{\dot{m}}) + \mathbf{1} \tag{A4}
$$

Therefore, the steady concentration divided by the source strength equals the local mean age in the space:

$$\dots \int\_0^\infty \frac{\mathbb{C}}{\mathbb{C}(0)} dt = \frac{\overline{\mathbb{C}}(\Leftrightarrow)}{\dot{m}} \tag{A5}$$

#### **9. References**

64 Fluid Dynamics, Computational Modeling and Applications

The purpose of ventilation is to supply fresh air to an occupied space and to effectively remove contaminants generated within the space. Ventilation performance is determined

This study dealt with ventilation effectiveness based on the concept of the age of air. Ventilation effectiveness was categorized into supply effectiveness and exhaust effectiveness. The local supply index was represented by the local mean age of supply air; similarly, the local exhaust index was represented by local mean residual-life-time of contaminant. Overall room ventilation effectiveness was expressed as one value, regardless of supply and exhaust, because the room average of the local supply index was found to be

The age concept has been extended to a space with multiple inlet and outlet openings. Theoretical derivations were made to obtain the relations between the LMAs from individual inlets and the combined LMA of total supply air, as well as the relations between the LMRs toward individual outlets and the overall LMR of the total exhaust air. Those relations can be used to investigate the effect of each supply inlet among many inlets, and the contribution of each exhaust outlet among many outlets in a space with

The tracer gas technique provided a powerful tool in our ventilation studies for measuring the ventilation effectiveness of a conditioned space as well as to evaluate the performance of diffusers and exhaust grills. The ventilation theories provided in this chapter can be applied to various applications to provide good indoor air quality and to

It can be easily proved that the local mean age distribution in a space is equivalent to the steady concentration distribution with uniformly distributed sources of unit strength in

The general equation that governs the transient concentration distribution can be expressed as

( ) *<sup>C</sup> v C DC m*

where *D* is the diffusion coefficient of the contaminant in air. Consider the case of a stepdown procedure with no contaminant source in the space. By integrating Eq. (A1) from zero

*C C v Cdt D Cdt* ( ) (0) ( ) ∞ ∞

The expression in the bracket is the local mean age under a step-down procedure.

0 0 ( ) <sup>1</sup> (0) (0) *C C v dtD dt C C* ∞ ∞ ⋅ ∇ =∇⋅ ∇ + 

On the other hand, Eq. (A1) can be simplified for steady concentration with uniformly distributed sources. As *m* is constant through the space, the equation can be simplified as

0 0

(A1)

∞ − + ⋅∇ =∇⋅ ∇ (A2)

(A3)

<sup>∂</sup> + ⋅∇ =∇⋅ ∇ +

*t*

∂

The steady concentration is zero; thus, Eq. (A2) can be rewritten as

not only by the air change rate, but also by the ventilation effectiveness.

identical to that of the local exhaust index.

save ventilation energy use in buildings.

to infinity with *m* equal to zero, we obtain

multiple inlets and outlets.

**8. Appendix** 

the space.

**7. Conclusion** 


**4** 

*Slovenia* 

**Fluid Dynamic Models Application** 

Peter Vidmar, Stojan Petelin and Marko Perkovič

*University of Ljubljana, Faculty of Maritime Studies and Transport* 

Risk is a common name for the probability of a hazard turning into a disaster. Vulnerability and hazard are not dangerous in and of themselves, but if they come together, they generate a risk. However, risk can be reduced and managed. If we are careful about how we treat the environment, and if we are aware of our weaknesses and vulnerabilities to existing hazards,

Hazard from LNG (Liquefied Natural Gas) cargo begins in the first processing stage of natural gas liquefaction and loading the substance into LNG tankers. The transport itself over the sea is the safest part of the distribution process, as is demonstrated by the statistic of nautical accidents in the past 40 years (DNV, 2007, Perkovic et al., 2010 & Gucma, 2007). A review of a Rand Corporation document (Murray et al.) published in 1976 indicates a high level of safety for LNG tankers. The document indicates that in the initial 16-year history (from 1959 up to 1974) there had been no significant accidents. It should be noted, though, that in 1974 the world LNG fleet included only 14 vessels; by November, 2009, there were 327 vessels,

The DNV (Det Norske Veritas) counts 185 nautical accidents involving LNG tankers, all without severe consequences for the crew. The frequency of LNG tanker accidents is therefore 5.6 x 10-2 per ship year. The findings of the DNV (2007) furthermore demonstrate that the potential loss of life for the LNG crew member is 9.74 x 10-3 or less, considering the occupational fatality rate onboard gas tankers is 4.9 x 10-4. The analysis of the northern Adriatic Sea (Petelin et al. 2009), or, precisely, the gulf of Trieste, demonstrates that nautical accidents should occur with a frequency of 1.25 x 10-2 per year, assuming the current traffic

The hazard associated with LNG is mainly in its potential to cause severe fires resulting in heat radiation. If a large quantity of LNG is spilled into a pool, the cloud that is formed as it evaporates is a mixture of natural gas, water vapour, and air. Initially the cloud is heavier than air (due to its low storage temperature) and remains close to the ground. The buoyancy moves the natural gas upward at a gas temperature of around 170 K (-1030C), as experimentally demonstrated by ioMosaic (2007). The major influences on natural gas diffusion are environmental conditions. The cloud moves in the direction of the wind and the wind causes the cloud to mix with more air. If the concentration of gas in the air is between 5% and 15% it is flammable and burns if it contacts any ignition source. A concentration of gas smaller than 5% will not ignite and if the concentration is over 15% the air becomes saturated. The explosion of natural gas is not possible in open spaces because

then we can take measures to make sure that hazards do not turn into disasters.

a figure expected to increase to 350 vessels sometime in 2010 (LNG Journal, 2008).

density, and increases to 2.62 x 10-2 if the ship traffic increases by 100%.

**1. Introduction** 

**in Risk Assessment** 

