**Fluid Dynamic Models Application in Risk Assessment**

Peter Vidmar, Stojan Petelin and Marko Perkovič *University of Ljubljana, Faculty of Maritime Studies and Transport Slovenia* 

#### **1. Introduction**

66 Fluid Dynamics, Computational Modeling and Applications

Sandberg, M. & Sjoberg, M. (1983). The Use of Moments for Assessing Air Quality in Ventilated Rooms, *Building and Environment*, Vol. 18, No. 4, pp. 181-197. Shaw, C. Y., Zhang, J. S., Said, M. N. A., Vaculik, F. & Magee, R. J. (1992). Effect of Air

Spalding, D. B. (1958). A Note on Mean Residence-Times in Steady Flows of Arbitrary

Yaglou, C. P. & Witheridge W. N. (1937). Ventilation Requirements, *ASHVE Trans*., Vol. 42,

Xing, H., Hatton, A. & Awbi, H. B. (2001). A Study of the Air Quality in the Breathing Zone

in a Room with Displacement Ventilation, *Building and Environment*, Vol. 36, pp.

*Environmental International*, Vol. 12, Issues 1-4, pp. 419-427.

Complexity, *Chemical Engineering Science*, Vol. 9, pp. 74-77.

pp. 423-436.

809-820.

Diffuser Layout on the Ventilation Conditions of a Workstation-Part II: Air Change Efficiency and Ventilation Efficiency, *ASHRAE Trans.,* Vol. 99, Pt. 2, pp. 133-143. Skaaret, E. (1986). Contaminant Removal Performance in Terms of Ventilation Effectiveness,

> 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, then we can take measures to make sure that hazards do not turn into disasters.

> 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, a figure expected to increase to 350 vessels sometime in 2010 (LNG Journal, 2008).

> 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 density, and increases to 2.62 x 10-2 if the ship traffic increases by 100%.

> 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

Fluid Dynamic Models Application in Risk Assessment 69

Societal risk posed by an LNG terminal facility or hazardous activity is measured by the probability that a group of persons would be exposed to a hazardous level of harm (fatality) due to all types of accidents at the facility or its hazardous activity. The societal risk is dependent on both the density of people in the vicinity of an LNG terminal and the location of the population with respect to the facility. The societal risk is generally presented in the form of a curve, expressing the relationship between the annual probability (F) of exceeding a given number of fatalities and the number N (Trbojevic, 2005). In most countries the risk assessment is performed on the basis of potential fatalities to the exposed population. Different countries use slightly different criteria for risk acceptability. In the UK, the Health and Safety Executive (HSE) guidelines are to use the individual risk as the principal measure, but also to use the societal risk criteria for land use planning. The acceptability criteria levels for risks for facilities in the UK are specified by HSE (1989). Facilities are permitted only when these (published) criteria are met. In the Netherlands, however, both the individual risk criteria and the societal risk criteria have to be met when considering (in risk assessment) those events whose hazardous effects extend to such distances at which the conditional probability for lethality is higher than 1% (Bottelberghs, 2000 & Raj et al., 2009) . The risk tolerability criteria for fatalities established in various countries for societal risks are summarized in Fig. 1 (Trbojević, 2005).

The risk calculations have to determine the number of people killed by a particular event and attach this to the associated frequency of the event to form the F-N curve. Probability of death is calculated from the consequence model and so it only remains for the risk model to integrate the probability of death for each event over the specified population. *N f f* / *<sup>e</sup>* represents the number of people killed by a given event, particular considering the weather category/wind direction combination. It is calculated for the population on a grid according

is the individual risk as a sum of all events under all weather conditions.

*po x* int\_ ,*<sup>y</sup> IR* is the individual risk over all weather conditions for a specific event

*Pw* is the probability of specific weather and

Fig. 1. Comparison of F-N Criteria (Trbojević, 2005)

to Equation ( 4 ) (DNV, 2007).

where

*total x* \_ ,*<sup>y</sup> IR*

**2.2 Societal risk** 

the low velocity of flame spread, around 0.4 m/s, is not enough to produce a pressure wave (Fells et al, 1969). The burnout of gas/air mixture in open air could result only in a flash fire. Explosions could only result in enclosed spaces where gas is going to accumulate.

#### **2. Risk definition and acceptance criteria**

#### **2.1 Individual risk**

The individual risk in regard to LNG is calculated as location specific, to a person exposed outside 24 hours per day. In several countries the authorities have defined criteria which have to be met in order to assess a level of societal risk as acceptable. Criteria for some countries were analysed and discussed by Trbojević (2005). Some of these are used internationally. In most countries an individual risk of 1E10-3 per year is taken as the upper bound criteria to assess the acceptability of an activity for employers working in an industrial installation. The upper risk criterion for the public is therefore 1\*10-4 per year. Also the risk at locations where vulnerable objects are situated is taken into account. Vulnerable objects are those where people are present who react in a different way to a threat posed upon them as ordinary people. This difference can be caused by differences in state of health or the possibility of evacuating the location in case of danger. The calculation of the individual risk for a specific failure event is influenced by three main parameters: the two coordinate location *x,y* and weather conditions (wind speed and stability). The individual risk *<sup>x</sup>*, /*<sup>y</sup> <sup>w</sup> IR* is therefore the function of the frequency *Ffe* of the failure event

occurring in a year, the probability of the event in a particular direction (influenced by weather) and the probability of people dying due to exposure (DNV, 2008).

$$IR\_{x,y/w} = F\_{fe} \int\_{\theta\_1}^{\theta\_2} \left[ P\_{\theta/w} P\_{d/\theta w} \right] d\theta \tag{1}$$

θis the direction of the release,

θ*1* is the lower value of θthat influences the computation point,

θ*2* is the upper value of θthat influences the computation point,

*P*θ/*<sup>w</sup>* is the probability of the release occurring in the direction of the wind and

*Pd w* /θis the probability of death considering the direction of the release and weather.

This is the contribution to the individual risk from a single event under specific weather conditions. To obtain the overall individual risk at a given point all possible events must be taken into account. However in the sense of the order of magnitude the worst events (low probability/high consequences) have the most influence on individual risk. In the calculation process, low consequence events increase the size of high risk isolines, but usually do not influence the low risk isolines, risk lower than 1\*10-6 that reaches greatest distances. The strong dependence of the individual risk on weather conditions influences the total risk, calculated using the following equations:

$$IR\_{\text{point\\_x},y} = \sum\_{\text{Ven\\_there}} P\_w \cdot IR\_{x,y/w} \tag{2}$$

and

$$IR\_{\text{total\\_x\\_y}} = \sum\_{all\\_x, y} IR\_{pv\,\text{int\\_x}, y} \tag{3}$$


#### **2.2 Societal risk**

68 Fluid Dynamics, Computational Modeling and Applications

the low velocity of flame spread, around 0.4 m/s, is not enough to produce a pressure wave (Fells et al, 1969). The burnout of gas/air mixture in open air could result only in a flash fire.

The individual risk in regard to LNG is calculated as location specific, to a person exposed outside 24 hours per day. In several countries the authorities have defined criteria which have to be met in order to assess a level of societal risk as acceptable. Criteria for some countries were analysed and discussed by Trbojević (2005). Some of these are used internationally. In most countries an individual risk of 1E10-3 per year is taken as the upper bound criteria to assess the acceptability of an activity for employers working in an industrial installation. The upper risk criterion for the public is therefore 1\*10-4 per year. Also the risk at locations where vulnerable objects are situated is taken into account. Vulnerable objects are those where people are present who react in a different way to a threat posed upon them as ordinary people. This difference can be caused by differences in state of health or the possibility of evacuating the location in case of danger. The calculation of the individual risk for a specific failure event is influenced by three main parameters: the two coordinate location *x,y* and weather conditions (wind speed and stability). The individual risk *<sup>x</sup>*, /*<sup>y</sup> <sup>w</sup> IR* is therefore the function of the frequency *Ffe* of the failure event occurring in a year, the probability of the event in a particular direction (influenced by

Explosions could only result in enclosed spaces where gas is going to accumulate.

weather) and the probability of people dying due to exposure (DNV, 2008).

θ

θ

the total risk, calculated using the following equations:

2

1 *x y w fe w d w* , / / / *IR F P P d* θ

 is the probability of death considering the direction of the release and weather. This is the contribution to the individual risk from a single event under specific weather conditions. To obtain the overall individual risk at a given point all possible events must be taken into account. However in the sense of the order of magnitude the worst events (low probability/high consequences) have the most influence on individual risk. In the calculation process, low consequence events increase the size of high risk isolines, but usually do not influence the low risk isolines, risk lower than 1\*10-6 that reaches greatest distances. The strong dependence of the individual risk on weather conditions influences

> *po x* int\_ ,*y w x*, /*y w Weathers*

\_ , int\_ , \_ , *total x y po x y all x y*

θ

/*<sup>w</sup>* is the probability of the release occurring in the direction of the wind and

θ

that influences the computation point,

that influences the computation point,

=

 θ θ

*IR* = ⋅ *P IR* (2)

*IR* <sup>=</sup> *IR* (3)

(1)

**2. Risk definition and acceptance criteria** 

is the direction of the release,

*1* is the lower value of

*2* is the upper value of

**2.1 Individual risk** 

θ

θ

θ

*P*θ

*Pd w* /θ

and

Societal risk posed by an LNG terminal facility or hazardous activity is measured by the probability that a group of persons would be exposed to a hazardous level of harm (fatality) due to all types of accidents at the facility or its hazardous activity. The societal risk is dependent on both the density of people in the vicinity of an LNG terminal and the location of the population with respect to the facility. The societal risk is generally presented in the form of a curve, expressing the relationship between the annual probability (F) of exceeding a given number of fatalities and the number N (Trbojevic, 2005). In most countries the risk assessment is performed on the basis of potential fatalities to the exposed population. Different countries use slightly different criteria for risk acceptability. In the UK, the Health and Safety Executive (HSE) guidelines are to use the individual risk as the principal measure, but also to use the societal risk criteria for land use planning. The acceptability criteria levels for risks for facilities in the UK are specified by HSE (1989). Facilities are permitted only when these (published) criteria are met. In the Netherlands, however, both the individual risk criteria and the societal risk criteria have to be met when considering (in risk assessment) those events whose hazardous effects extend to such distances at which the conditional probability for lethality is higher than 1% (Bottelberghs, 2000 & Raj et al., 2009) . The risk tolerability criteria for fatalities established in various countries for societal risks are summarized in Fig. 1 (Trbojević, 2005).

Fig. 1. Comparison of F-N Criteria (Trbojević, 2005)

The risk calculations have to determine the number of people killed by a particular event and attach this to the associated frequency of the event to form the F-N curve. Probability of death is calculated from the consequence model and so it only remains for the risk model to integrate the probability of death for each event over the specified population. *N f f* / *<sup>e</sup>* represents the number of people killed by a given event, particular considering the weather category/wind direction combination. It is calculated for the population on a grid according to Equation ( 4 ) (DNV, 2007).

Fluid Dynamic Models Application in Risk Assessment 71

before it ruptures. The characteristic of this material is that it remains elastic at a cargo storage temperature of -1620C. The second fact is that the reservoir is filled up to 96-97% and the rest is gas, about 800 m3. Because of this gas space the containment structure allows slow deformations that occur during an accident. The accident of the El Paso Paul Kayser (Bubbico et al., 2009) is a good example. The accident caused significant deformations of the ship's hull and ribbed construction yet the deformation of the cargo containment was limited to one meter and there was no release. Pitpablo (2004) suggests that the maximum possible rupture size is 250 mm in diameter when caused by

Research conducted by DNV based on accident statistics and ship damage are based on the IMO (International Maritime Organisation) Guideline for Alternative Tanker Design (IMO, 1995). DNV found that an LNG tanker could resist a hull deformation of 3 meters and the reservoirs could resist 2 additional meters of deformation. Considered and analyzed ruptures on LNG vessels are only caused by accidents from the traffic point of view. Terroristic attacks and attacks with weapons are not considered. However, the interpretation of DNV (2008) and Sandia National Laboratory (Hightower et al., 2006) is widely acceptable. Major damage to LNG tankers by weapons is less probable because the vessel's structural stiffness is much greater than any building, bridge or other land structure. In case a projectile breaks through the primary reservoir, there is a high probability of immediate ignition, a local fire or even the destruction of the ship, but not the formation of a flammable gas cloud and subsequent flash fire that is a danger to the neighbour population. On this basis DNV suggests the consideration of the largest rupture size of 1500 mm for a

Therefore the further scenarios analysis would consider the following LNG reservoir

• 10,000 m3/h Maximum possible leakage rate caused by sabotage for 60 min

As mentioned above, the spreading of the LNG pool is a dynamic process that is alimented by the spilling from the tank slowed down by evaporation or even vapours burning. Research conducted by Hissong (2007) and ioMosaic (Melhem, 2007) presented in Table 1

> Evaporation rate 20 sec. avg. kg/m2s

• 250 mm Maximum possible rupture size caused by grounding • 750 mm Maximum possible rupture size caused by collision • 1500 mm Maximum possible rupture size caused by terrorist attack

• 7,000 m3/h Maximum possible leakage rate for 10 min

explain that the evaporation rate could vary depending on water conditions.

Max value kg/m2s

Water 0.303 0.147 Water 0.245 0.191 Water 0.230 0.196 Ice 0.513 0.191 Ice 0.333 0.171 Brine 0.254 0.186

Evaporation rate

Table 1. Evaporation rates depending on water conditions

**3.3 Atmospheric dispersion of a vapour cloud** 

grounding or collision with the shore.

terroristic attack.

rupture size:

$$N\_{f/f^\circ} = \iint\limits\_{area} P\_{d,xy} n\_{xy} dx dy \tag{4}$$

#### **3. Consequences analyses background**

Different scenarios of LNG release on water surface are possible, depending on where a rupture occurred. Accident scenarios like a collision between ships, grounding, and collision at the coast could all lead to a rupture of primary cargo containment that could spill LNG into the environment. The collision of two ships at 900, in which the side collapsed ship is the LNG carrier is the worst scenario. The HAZID model of this is presented by (Pitblado, 2004). The scenarios of grounding and collision with the shore are less critical because the impact force is distributed over a wider area of the ship and the penetration is not as deep in the vessel's structure. During the unloading operation the transfer arm could break because of an unpredictable ship movement, even it is moored. The shutdown of the transfer operation occurs in 1 to 10 minutes depending on the emergency system used. Transfer lines with EMS (Emergency Shut Down System) could interrupt the transfer in 1 minute. During the calculation a tolerance of an additional minute is added.

The investigation conducted indicates that the consequence approach should be based on some conservative initial conditions. We found that the best idea is to assume the greatest possible rupture size on primary containment or the transfer pipe, resulting from an accident. In this way a separate investigation could be conducted. The first is the calculation or assumption of possible rupture sizes and the second is the analyses of the resulting spills dynamics. In the paper the rupture size is assumed and taken from the literature (Pitblado, 2004), but the spill dynamics are computed.

#### **3.1 Analyses approaches**

The computation of potential consequences is the first step toward understanding the severity of an accident. Several commercial computer programs apply a lumped parameter modelling technique. Empirical equations derived from physical equations consider different release scenarios, meteorological conditions, obstacle density, etc. A zone modelling technique divides the physical space into zones. Within each zone the uniform physical phenomena are computed using physical and empirical equations. With this approach several approximations and simplifications within models are assumed to reduce the complexity of formulations and to reduce the computation time. In some applications, like gas dispersion dynamics, results from empirical models could lead to underestimated results. For this reason the empirical models usually apply corrective parameters or factors to compute more conservative, corrective, results. Overestimated consequences results could lead to overestimated risk results and falsely estimate the risk level to exceed allowed limits (Fig. 1). When this is suspected, further consequences analysis should be conducted, applying more accurate methods. The case of gas dispersion is computed using a CFD (Computational Fluid Dynamic) code and analysed for improved results.

#### **3.2 Rupture size**

The estimation of the rupture size caused by an accident involving an LNG carrier is not easy, because of the variables: a complex structure, the type of accident and the location of the primary rupture. The LNG carrier has four to five barriers that would have to fail before the LNG cargo is spilled. The shield of the LNG containment allows deformation

*N P n dxdy* = (4)

*f fe d xy xy* / , *area*

Different scenarios of LNG release on water surface are possible, depending on where a rupture occurred. Accident scenarios like a collision between ships, grounding, and collision at the coast could all lead to a rupture of primary cargo containment that could spill LNG into the environment. The collision of two ships at 900, in which the side collapsed ship is the LNG carrier is the worst scenario. The HAZID model of this is presented by (Pitblado, 2004). The scenarios of grounding and collision with the shore are less critical because the impact force is distributed over a wider area of the ship and the penetration is not as deep in the vessel's structure. During the unloading operation the transfer arm could break because of an unpredictable ship movement, even it is moored. The shutdown of the transfer operation occurs in 1 to 10 minutes depending on the emergency system used. Transfer lines with EMS (Emergency Shut Down System) could interrupt the transfer in 1 minute. During

The investigation conducted indicates that the consequence approach should be based on some conservative initial conditions. We found that the best idea is to assume the greatest possible rupture size on primary containment or the transfer pipe, resulting from an accident. In this way a separate investigation could be conducted. The first is the calculation or assumption of possible rupture sizes and the second is the analyses of the resulting spills dynamics. In the paper the rupture size is assumed and taken from the literature (Pitblado,

The computation of potential consequences is the first step toward understanding the severity of an accident. Several commercial computer programs apply a lumped parameter modelling technique. Empirical equations derived from physical equations consider different release scenarios, meteorological conditions, obstacle density, etc. A zone modelling technique divides the physical space into zones. Within each zone the uniform physical phenomena are computed using physical and empirical equations. With this approach several approximations and simplifications within models are assumed to reduce the complexity of formulations and to reduce the computation time. In some applications, like gas dispersion dynamics, results from empirical models could lead to underestimated results. For this reason the empirical models usually apply corrective parameters or factors to compute more conservative, corrective, results. Overestimated consequences results could lead to overestimated risk results and falsely estimate the risk level to exceed allowed limits (Fig. 1). When this is suspected, further consequences analysis should be conducted, applying more accurate methods. The case of gas dispersion is computed using a CFD

The estimation of the rupture size caused by an accident involving an LNG carrier is not easy, because of the variables: a complex structure, the type of accident and the location of the primary rupture. The LNG carrier has four to five barriers that would have to fail before the LNG cargo is spilled. The shield of the LNG containment allows deformation

(Computational Fluid Dynamic) code and analysed for improved results.

**3. Consequences analyses background** 

the calculation a tolerance of an additional minute is added.

2004), but the spill dynamics are computed.

**3.1 Analyses approaches** 

**3.2 Rupture size** 

before it ruptures. The characteristic of this material is that it remains elastic at a cargo storage temperature of -1620C. The second fact is that the reservoir is filled up to 96-97% and the rest is gas, about 800 m3. Because of this gas space the containment structure allows slow deformations that occur during an accident. The accident of the El Paso Paul Kayser (Bubbico et al., 2009) is a good example. The accident caused significant deformations of the ship's hull and ribbed construction yet the deformation of the cargo containment was limited to one meter and there was no release. Pitpablo (2004) suggests that the maximum possible rupture size is 250 mm in diameter when caused by grounding or collision with the shore.

Research conducted by DNV based on accident statistics and ship damage are based on the IMO (International Maritime Organisation) Guideline for Alternative Tanker Design (IMO, 1995). DNV found that an LNG tanker could resist a hull deformation of 3 meters and the reservoirs could resist 2 additional meters of deformation. Considered and analyzed ruptures on LNG vessels are only caused by accidents from the traffic point of view. Terroristic attacks and attacks with weapons are not considered. However, the interpretation of DNV (2008) and Sandia National Laboratory (Hightower et al., 2006) is widely acceptable. Major damage to LNG tankers by weapons is less probable because the vessel's structural stiffness is much greater than any building, bridge or other land structure. In case a projectile breaks through the primary reservoir, there is a high probability of immediate ignition, a local fire or even the destruction of the ship, but not the formation of a flammable gas cloud and subsequent flash fire that is a danger to the neighbour population. On this basis DNV suggests the consideration of the largest rupture size of 1500 mm for a terroristic attack.

Therefore the further scenarios analysis would consider the following LNG reservoir rupture size:


#### **3.3 Atmospheric dispersion of a vapour cloud**

As mentioned above, the spreading of the LNG pool is a dynamic process that is alimented by the spilling from the tank slowed down by evaporation or even vapours burning. Research conducted by Hissong (2007) and ioMosaic (Melhem, 2007) presented in Table 1 explain that the evaporation rate could vary depending on water conditions.


Table 1. Evaporation rates depending on water conditions

Fluid Dynamic Models Application in Risk Assessment 73

The first simulations conducted are used to compare results obtained from the lumped model (fast computation) and CFD model like FDS (Fire Dynamics Simulator). The scenarios include a simple geometry with a domain size 500m x 200m x 100 m in *x*, *y* and *z* directions. The diameter of the spilled pool is related with the diameter of the hole on a reservoir shield. Possible or applied scenarios are refined in a previous section, the resulting pool sizes are further obtained from functions explained in (ioMosaic, 2007). The following scenarios assume the hole at the bottom of the reservoir and 9000 kg of LNG content is spilled in 10 seconds and forms the pool. The evaporation rate of natural gas from the pool surface is calculated by the lumped model and assumed to be 0.16 kg/m2s for the CFD

Because of the very similar pool sizes scenarios 2 and 4 are calculated once as are scenarios 3

The course of events after an accident could follow different directions before the course of the accident is definitive. Depending on the accident dynamics, and the time and place of occurrence it could lead to different levels of risk. An example is the immediate ignition of a vapour cloud. If this occurs in a populated area, it represents a high risk event, but if it occurs outside the populated area it would be a low risk event. The reason is in the consequences of the event. The second event mentioned would be high risk if the ignition of the cloud would be delayed and the cloud would be transported downwind into a

The first scenario assumes the evaporation of natural gas from the LNG pool on a water surface of an area of 400 m2. The evaporating gas is moved downwind at 2 m/s; the surrounding temperature is 20oC. The model does not include complex geometry, buildings,

> Material ID Methane Specific heat 2.2 kJ/kg K Conductivity 0.03281 W/m K

Reference temperature -162 oC Absorption coefficient 40 1/m Boiling temperature -162 oC Fuel rate 1.0 Density (gas) 1.8 kg/m3 Gas rate 1.0 Mass flux max: 0.16 kg/m2s

Gas mass fraction 1.0

and relief that would slightly change the flow dynamics of the vapour cloud.

**4. Scenarios and results** 

and 5.

populated area.

model in a stable sea and stable weather conditions.

**4.1 Pool size related to rupture size** 

Table 2. Mandatory input data for Scenario 1:

1. D= 50 m Calculated pool diameter for the hole diameter 250 mm 2. D= 200 m Calculated pool diameter for the hole diameter 750 mm 3. D= 450 m Calculated pool diameter for the hole diameter 1500 mm 4. D= 196 m Calculated pool diameter for leaking 7,000 m3/h for 10 min 5. D = 400 m Calculated pool diameter for leaking 10,000 m3/h for 60 min

Leak rates from a tank rupture were calculated; initial gas, "flashing" and an aerosol are formed. An aerosol is a cloud of tiny liquids droplets or fine solid droplets suspended in air. Calculating the droplet evaporation along the cloud trajectory, the overall vapour generation rate is obtained. In this section the dispersion model of the vapour cloud is described. The model considers different types of release: instantaneous, steady continuous and transient for dense (active) and lean (passive) gases.

Our case considers a ground level instantaneous release from a tank rupture. That is the most catastrophic scenario, one which can lead to fast vapour cloud formation near ground level and with a high gas concentration. Because of a very fast transient and changes of variables, it is difficult to predict the course of events, especially close to the source of dispersion. Once the cloud, modelled as a cylinder, is formed, it begins to slump under the effect of gravity. The velocity of the edge of the cloud can be described as (Safer Sysetm, 1996):

$$\frac{dR}{dt} = k\_1 \left[ \left( \frac{\rho\_{\rm{CLOID}} - \rho\_A}{\rho\_{\rm{CLOID}}} \right) \cdot g \cdot h \right]^{1/2} \tag{5}$$

*k1* is a slumping constant that depends on the characteristic of released gases and weather conditions. The most important factor is wind speed. It is important to note that the dispersion model does not assume turbulent flow. Dispersion of the cloud is therefore a suitable element to be computed with the CFD model. The rate at which liquid fuel evaporates when burning is a function of the liquid temperature and the concentration of fuel vapour above the pool surface. According to the Clausius-Clapeyron relation, the volume fraction of the fuel vapour above the surface is a function of the liquid boiling temperature;

$$V\_f = \exp\left[-\frac{h\_v \mathcal{W}\_f}{R\_f} \left(\frac{1}{T\_s} - \frac{1}{T\_b}\right)\right] \tag{6}$$

where *<sup>v</sup> h* is the heat of vaporization, *Wf* is the molecular weight, *Ts* is the surface temperature, and *Tb* is the boiling temperature of the fuel.

The cloud transport due to wind is modelled as (Vidmar et al., 2003 & Safer System, 1996):

$$\frac{d\mathbf{x}}{dt} = \mathcal{U}\_{\text{CLOUD}} \left( \mathbf{z} = \mathbf{0}.4 \; h\_{\text{CLOUD}} \right), \tag{7}$$

where *x* is the downwind distance variable and U*CLOUD* is the cloud speed, assumed to be equal to the wind speed. The value0.4 times the cloud height is assumed to be a reference height or the centre of gravity of wind force in the direction of cloud movement.

In most dense gas releases, it is expected that there will be a central core region of uniform concentration along with edges at which the concentration decreases. It is expected that close to the source, the edges will be sharp and as the cloud disperses downwind the edges will become less steep. The concentration field is calculated considering a Gaussian distribution. That means that the isopleth limits, or the edge of observed concentration limits, take on a typical Gaussian distribution. Gaussian correlations for the atmospheric cloud dispersion are not proper for any initial conditions and release type, but are used as a conservative model (McGrattan, 1997).

#### **4. Scenarios and results**

72 Fluid Dynamics, Computational Modeling and Applications

Leak rates from a tank rupture were calculated; initial gas, "flashing" and an aerosol are formed. An aerosol is a cloud of tiny liquids droplets or fine solid droplets suspended in air. Calculating the droplet evaporation along the cloud trajectory, the overall vapour generation rate is obtained. In this section the dispersion model of the vapour cloud is described. The model considers different types of release: instantaneous, steady continuous

Our case considers a ground level instantaneous release from a tank rupture. That is the most catastrophic scenario, one which can lead to fast vapour cloud formation near ground level and with a high gas concentration. Because of a very fast transient and changes of variables, it is difficult to predict the course of events, especially close to the source of dispersion. Once the cloud, modelled as a cylinder, is formed, it begins to slump under the effect of gravity. The velocity of the edge of the cloud can be described as (Safer

> *CLOUD A CLOUD dR <sup>k</sup> g h dt* ρ

<sup>−</sup> <sup>=</sup> ⋅ ⋅

*k1* is a slumping constant that depends on the characteristic of released gases and weather conditions. The most important factor is wind speed. It is important to note that the dispersion model does not assume turbulent flow. Dispersion of the cloud is therefore a suitable element to be computed with the CFD model. The rate at which liquid fuel evaporates when burning is a function of the liquid temperature and the concentration of fuel vapour above the pool surface. According to the Clausius-Clapeyron relation, the volume fraction of the fuel vapour above the surface is a function of the liquid boiling

ρ

exp *v f*

*h W*

where *<sup>v</sup> h* is the heat of vaporization, *Wf* is the molecular weight, *Ts* is the surface

where *x* is the downwind distance variable and U*CLOUD* is the cloud speed, assumed to be equal to the wind speed. The value0.4 times the cloud height is assumed to be a reference

In most dense gas releases, it is expected that there will be a central core region of uniform concentration along with edges at which the concentration decreases. It is expected that close to the source, the edges will be sharp and as the cloud disperses downwind the edges will become less steep. The concentration field is calculated considering a Gaussian distribution. That means that the isopleth limits, or the edge of observed concentration limits, take on a typical Gaussian distribution. Gaussian correlations for the atmospheric cloud dispersion are not proper for any initial conditions and release type, but are used as a

( 0.4 ) *CLOUD CLOUD*

The cloud transport due to wind is modelled as (Vidmar et al., 2003 & Safer System, 1996):

*dx U zh*

height or the centre of gravity of wind force in the direction of cloud movement.

*f*

*V*

temperature, and *Tb* is the boiling temperature of the fuel.

conservative model (McGrattan, 1997).

 ρ

1 1

*dt* = = , (7)

*f sb*

*R TT* =− − ,

1/2

(5)

(6)

and transient for dense (active) and lean (passive) gases.

1

Sysetm, 1996):

temperature;

The first simulations conducted are used to compare results obtained from the lumped model (fast computation) and CFD model like FDS (Fire Dynamics Simulator). The scenarios include a simple geometry with a domain size 500m x 200m x 100 m in *x*, *y* and *z* directions. The diameter of the spilled pool is related with the diameter of the hole on a reservoir shield. Possible or applied scenarios are refined in a previous section, the resulting pool sizes are further obtained from functions explained in (ioMosaic, 2007). The following scenarios assume the hole at the bottom of the reservoir and 9000 kg of LNG content is spilled in 10 seconds and forms the pool. The evaporation rate of natural gas from the pool surface is calculated by the lumped model and assumed to be 0.16 kg/m2s for the CFD model in a stable sea and stable weather conditions.


Because of the very similar pool sizes scenarios 2 and 4 are calculated once as are scenarios 3 and 5.

#### **4.1 Pool size related to rupture size**

The course of events after an accident could follow different directions before the course of the accident is definitive. Depending on the accident dynamics, and the time and place of occurrence it could lead to different levels of risk. An example is the immediate ignition of a vapour cloud. If this occurs in a populated area, it represents a high risk event, but if it occurs outside the populated area it would be a low risk event. The reason is in the consequences of the event. The second event mentioned would be high risk if the ignition of the cloud would be delayed and the cloud would be transported downwind into a populated area.

The first scenario assumes the evaporation of natural gas from the LNG pool on a water surface of an area of 400 m2. The evaporating gas is moved downwind at 2 m/s; the surrounding temperature is 20oC. The model does not include complex geometry, buildings, and relief that would slightly change the flow dynamics of the vapour cloud.


Table 2. Mandatory input data for Scenario 1:

Fluid Dynamic Models Application in Risk Assessment 75

(Fells et al., 1969, ioMosaic, 2007 & Hissong, 2007) does not devote special attention to this

100 200 300 400 500

CFD 150s

Lumped model

Time (s)

CFD 100s

Le n gh t (m )

0 80 160 100 200 240 325 410

Fig. 3. Methane concentration of 5% computed with a lumped model and CFD model for

The assessment of the installation presented is focused on the accidental events of the LNGC approaching territorial waters though the northern Adriatic separation zone, entering the port and the unloading operation in the sense of individual risk and social risk. Therefore the appropriate approach is to discover and analyse processes and installation elements that could lead to undesired consequences. There are several widely approved methods, like HAZID (Hazard Identification) and HAZOP (Hazard and Operability), described in Macdonald (2004) and in the Guidelines for the management of change for process safety (2008) that are used for structured and systematic examination of a planned or existing process or operation in order to identify and evaluate problems that may represent risks to personnel or equipment. The HAZOP study and report are not presented in the paper, but the evaluated accident events

The population is exposed to risk depending on the reliability of the installation and its components. This is done according to recommendation of the Manual Risk Calculations (2008), Guidelines for the management and change for process safety (2008), DNV reports (2008) and ioMosaic reports (2007). In general the following scenarios are modelled for

The failure frequencies for the scenarios are derived from the Dutch guidelines (2008). The study leads to several accident scenarios with elevated risk and takes into account accidents

• a large leak resulting in a continuous release with a fixed duration (10 minutes);

phenomenon because it is would require overly complex stochastic machinations.

CFD 400s CFD 300s

0

**5. Risk assessment of an LNG terminal** 

frequencies are used for the evaluation of individual and social risk.

• catastrophic rupture with instantaneous release of full inventory;

10

He igh t (m )

**5.1 Scenario based models** 

• a small leak (10 mm)

Scenario 1

vessels:

20

30

The table above explains how the cloud formation is modelled with FDS. There is no phase transition modelled from liquid to gas. Instead of this the gas release with a specified mass flux dynamics or evaporation rate (Table 1) is applied.

The first scenario treats a minor leak and consequently a smaller spilling pool. The scenario is primary used for model testing and for the selection of adequate boundary and initial conditions. It is also used for the comparison of results with further models where the differences of consequences could be compared for different releases.

#### **4.2 Results**

The comparison of results from the lumped model and CFD model explain the adequacy of the analysis. Fig. 2 shows the results of a Methane concentration field. The bold polyline represent the concentration 3000 ppm and the shape of the cloud downwind. The overlaid curves (continuous and dashed) represent the cloud height, calculated by the lumped model described above. Instead of the concentration field the dotted red line represents the cloud height obtained from averaging local heights of the Methane cloud.

Fig. 2. Methane cloud height 8 min after release, concentration CH4 [ppm]

The following in Fig. 3 is the comparison of evaporated gas concentration between two models. 5% concentration is selected because it represents the LFL (lower flammable limit) for methane. The time scale, presented on a graph requires a second axis only for the lumped model results. However, the comparison of CFD fields and the cloud heights of the lumped model show close results and time dynamics.

The obtained results are a good guideline for the setup of further model scenarios. One finding is that the FDS program does not allow for a satisfying definition of the pool evaporation process; therefore the dispersion is overestimated at a later time. The pool evaporation process does not have equal dynamics in both models. The solution found is the combination, where the pool evaporation process is modelled by a lumped model and the cloud dispersion by CFD. This technique is also applied in further scenarios.

A small release of LNG on a water surface reaches the thermal equilibrium fairly rapidly and the evaporated gas quickly dilutes to concentrations below the LFL (low flammable limit) that is 5% in the air. The area at risk of fire is about 150 to 200 m downwind.

The model does not include complex geometry, buildings, or relief that would slightly change the flow dynamics of the vapour cloud. In real situations local gas accumulations could occur in bounded spaces where the gas could exceed 5% concentrations. The literature

The table above explains how the cloud formation is modelled with FDS. There is no phase transition modelled from liquid to gas. Instead of this the gas release with a specified mass

The first scenario treats a minor leak and consequently a smaller spilling pool. The scenario is primary used for model testing and for the selection of adequate boundary and initial conditions. It is also used for the comparison of results with further models where the

The comparison of results from the lumped model and CFD model explain the adequacy of the analysis. Fig. 2 shows the results of a Methane concentration field. The bold polyline represent the concentration 3000 ppm and the shape of the cloud downwind. The overlaid curves (continuous and dashed) represent the cloud height, calculated by the lumped model described above. Instead of the concentration field the dotted red line represents the cloud

1,5E3ppm

0 100 200 300 400 500 600

Lenght (m)

The following in Fig. 3 is the comparison of evaporated gas concentration between two models. 5% concentration is selected because it represents the LFL (lower flammable limit) for methane. The time scale, presented on a graph requires a second axis only for the lumped model results. However, the comparison of CFD fields and the cloud heights of the

The obtained results are a good guideline for the setup of further model scenarios. One finding is that the FDS program does not allow for a satisfying definition of the pool evaporation process; therefore the dispersion is overestimated at a later time. The pool evaporation process does not have equal dynamics in both models. The solution found is the combination, where the pool evaporation process is modelled by a lumped model and the

A small release of LNG on a water surface reaches the thermal equilibrium fairly rapidly and the evaporated gas quickly dilutes to concentrations below the LFL (low flammable

The model does not include complex geometry, buildings, or relief that would slightly change the flow dynamics of the vapour cloud. In real situations local gas accumulations could occur in bounded spaces where the gas could exceed 5% concentrations. The literature

flux dynamics or evaporation rate (Table 1) is applied.

**4.2 Results** 

differences of consequences could be compared for different releases.

height obtained from averaging local heights of the Methane cloud.

1,5E3ppm

3,0E3ppm

Fig. 2. Methane cloud height 8 min after release, concentration CH4 [ppm]

cloud dispersion by CFD. This technique is also applied in further scenarios.

limit) that is 5% in the air. The area at risk of fire is about 150 to 200 m downwind.

lumped model show close results and time dynamics.

5,0E3ppm 1,0E4ppm

0

20

40

Height (m)

60

80

100

(Fells et al., 1969, ioMosaic, 2007 & Hissong, 2007) does not devote special attention to this phenomenon because it is would require overly complex stochastic machinations.

Fig. 3. Methane concentration of 5% computed with a lumped model and CFD model for Scenario 1

### **5. Risk assessment of an LNG terminal**

The assessment of the installation presented is focused on the accidental events of the LNGC approaching territorial waters though the northern Adriatic separation zone, entering the port and the unloading operation in the sense of individual risk and social risk. Therefore the appropriate approach is to discover and analyse processes and installation elements that could lead to undesired consequences. There are several widely approved methods, like HAZID (Hazard Identification) and HAZOP (Hazard and Operability), described in Macdonald (2004) and in the Guidelines for the management of change for process safety (2008) that are used for structured and systematic examination of a planned or existing process or operation in order to identify and evaluate problems that may represent risks to personnel or equipment. The HAZOP study and report are not presented in the paper, but the evaluated accident events frequencies are used for the evaluation of individual and social risk.

#### **5.1 Scenario based models**

The population is exposed to risk depending on the reliability of the installation and its components. This is done according to recommendation of the Manual Risk Calculations (2008), Guidelines for the management and change for process safety (2008), DNV reports (2008) and ioMosaic reports (2007). In general the following scenarios are modelled for vessels:


The failure frequencies for the scenarios are derived from the Dutch guidelines (2008). The study leads to several accident scenarios with elevated risk and takes into account accidents

Fluid Dynamic Models Application in Risk Assessment 77

The societal risk is further presented in Fig. 5 with an F-N curve. The maximum expected number of fatalities is 20 with the probability of 4\*10-7 that is found at the end of F-N curve.

**Societal risk ranking per scenarios** 

The entire societal risk curve (F-N) is located within the ALARP (As Low as Reasonably Practicable) area delimited by the Max and Min Risk Criteria straight lines. This could be taken as an argument that the LNG terminal operations are not harmful to neighbouring populations. Employed personnel on the jetty and pier 2 are the only people who contribute to societal risk, excluding them, the F-N curve would be lower than the Min Risk Criteria. The summary of risk assessment shows that the 1E-6/ avg year individual risk contours do not reach any vulnerable location outside the port area and therefore the level of societal risk is not unacceptable, since the presence of people in the direct vicinity of the jetty is limited and the distance to populated areas like the town centre is too far to pose a societal risk. The most important factor in risk assessment is determining the level of assessment, which then

**Number of fatalities**

P1 P2 P3 P8 P9 Cumulative risk

Societal risk/Average Year (F-N Curve)

Max Risk criteria Min Risk Criteria

0 5 10 15 20 25 30 35 40 45 50 55

Fig. 4. Individual risk for LNG terminal

1,E-09

Fig. 5. Social risk of LNG terminal

1,E-08

1,E-07

**F r e q u e n c y ( / a v g e y e a r )**

1,E-06

1,E-05

1,E-04

1,E-03

that occur due to technical failures. The catastrophic rupture with full lost of inventory is left out of possible technical failure scenarios. This is an important assumption because technical failures could be identified and evaluated by the above methods. The conducted QRA (Quantitative Risk Assessment) takes into account scenarios listed in Table 3.


Table 3. Scenarios modelled in QRA (DNV, 2008)

The worst consequences would result from the largest releases. The rupture of a loading arm during the unloading operation could lead to large spills and influence the risk level. The difference between the first and second scenario is the installation of the emergency shutdown valve on the transfer line. This valve restricts the release time for ruptures of the arms. A value of 120 seconds as maximum release duration is chosen as it is widely accepted to be a conservative value for pipe isolation. Also the scenario that the ESD valves can fail has been taken into account (probability of 0.001 is taken into account). In that case a maximum release time of 30 min is applied.

The defined scenarios are assessed for their impact on external safety by using criteria for individual risk and for societal risk. Fig. 4 shows curves of the individual risk for scenarios P1 to P8 from Table 3, which includes unloading operations on moored tankers. The criterion for individual risk is 1E-6/ avg year for population, a commonly accepted risk throughout the world (Trbojević, 2005). The area with this risk is contained within the installation boundary. The personnel of the installation and of the terminal are the only ones exposed to risk higher than 1E-6/avge year. However, the employed personnel are adequately educated and trained for possible threats and have defined procedures to overcome possible emergencies.

Scenarios with a major influence on risk are larger spills (scenarios P1 and P2) like rupture of a loading arm at lower wind speeds. The reason is that the vapour cloud remains compact for a longer time and moves downwind, possibly into the direction of a populated area.

The populated neighbourhoods located south and north of the port are exposed to a much lower risk than 1E-6/avge year. In this case the LNG unloading operation, considering the eight most likely accident events, is safe for people living or passing through the LNG terminal neighbourhood. Personnel working on the LNG terminal on pier 2 (Fig. 4) and on neighbouring terminals on pier 1 are covered by the risk zone 1E-6/ avg year.

that occur due to technical failures. The catastrophic rupture with full lost of inventory is left out of possible technical failure scenarios. This is an important assumption because technical failures could be identified and evaluated by the above methods. The conducted QRA

P3- Loading arm – 10% leak (max. 50mm) 1.6" 16 kg/s 7.5E-04

rupture with ESD 16" 9 kg/s 2.5E-08

(max. 50mm) 1.6" 0.1 kg/s 2.5E-04 P6- Jetty drum -catastrophic rupture 30 m3 4.8E-07 P7- Jetty drum - release in 10 min 30 m3 4.8E-07 P8- Jetty drum -10 mm leak 10 mm 30m3 9.6E-06

The worst consequences would result from the largest releases. The rupture of a loading arm during the unloading operation could lead to large spills and influence the risk level. The difference between the first and second scenario is the installation of the emergency shutdown valve on the transfer line. This valve restricts the release time for ruptures of the arms. A value of 120 seconds as maximum release duration is chosen as it is widely accepted to be a conservative value for pipe isolation. Also the scenario that the ESD valves can fail has been taken into account (probability of 0.001 is taken into account). In that case a

The defined scenarios are assessed for their impact on external safety by using criteria for individual risk and for societal risk. Fig. 4 shows curves of the individual risk for scenarios P1 to P8 from Table 3, which includes unloading operations on moored tankers. The criterion for individual risk is 1E-6/ avg year for population, a commonly accepted risk throughout the world (Trbojević, 2005). The area with this risk is contained within the installation boundary. The personnel of the installation and of the terminal are the only ones exposed to risk higher than 1E-6/avge year. However, the employed personnel are adequately educated and trained for possible threats and have defined procedures to

Scenarios with a major influence on risk are larger spills (scenarios P1 and P2) like rupture of a loading arm at lower wind speeds. The reason is that the vapour cloud remains compact for a longer time and moves downwind, possibly into the direction of a populated area. The populated neighbourhoods located south and north of the port are exposed to a much lower risk than 1E-6/avge year. In this case the LNG unloading operation, considering the eight most likely accident events, is safe for people living or passing through the LNG terminal neighbourhood. Personnel working on the LNG terminal on pier 2 (Fig. 4) and on

neighbouring terminals on pier 1 are covered by the risk zone 1E-6/ avg year.

Volume (m³) or release rate (kg/s) - Release time



Failure frequency (per year)

(Quantitative Risk Assessment) takes into account scenarios listed in Table 3.

Scenario Diameter

(Emergency Shut Down) 16" 700 kg/s

P2- Loading arm – rupture without ESD 16" 705 kg/s

P1- Loading arm – rupture with ESD

P5- Vapour return arm – 10% leak

Table 3. Scenarios modelled in QRA (DNV, 2008)

maximum release time of 30 min is applied.

overcome possible emergencies.

P4- Vapour return arm –

Fig. 4. Individual risk for LNG terminal

The societal risk is further presented in Fig. 5 with an F-N curve. The maximum expected number of fatalities is 20 with the probability of 4\*10-7 that is found at the end of F-N curve.

Fig. 5. Social risk of LNG terminal

The entire societal risk curve (F-N) is located within the ALARP (As Low as Reasonably Practicable) area delimited by the Max and Min Risk Criteria straight lines. This could be taken as an argument that the LNG terminal operations are not harmful to neighbouring populations. Employed personnel on the jetty and pier 2 are the only people who contribute to societal risk, excluding them, the F-N curve would be lower than the Min Risk Criteria. The summary of risk assessment shows that the 1E-6/ avg year individual risk contours do not reach any vulnerable location outside the port area and therefore the level of societal risk is not unacceptable, since the presence of people in the direct vicinity of the jetty is limited and the distance to populated areas like the town centre is too far to pose a societal risk. The most important factor in risk assessment is determining the level of assessment, which then

Fluid Dynamic Models Application in Risk Assessment 79

variables, soot density and temperature, are presented in minutes time steps trough the entire tunnel length. Comparing the obtained data in a table allows the analyses of the ventilation conditions for different heat releases from fires. The second step is to add additional criteria of human behaviour inside the tunnel (evacuation) and human resistance

In order to identify the interactive and uniting relationships in a system, analysis is necessary to replace the apparent structure of individual statements on the components of a system and their relationships with their underlying common logical structure (system analysis). For example, if we are dealing with a system which we call "a chemical process plant", we get at its various components successively, by means of deductive analysis: the buildings, the operators, the storage tanks, the control systems, the operating procedures, etc.. Each such component is thrown into the modelling reality by a distinct act of noticing, and is steadily held together with those components already segregated. The aim of system analysis is to investigate the system's behaviour (i.e. the succession of its states over time) on the basis of its components' changes with time. The results of system analysis can be expressed in qualitative and quantitative terms (statements resulting from "qualitative

The deterministic approach breathe into the analysis of the greater part of physical events like fire source characteristic and its dynamics, the operation of the ventilation system and other conditions as well as their reciprocal interactions. The approach leads also to the definition of the technical system "safety efficiency" in the range of possibilities that exist in a real word and are functionally descriptive. When the approach is used in practice, we should define a number of "safety categories" base on events probability and consequences

> A Probability once in a year 0.3 – 3 B Possible but not likely 0.03 – 0.3 C Unlikely 0.003 – 0.03 D Very unlikely 0.0003 – 0.003 E Remote 0.00003 – 0.0003

1 Catastrophic Multiple fatalities 2 Major Single fatality, multiple injuries 3 Very serious Permanently disable injuries 4 Serious Serious injury, full recovery 5 Minor Lost time injury, short absence

Table 4. Deterministic safety analysis – supposed safety categories (Kirchsteiger, 1999)

Qualitative definition Underlying quantitative

Qualitative definition Underlying semi-quantitative

definition (times per year)

definition

from work

to the elevated gas concentrations and temperature (Haack, 1998 & 2002).

analysis" and numbers resulting from "quantitative analysis").

for the individual risk. The example in presented in a Table 4.

(a) Likelihood categories Severity category (frequencies)

> (b) Consequence categories Severity category (consequences)

**7.1 Methodological approach on tunnel safety** 

determines which scenarios and types of accident are investigated. Here we dealt only with technical failures during the unloading operation with major predictable consequences.

Modern industrial installations are no longer problematic regarding safety aspects mainly because of the application of state of the art standards for construction, materials, and operations and so on. During last decade a lot of attention has been focused on unpredictable risks that are human and societal related; that is to say, terrorism, which has not yet been well defined. A lot of work is done by developing countries to overcome this risk and to include it in different standards and procedures. But up to now no one standard regarding construction or materials of industrial installations has been changed because of the risk posed by terrorism, though several new operational procedures have been applied in practice.

#### **6. Conclusion on LNG risk assessment**

The investigation conducted is focused on the discovery of an appropriate approach to risk assessment for an LNG delivery terminal located close to a populated area. The assessment should focus mainly on the identification of accident scenarios which results in individual risk higher than 1E-3 / avg year for the neighbouring population. Because the risk calculations are most influenced by large accident consequences, these accidents should be analysed in detail to avoid over or under estimation of fast computing consequence (dispersion) models, which are usually used in risk assessment software. The comparison of a lumped model and CFD model is therefore conducted and the differences are analysed and discussed. After obtaining satisfying results the individual and societal risk is computed for the specified number of accident scenarios. The lumped model approach for consequence calculation is the best choice in the initial phase of risk analyses when several simulations need to be done. Additional detailed analyses are required when the risk is near the limit of that ALARP region. The numerical simulation of a problematic scenario could explain whether the consequences and therefore the risk is really too high, or whether it should fall within an acceptable area. Usually this is not a satisfying solution, so regarding this some modifications of the project are commonly proposed.

A possibility for further development is on a risk model based on CFD that relies on several overlaying CFD scenarios in combination with risk criteria functions and obtained individual risk levels. Depending on the accuracy of analysed scenarios the risk range could be much more realistic than in lumped models where the range is intentionally more conservative. In any case, the method would be very time consuming yet could be quite useful in the last stage of risk assessment.

#### **7. Fluid dynamic models for road tunnels risk assessment**

The definition of the deterministic approach in safety analyses arises from the need to understand the conditions that emerge during a fire accident in a road tunnel. The key factor of the tunnel operations during the fire is the ventilation, which during the initial phases of the fire have a strong impact on the evacuation of people and later on the access of the intervention units in the tunnel. The text presents the use of the CFD model in the tunnel safety assessment process. The set-up of the initial and boundary conditions and the requirement for grid density found from validation tests of an FDS (Fire Dynamics Simulator) is used to prepare different kinds of fire scenarios in different ventilation conditions; natural, semi transverse, transverse and longitudinal ventilation. The observed

determines which scenarios and types of accident are investigated. Here we dealt only with

Modern industrial installations are no longer problematic regarding safety aspects mainly because of the application of state of the art standards for construction, materials, and operations and so on. During last decade a lot of attention has been focused on unpredictable risks that are human and societal related; that is to say, terrorism, which has not yet been well defined. A lot of work is done by developing countries to overcome this risk and to include it in different standards and procedures. But up to now no one standard regarding construction or materials of industrial installations has been changed because of the risk posed by terrorism,

The investigation conducted is focused on the discovery of an appropriate approach to risk assessment for an LNG delivery terminal located close to a populated area. The assessment should focus mainly on the identification of accident scenarios which results in individual risk higher than 1E-3 / avg year for the neighbouring population. Because the risk calculations are most influenced by large accident consequences, these accidents should be analysed in detail to avoid over or under estimation of fast computing consequence (dispersion) models, which are usually used in risk assessment software. The comparison of a lumped model and CFD model is therefore conducted and the differences are analysed and discussed. After obtaining satisfying results the individual and societal risk is computed for the specified number of accident scenarios. The lumped model approach for consequence calculation is the best choice in the initial phase of risk analyses when several simulations need to be done. Additional detailed analyses are required when the risk is near the limit of that ALARP region. The numerical simulation of a problematic scenario could explain whether the consequences and therefore the risk is really too high, or whether it should fall within an acceptable area. Usually this is not a satisfying solution, so regarding

A possibility for further development is on a risk model based on CFD that relies on several overlaying CFD scenarios in combination with risk criteria functions and obtained individual risk levels. Depending on the accuracy of analysed scenarios the risk range could be much more realistic than in lumped models where the range is intentionally more conservative. In any case, the method would be very time consuming yet could be quite

The definition of the deterministic approach in safety analyses arises from the need to understand the conditions that emerge during a fire accident in a road tunnel. The key factor of the tunnel operations during the fire is the ventilation, which during the initial phases of the fire have a strong impact on the evacuation of people and later on the access of the intervention units in the tunnel. The text presents the use of the CFD model in the tunnel safety assessment process. The set-up of the initial and boundary conditions and the requirement for grid density found from validation tests of an FDS (Fire Dynamics Simulator) is used to prepare different kinds of fire scenarios in different ventilation conditions; natural, semi transverse, transverse and longitudinal ventilation. The observed

technical failures during the unloading operation with major predictable consequences.

though several new operational procedures have been applied in practice.

this some modifications of the project are commonly proposed.

**7. Fluid dynamic models for road tunnels risk assessment** 

useful in the last stage of risk assessment.

**6. Conclusion on LNG risk assessment** 

variables, soot density and temperature, are presented in minutes time steps trough the entire tunnel length. Comparing the obtained data in a table allows the analyses of the ventilation conditions for different heat releases from fires. The second step is to add additional criteria of human behaviour inside the tunnel (evacuation) and human resistance to the elevated gas concentrations and temperature (Haack, 1998 & 2002).

#### **7.1 Methodological approach on tunnel safety**

In order to identify the interactive and uniting relationships in a system, analysis is necessary to replace the apparent structure of individual statements on the components of a system and their relationships with their underlying common logical structure (system analysis). For example, if we are dealing with a system which we call "a chemical process plant", we get at its various components successively, by means of deductive analysis: the buildings, the operators, the storage tanks, the control systems, the operating procedures, etc.. Each such component is thrown into the modelling reality by a distinct act of noticing, and is steadily held together with those components already segregated. The aim of system analysis is to investigate the system's behaviour (i.e. the succession of its states over time) on the basis of its components' changes with time. The results of system analysis can be expressed in qualitative and quantitative terms (statements resulting from "qualitative analysis" and numbers resulting from "quantitative analysis").

The deterministic approach breathe into the analysis of the greater part of physical events like fire source characteristic and its dynamics, the operation of the ventilation system and other conditions as well as their reciprocal interactions. The approach leads also to the definition of the technical system "safety efficiency" in the range of possibilities that exist in a real word and are functionally descriptive. When the approach is used in practice, we should define a number of "safety categories" base on events probability and consequences for the individual risk. The example in presented in a Table 4.



Fluid Dynamic Models Application in Risk Assessment 81

D is a molecular diffusivity. *p* is the perturbation pressure caused by pressure differences,

The effect of the flow field turbulence is modelled using LES (Large Eddy Simulation), in which the large scale eddies are computed directly and the sub-grid scale dissipative

Further the combustion model is based on the assumption that the combustion in mixingcontrolled. This implies that all species of interest can be described in terms of the mixture fraction Z. Heat from the reaction of fuel and oxygen is released along an infinitely thin sheet where Z takes on its stoichiometric value as determined by the solution of the transport equation for Z. The state relations are calculated for a stoichiometric reaction of C7H16 (Oil), which is proposed by (McGrattan, 2001, Heskestad, 1995 & Mingchung, 1999)

The idea is based on the creation of a deterministic risk matrix as it is showed in the Table 5. The safety category is represented by the power of the fire and the type of ventilation at different strengths and on other side the consequences are evaluated in the time during the progress of the fire. The risk criteria are defined as a relation between the hot smoke layer height, the distance from the fire position and the evacuation time of the users. In case the speed of the smoke is higher than the speed of the evacuation and in case the height of the

All together 12 tunnel fire scenarios are presented. Three levels of fire force are simulated, each with four different types of ventilation. The span of the fire force is between 20 MW, 50 MW and 100 MW whilst the ventilation is sorted from the less to the more effective: 1 – natural, 2 – longitudinal, 3 – semi transverse and 4 – transverse or improved transverse ventilation. The whole section of the simulated tunnel is 650 m long, the other dimensions are width 10 m and height 8 m or 6 m when the roof is lowered. Though the dimensions and shape of the tunnel tube partly differentiate among them that does not influence what happens during the fire. That is why ordinary skeleton measurements are chosen. The geometry of the

tunnel model, the type of ventilation and the location of the fire are shown on Fig. 6.

The fire is placed on a distance of 300 m in all the models, it differs only in the size of the burning area. The focus point is defined as the heat source to which the combusted model calculates the mass transfer on the base of the accorded combusting reaction and the oxygen consumption. The ignition point is shown in Fig. 6. When we define the base igniting temperature, heat conductivity, calorific value, etc. (depending on the models demands) it is treated in the model as combustible substance and it cooperates with the generation of heat

processes are modelled (Sagaut, 2002). The unknown sub-grid stress tensor

hot layer is higher than the speed of the evacuation, the risk is high.

 the viscosity stress tensor and *k* the thermal conductivity. *<sup>c</sup> q*′′′ and ∇**q**r are the source terms of chemical reaction and radiation, respectively. The radiation term has a negative sign

is a density, **u** is a velocity vector, Z is the mixture fraction, T the temperature and

<sup>∂</sup> + ⋅∇ = −∇⋅ +∇⋅ ∇ ′′′ <sup>∂</sup> (11)

τ

is modelled

*<sup>p</sup>* **u q** *c R <sup>T</sup> c T <sup>q</sup> k T*

*t*

ρ

because it represents a heat sink.

and called a Crude oil reaction.

**8. Tunnel fire analysis** 

**8.1 Tunnel fire scenario** 

by Smagorinsky model (Lesieur, 1997).

Where ρ

τ

Note that, in these schemes, a quantitative definition is often given in addition to the qualitative definition, mainly to ensure consistency in the course of the analysis and provide benchmarks ("semi-quantitative analysis"). In schemes of this type, the assessment team, usually comprising members of line management, safety engineers and operations personnel, will first identify all hazards, using HAZOP or similar approaches, and then assigns a severity category to each of these, for both likelihood and consequences (PIARC, 2003 & Brussaard, 2004).

Following the assumptions of (Kirchsteiger, 1999), a "risk matrix" would then be defined as a 5 x 5 matrix with each side corresponding to one severity category.


Table 5. Deterministic safety analysis – example of risk matrix (Kirchsteiger, 1999)

Different shading in a table indicates different risk levels. Hazards with high assessments, such as A1, B1 and A2 in the black squares, are thought of as being very severe and requiring immediate action to reduce. Hazards with low assessments, such as E5, E4 and D5 in the white squares, are considered to require no further action. Hazards between these two (grey squares) are considered worthy of some improvement if a cost-effective solution can be found.

#### **7.2 Computer models and simulations**

Deterministic models that would consider all physical parameters are almost unfeasible in practice and if feasible would require very complex and time consuming computations. The application of deterministic analyses results in practice is conditioned by the simplification of some physical phenomena (like turbulence) (Gasser et al., 2002 & Floyd et al., 2001).

The fluid flow is modelled by solving the basic conservation equations. Those are conservation of mass (8), conservation of mixture fraction (9), conservation of momentum (10) and conservation of energy (11) using a form for low Mach number (McGrattan, 2001 & Fletcher, 1991). The approximation involves the filtering out of acoustic waves.

$$\frac{\partial \rho}{\partial t} + \nabla \cdot \rho \mathbf{u} = 0 \tag{8}$$

$$\frac{\partial \rho}{\partial t}(\rho Z) + \nabla \cdot \rho Z \mathbf{u} = \nabla \cdot \rho \text{DVD} \tag{9}$$

$$\rho \left( \frac{\partial \mathbf{u}}{\partial t} + \frac{1}{2} \nabla \left| \mathbf{u} \right|^2 - \mathbf{u} \times \boldsymbol{\omega} \right) + \nabla \tilde{p} = (\rho - \rho\_{\text{os}}) \mathbf{g} + \nabla \cdot \boldsymbol{\tau} \tag{10}$$

Note that, in these schemes, a quantitative definition is often given in addition to the qualitative definition, mainly to ensure consistency in the course of the analysis and provide benchmarks ("semi-quantitative analysis"). In schemes of this type, the assessment team, usually comprising members of line management, safety engineers and operations personnel, will first identify all hazards, using HAZOP or similar approaches, and then assigns a severity category to each of these, for both likelihood and consequences (PIARC,

Following the assumptions of (Kirchsteiger, 1999), a "risk matrix" would then be defined as

Table 5. Deterministic safety analysis – example of risk matrix (Kirchsteiger, 1999)

Different shading in a table indicates different risk levels. Hazards with high assessments, such as A1, B1 and A2 in the black squares, are thought of as being very severe and requiring immediate action to reduce. Hazards with low assessments, such as E5, E4 and D5 in the white squares, are considered to require no further action. Hazards between these two (grey squares) are considered worthy of some improvement if a cost-effective solution can

Deterministic models that would consider all physical parameters are almost unfeasible in practice and if feasible would require very complex and time consuming computations. The application of deterministic analyses results in practice is conditioned by the simplification of some physical phenomena (like turbulence) (Gasser et al., 2002 & Floyd et al., 2001). The fluid flow is modelled by solving the basic conservation equations. Those are conservation of mass (8), conservation of mixture fraction (9), conservation of momentum (10) and conservation of energy (11) using a form for low Mach number (McGrattan, 2001 &

**u** 0

 ρ

> ρ ρ

<sup>∂</sup> (10)

∂ (8)

∂ (9)

+

∞

τ

ρ

( ) *Z Z DZ* **u**

**u u g** <sup>1</sup> <sup>2</sup> ( ) <sup>2</sup> *<sup>p</sup> <sup>t</sup>*

<sup>∂</sup> +∇⋅ =

<sup>∂</sup> +∇⋅ =∇⋅ ∇

ω

<sup>∂</sup> + ∇ − × +∇ = − ∇⋅

Fletcher, 1991). The approximation involves the filtering out of acoustic waves.

*t* ρ

ρρ

*t* ρ

**u**

ρ

Consequences ''Severity category'' **5 4 3 2 1** 

a 5 x 5 matrix with each side corresponding to one severity category.

2003 & Brussaard, 2004).

''Likelihood'' ''Severity category ''

> **A B C D E**

be found.

**7.2 Computer models and simulations** 

$$
\rho c\_p \left( \frac{\partial T}{\partial t} + \mathbf{u} \cdot \nabla T \right) = \dot{q}\_c''' - \nabla \cdot \mathbf{q}\_R + \nabla \cdot k \nabla T \tag{11}
$$

Where ρ is a density, **u** is a velocity vector, Z is the mixture fraction, T the temperature and D is a molecular diffusivity. *p* is the perturbation pressure caused by pressure differences, τ the viscosity stress tensor and *k* the thermal conductivity. *<sup>c</sup> q*′′′ and ∇**q**r are the source terms of chemical reaction and radiation, respectively. The radiation term has a negative sign because it represents a heat sink.

The effect of the flow field turbulence is modelled using LES (Large Eddy Simulation), in which the large scale eddies are computed directly and the sub-grid scale dissipative processes are modelled (Sagaut, 2002). The unknown sub-grid stress tensor τ is modelled by Smagorinsky model (Lesieur, 1997).

Further the combustion model is based on the assumption that the combustion in mixingcontrolled. This implies that all species of interest can be described in terms of the mixture fraction Z. Heat from the reaction of fuel and oxygen is released along an infinitely thin sheet where Z takes on its stoichiometric value as determined by the solution of the transport equation for Z. The state relations are calculated for a stoichiometric reaction of C7H16 (Oil), which is proposed by (McGrattan, 2001, Heskestad, 1995 & Mingchung, 1999) and called a Crude oil reaction.

#### **8. Tunnel fire analysis**

The idea is based on the creation of a deterministic risk matrix as it is showed in the Table 5. The safety category is represented by the power of the fire and the type of ventilation at different strengths and on other side the consequences are evaluated in the time during the progress of the fire. The risk criteria are defined as a relation between the hot smoke layer height, the distance from the fire position and the evacuation time of the users. In case the speed of the smoke is higher than the speed of the evacuation and in case the height of the hot layer is higher than the speed of the evacuation, the risk is high.

#### **8.1 Tunnel fire scenario**

All together 12 tunnel fire scenarios are presented. Three levels of fire force are simulated, each with four different types of ventilation. The span of the fire force is between 20 MW, 50 MW and 100 MW whilst the ventilation is sorted from the less to the more effective: 1 – natural, 2 – longitudinal, 3 – semi transverse and 4 – transverse or improved transverse ventilation.

The whole section of the simulated tunnel is 650 m long, the other dimensions are width 10 m and height 8 m or 6 m when the roof is lowered. Though the dimensions and shape of the tunnel tube partly differentiate among them that does not influence what happens during the fire. That is why ordinary skeleton measurements are chosen. The geometry of the tunnel model, the type of ventilation and the location of the fire are shown on Fig. 6.

The fire is placed on a distance of 300 m in all the models, it differs only in the size of the burning area. The focus point is defined as the heat source to which the combusted model calculates the mass transfer on the base of the accorded combusting reaction and the oxygen consumption. The ignition point is shown in Fig. 6. When we define the base igniting temperature, heat conductivity, calorific value, etc. (depending on the models demands) it is treated in the model as combustible substance and it cooperates with the generation of heat

Fluid Dynamic Models Application in Risk Assessment 83

The simulation results are presented on levels of fire force and types of tunnel ventilation. The consequences of the distance of the smoke and the temperature are qualitatively evaluated from the current and temperature field. With this, it must be noted that mistakes are possible in calculating the average value in different time and space steps, which are limited with the unified way of average calculating. With that, it is true that the risk of exposure to smoke is that the participant is exposed in the moment when the smoke reaches him. The most risky examples are the ones when the participant does not start with the immediate self-rescue procedure after the start of the fire and the second when the spreading speed of the smoke is higher than the self-rescue procedure speed of the participants in the tunnel. The other risk criterion is high temperature that usually has a lower contribution to the risk than smoke. In most cases this depends on the way of

The limit value of the concentration of smoke particles (PM10 heavy particles with the diameter up to 10 μm) is 1000 mg/m3 (Vidmar et al., 2003) and the limit temperature is 50oC. Though the smoke particles are less problematic from a poisonous point of view, than other combustible products (CO2-carbon dioxide, CO-carbon monoxide, HCNhydrogen cyanide, HCl-hydrogen chloride, etc.) their relation to the concentration is conditional and often very similar. From different experiments in the Memorial Tunnel (1996) it can be found for example concentrations of smoke particles and CO in relation around 10:1. A similar relation can be also found on toxic levels of these products. LC50 (Lethal Concentration 50 %) for soot particles is 30g/m3 in a 30 min exposure or 1-3g min/m3 LC50, for CO is 2000-3500 ppm, which is 2300-4000 mg/m3 in a 30-60 min exposure. The limit temperature values of human endurance are according to Gann (1994) 1000C for 30 min and 750C for 60 min of exposure. Because this information is true for an adult man it is the most optimal. But within the same research there are difficulties in breathing already at 65oC of air temperature. Taking this into account there are two values that are used in the result analysis. The chosen limit concentration of smoke particles is

The risk or consequences are divided in five categories that are shown in the Table 6, these are:

In the result analysis each category matches a logical inscription and it conditions with the time from the start of the simulation, any distance from the fire area, fire force, way of ventilation, limit value of the concentration of smoke particles and the limit temperature.

**8.3 Parameters and approach to the result analysis** 

1000 m g/m3 and the limit temperature is 50oC.

1. LR – low risk: smaller injury

3. SR – serious risk: permanent injury

Then follow the conditional clauses of each category:

ASD - Average smoke density value in profile [mg/m3]

MR: ASDL > 500. ∧ SLH > ASLH

VHR: ASDL > 500. ∧ SLH < ASLH

Where the abbreviations mean:

EHR: ((TH ∨ VHT) ∧ AT > 50.) ∨ ATL > 50.

LR: ASD < 500.

SR: ASD> 500.

5. EHR – extremely high risk: numerous casualties

2. MR – medium risk: serious injury with full recovery

4. VHR – very high risk: low casualty number (1-3), numerous injured

ventilation.

in the combusting model. In case of the described scenarios, the base is relatively small or of small volume, that is why the heat contribution of the burning base is only a few percent of the defined freed heat of the boundary condition.

#### **8.2 Initial, boundary conditions and discretization**

The definition of the initial and boundary conditions is a peculiarity of each model. Four elementary ways of ventilation are discussed: natural, longitudinal, semi transverse and transverse. In definition of the geometry of the tunnel, the natural and longitudinal ventilation are discussed together and the semi transverse and transverse ones also in the same way (Woodburn, 1996). The comparison is shown on the lower picture:

d) Transverse ventilation

Fig. 6. Model geometry in different ways of ventilation

It is clear that the tunnel models with natural and longitudinal ventilation take the whole section of the tunnel, however the tunnel models with semi transverse and transverse ventilation consider only the light section of the tunnel (without the ventilation drains). The suction flaps are defined with the speed margin condition on the limit of the calculating domain (Cheng et al., 2001).

in the combusting model. In case of the described scenarios, the base is relatively small or of small volume, that is why the heat contribution of the burning base is only a few percent of

The definition of the initial and boundary conditions is a peculiarity of each model. Four elementary ways of ventilation are discussed: natural, longitudinal, semi transverse and transverse. In definition of the geometry of the tunnel, the natural and longitudinal ventilation are discussed together and the semi transverse and transverse ones also in the

a) b)

c) d)

It is clear that the tunnel models with natural and longitudinal ventilation take the whole section of the tunnel, however the tunnel models with semi transverse and transverse ventilation consider only the light section of the tunnel (without the ventilation drains). The suction flaps are defined with the speed margin condition on the limit of the calculating

a) Natural ventilation b) Longitudinal ventilation c) Semi transverse ventilation

Fig. 6. Model geometry in different ways of ventilation

d) Transverse ventilation

domain (Cheng et al., 2001).

same way (Woodburn, 1996). The comparison is shown on the lower picture:

the defined freed heat of the boundary condition.

**8.2 Initial, boundary conditions and discretization** 

#### **8.3 Parameters and approach to the result analysis**

The simulation results are presented on levels of fire force and types of tunnel ventilation. The consequences of the distance of the smoke and the temperature are qualitatively evaluated from the current and temperature field. With this, it must be noted that mistakes are possible in calculating the average value in different time and space steps, which are limited with the unified way of average calculating. With that, it is true that the risk of exposure to smoke is that the participant is exposed in the moment when the smoke reaches him. The most risky examples are the ones when the participant does not start with the immediate self-rescue procedure after the start of the fire and the second when the spreading speed of the smoke is higher than the self-rescue procedure speed of the participants in the tunnel. The other risk criterion is high temperature that usually has a lower contribution to the risk than smoke. In most cases this depends on the way of ventilation.

The limit value of the concentration of smoke particles (PM10 heavy particles with the diameter up to 10 μm) is 1000 mg/m3 (Vidmar et al., 2003) and the limit temperature is 50oC. Though the smoke particles are less problematic from a poisonous point of view, than other combustible products (CO2-carbon dioxide, CO-carbon monoxide, HCNhydrogen cyanide, HCl-hydrogen chloride, etc.) their relation to the concentration is conditional and often very similar. From different experiments in the Memorial Tunnel (1996) it can be found for example concentrations of smoke particles and CO in relation around 10:1. A similar relation can be also found on toxic levels of these products. LC50 (Lethal Concentration 50 %) for soot particles is 30g/m3 in a 30 min exposure or 1-3g min/m3 LC50, for CO is 2000-3500 ppm, which is 2300-4000 mg/m3 in a 30-60 min exposure. The limit temperature values of human endurance are according to Gann (1994) 1000C for 30 min and 750C for 60 min of exposure. Because this information is true for an adult man it is the most optimal. But within the same research there are difficulties in breathing already at 65oC of air temperature. Taking this into account there are two values that are used in the result analysis. The chosen limit concentration of smoke particles is 1000 m g/m3 and the limit temperature is 50oC.

The risk or consequences are divided in five categories that are shown in the Table 6, these are:


In the result analysis each category matches a logical inscription and it conditions with the time from the start of the simulation, any distance from the fire area, fire force, way of ventilation, limit value of the concentration of smoke particles and the limit temperature. Then follow the conditional clauses of each category:

$$\text{LR} \colon \qquad \text{ASD} \le 500.$$

MR: ASDL > 500. ∧ SLH > ASLH

SR: ASD> 500.

VHR: ASDL > 500. ∧ SLH < ASLH

EHR: ((TH ∨ VHT) ∧ AT > 50.) ∨ ATL > 50.

Where the abbreviations mean:

ASD - Average smoke density value in profile [mg/m3]

Fluid Dynamic Models Application in Risk Assessment 85

With the presented model the possibilities of the analysis or the following of the movement of the users in the tunnel increase additionally. The calculated locations are then used for checking the temperature and the smoke concentration on the ground in these places and

The first level of risk is presented by the presence of smoke that includes the first four risk stages, the presence of high temperature contributes additional (the highest) risk stage. The Table 6 presents a deterministic register of risk for a constant location in a tunnel during a fire that is 252 m north of the fire. The picture that we get whit this is very representative because it confirms the theory on safety analyses. From the table the safety categories can be seen and appropriate consequences can be allocated. Table 6 has especially a comparative purpose for finding the influence of different ways of ventilation on the fire dynamic, smoke and temperature development. We can logically assume that the risk in low fire force is lower in comparison with bigger fires. Following the same logic along with the consideration of different ways of ventilation and manner of management it soon becomes difficult. One of the noticeable differences is the level of the calculated risk (MR – medium risk) in longitudinal ventilation of a 20 MW fire. It is expected that the increased risk also appears in the 50 MW fire but it is not so. The search for a cause is difficult because this is hidden in the fluid dynamics during the fire, taking into account that the geometry, the discreetness and the initial and boundary conditions

The second important result in the table is the possibility of analysing the influence of turning on the ventilators on the forming of the smoke curtain. It is especially noticeable in the transverse ventilation of 50 MW and 100MW fires where on the turning on, local increased temperatures and smoke concentrations occur. Table 6 is made as a functionally dependent dynamic register, which chooses the calculated values from the data base with the changing of the observed location and on the basis of conditions from the chapter above,

With the evacuation model is possible to define the tunnel user's location in time intervals of one minute on the bases of the starting user's position, delay with the self-rescue procedure and the walking speed. In this way it is possible to predict the tunnel users movement for the following 15 min and check the smoke concentration and the temperature height to which they will be exposed or establish the risk level. The table represents a conceptual

*SJt* <sup>1</sup> *X* <sup>−</sup> - Position of the user regarding the location of the fire (N or S) in precedent time

*XL*<sup>0</sup> - Starting location of the self-rescue observed from a starting portal [m]

*XLt* - Users location in the observed time period [m]

*SR t* - Delay of the self-rescue after the start of the fire [s]

*t* - Momentarily observed time period [s]

(except the force of the fire) are unaltered.

calculates the risk.

*Xpoz* - Locations of the fire observed form a starting portal [m]

The marks in the note represent:

*vH* - Walking speed [m/s]

consequently the level of risk.

step [m]

**8.5 Results** 

ASDL - Average smoke density value in layer [mg/m3] SLH – Smoke layer height [m] ASLH – Allowed smoke layer height [m] AT – Average temperature in profile [oC] ATL – Average temperature in layer [oC] TLH – Temperature layer height [m]

## ATLH – Allowed temperature layer height [m]

#### **8.4 Evacuation model**

The easier discussing of results is enforced with the understanding of people behaviour during a fire in the tunnel which after the spotted fire begins with the self-rescue procedure. This is a withdrawal from the tunnel or to the first transverse passage in two tube tunnel scenarios. The movement of the people in similar conditions is very unpredictable, some become immediately aware of the danger and begin with the selfrescue procedure others do not perceive the danger in time and start with the self-rescue procedure too late. On self-rescue there is a simplified model of people movement in the tunnel. The model takes into consideration the elementary movement parameters as: start of the self-rescue, walking speed, tunnel length and logical curiosity that in the initial location north or south arranges the movement direction north or south. With this the possibility of a tunnel user approaching the fire during the self-rescue procedure is excluded in the model. With a program the self-rescue procedure is defined with the following conditional note:

$$\begin{aligned} \text{IFF}\left(X\_{S\_{l-1}}=\text{'S'}\land\text{XL}\_{0}-v\_{H}\cdot\left(t-t\_{SR}\right)\leq 0\right) &\text{ then \text{'UUT'}}\\else\\ \text{else} \\ \text{IFF}\left(X\_{S\_{l-1}}=\text{'S'}\right) &\text{then } \text{XL}\_{0}-v\_{H}\cdot\left(t-t\_{SR}\right) &\text{If } \text{F}\left(\text{XL}\_{t}=\text{'OLIT'}\right) &\text{then \text{'UUT'}}\\else\\ \text{else} &\text{IFF}\left(X\_{S\_{l-1}}=\text{'I'}\land\text{XL}\_{0}+v\_{H}\cdot\left(t-t\_{SR}\right)\geq \text{Lp}\right) &\text{then \text{'UUT'}} &\text{else}\\else\\ \text{else} &\\ \text{XL}\_{0}+v\_{H}\cdot\left(t-t\_{SR}\right) &\text{END }\text{IF}\\\ \text{END IF} &\\\ \text{END IF} &\\\ \text{END IF} &\\\ \text{END IF} &\\\ \text{Odd}\downarrow\text{jonot}\ \stackrel{\text{\textquotedblleft}}{\text{\textquotedblright}} &\text{Odd}\downarrow\text{jonot}\ \stackrel{\text{\textquotedblleft}}{\text{\textquotedblright}} &\text{Odd}\downarrow\text{jonot}\ \stackrel{\text{\textquotedblleft}}{\text{\textquotedblleft}} &\text{if}\downarrow\text{jonot}\ \stackrel{\text{\textquotedblleft}}{\text{\textquotedblleft}} \\\ \text{U} &\\\ \text{U}\text{\uquotedblright} &\text{N-S'} &\\\ \text{U}\text{\uquotedblleft} &\text{N-S'} &\\\ \text{U}\text{\uquotedblright} &\text{N-S'} &\\\ \text{U}\text{\uquotedblleft} &\text{N-S'} &\\\ \text{U}\text{\uquotedbl} &\text$$

Fig. 7. Evacuation model

The marks in the note represent:


*XLt* - Users location in the observed time period [m]

*Xpoz* - Locations of the fire observed form a starting portal [m]

*vH* - Walking speed [m/s]

*t* - Momentarily observed time period [s]

*SR t* - Delay of the self-rescue after the start of the fire [s]

With the presented model the possibilities of the analysis or the following of the movement of the users in the tunnel increase additionally. The calculated locations are then used for checking the temperature and the smoke concentration on the ground in these places and consequently the level of risk.

#### **8.5 Results**

84 Fluid Dynamics, Computational Modeling and Applications

The easier discussing of results is enforced with the understanding of people behaviour during a fire in the tunnel which after the spotted fire begins with the self-rescue procedure. This is a withdrawal from the tunnel or to the first transverse passage in two tube tunnel scenarios. The movement of the people in similar conditions is very unpredictable, some become immediately aware of the danger and begin with the selfrescue procedure others do not perceive the danger in time and start with the self-rescue procedure too late. On self-rescue there is a simplified model of people movement in the tunnel. The model takes into consideration the elementary movement parameters as: start of the self-rescue, walking speed, tunnel length and logical curiosity that in the initial location north or south arranges the movement direction north or south. With this the possibility of a tunnel user approaching the fire during the self-rescue procedure is excluded in the model. With a program the self-rescue procedure is defined with the

( )

=

*else*

*END IF*

*IF XL OUT then*

' ' 'OUT' *<sup>t</sup>*

Users' location N-North or S-South from the fire

−

( )

*ABS X XL*

*poz t*

ASDL - Average smoke density value in layer [mg/m3]

SLH – Smoke layer height [m]

**8.4 Evacuation model** 

following conditional note:

1

*t*

−

' '

0

*END I*

*END IF*

Fig. 7. Evacuation model

*else*

1

*t*

−

*XL v t t*

*F*

+ ⋅−

1

*t*

−

*else*

*else*

 

*END IF*

ASLH – Allowed smoke layer height [m] AT – Average temperature in profile [oC] ATL – Average temperature in layer [oC] TLH – Temperature layer height [m]

ATLH – Allowed temperature layer height [m]

( ) ( )

*IF X S XL v t t then*

= ∧ − ⋅− ≤

*IF X S then XL v t t*

0

*SJ H SR*

( )

*H SR*

( ) ( )

*SJ H SR*

= − ⋅−

( ) ( )

'' 'OUT'

*IF X J XL v t t Lp then*

=∧ + ⋅− ≥

0

*SJ H SR*

0

' ' 0 'OUT'

The first level of risk is presented by the presence of smoke that includes the first four risk stages, the presence of high temperature contributes additional (the highest) risk stage. The Table 6 presents a deterministic register of risk for a constant location in a tunnel during a fire that is 252 m north of the fire. The picture that we get whit this is very representative because it confirms the theory on safety analyses. From the table the safety categories can be seen and appropriate consequences can be allocated. Table 6 has especially a comparative purpose for finding the influence of different ways of ventilation on the fire dynamic, smoke and temperature development. We can logically assume that the risk in low fire force is lower in comparison with bigger fires. Following the same logic along with the consideration of different ways of ventilation and manner of management it soon becomes difficult. One of the noticeable differences is the level of the calculated risk (MR – medium risk) in longitudinal ventilation of a 20 MW fire. It is expected that the increased risk also appears in the 50 MW fire but it is not so. The search for a cause is difficult because this is hidden in the fluid dynamics during the fire, taking into account that the geometry, the discreetness and the initial and boundary conditions (except the force of the fire) are unaltered.

The second important result in the table is the possibility of analysing the influence of turning on the ventilators on the forming of the smoke curtain. It is especially noticeable in the transverse ventilation of 50 MW and 100MW fires where on the turning on, local increased temperatures and smoke concentrations occur. Table 6 is made as a functionally dependent dynamic register, which chooses the calculated values from the data base with the changing of the observed location and on the basis of conditions from the chapter above, calculates the risk.

With the evacuation model is possible to define the tunnel user's location in time intervals of one minute on the bases of the starting user's position, delay with the self-rescue procedure and the walking speed. In this way it is possible to predict the tunnel users movement for the following 15 min and check the smoke concentration and the temperature height to which they will be exposed or establish the risk level. The table represents a conceptual


Table 6. The deterministic risk register for the chosen observer location

Fluid Dynamic Models Application in Risk Assessment 87

model for a general presentation of the risk in tunnels with different types of ventilation and

For a register of a fire in a tunnel and a safety evaluation, the probability accession is too general because a greater event number of physical legality is shown with a statistical probability. A relatively accurate fire dynamics register, which is possible with mathematical models, is often meant only for science. That is why the dissertation is ideally oriented in the use of mathematical CFD models for the making of a system of scenarios that can be further used for developing an effective fire plan or fire management, fire drills, etc. A complex of fire scenarios in different tunnel ventilations and fire forces is presented in the work. The work includes a qualitative analysis of the current circumstances in four different ventilation conditions; natural, longitudinal, transverse and semi transverse. In this way a comparison of individual types of ventilation systems, ventilation plans and their effectiveness in assuring sufficient evacuation times is possible. Also a possibility of usage on a singular tunnel is presented, for which a deterministic safety analysis within a selected number of scenarios would be made. Such an accession requires a lot of calculating time but it is changeable in the development of the safety analysis and fire plan. The geometry and some ventilation

plans are "constants" in this case and only the fire location can be changed.

accessible data of the situation during a fire in a tunnel.

facilities, Fire Safety Journal 36 (2001) 597-619.

Springer-VerlagBerlin Heidelberg, 1991.

Edinburgh, 17-19 September 2001:767-778.

Apart from the number and way of setting the scenarios, the simulation results are values of the selected variables. The discussed variables are mostly the smoke density and the temperature which define the different risk levels on the basis of the human endurance in increased values and conditional interacting dependence. On this basis the deterministic risk register is made which is the key element of the dissertation. The register presents a passage between a practical way of using the CFD model and the user who needs clear and fast

Brussaard, L.A., Kruiskamp, M.M., Oude Essink M.P., The Dutch model for the

Cheng, L.H., T.H. Ueng, C.W. Liu, Simulation of ventilation and fire in the underground

Commision of the European Communities, Directive of the European Parliament and of the

Fletcher, C.A.J., Computational Techniques of Fluid Dynamics-Volume 2, Second edition,

Floyd, J.E., Wieczorek, C.J., Vandsburger, U., 2001. Simulation of the Virginia tech fire

Quantitative Risk Analyses of roadd tunnel, Ministry of Transport, The

Council on minimum safety requirements for tunnels in the Trans-European Road

research laboratory using Large Eddy Simulations with mixture fraction chemistry and finite volume radiative heat transfer, INTERFLAM 2001; Proc. intern. Symp.,

different fire forces.

**10. References** 

Netherlands 2004

Network, Brussels, 30.04.2004.

**9. Conclusion on tunnel risk assessment** 

Table 6. The deterministic risk register for the chosen observer location

model for a general presentation of the risk in tunnels with different types of ventilation and different fire forces.

#### **9. Conclusion on tunnel risk assessment**

For a register of a fire in a tunnel and a safety evaluation, the probability accession is too general because a greater event number of physical legality is shown with a statistical probability. A relatively accurate fire dynamics register, which is possible with mathematical models, is often meant only for science. That is why the dissertation is ideally oriented in the use of mathematical CFD models for the making of a system of scenarios that can be further used for developing an effective fire plan or fire management, fire drills, etc. A complex of fire scenarios in different tunnel ventilations and fire forces is presented in the work. The work includes a qualitative analysis of the current circumstances in four different ventilation conditions; natural, longitudinal, transverse and semi transverse. In this way a comparison of individual types of ventilation systems, ventilation plans and their effectiveness in assuring sufficient evacuation times is possible. Also a possibility of usage on a singular tunnel is presented, for which a deterministic safety analysis within a selected number of scenarios would be made. Such an accession requires a lot of calculating time but it is changeable in the development of the safety analysis and fire plan. The geometry and some ventilation plans are "constants" in this case and only the fire location can be changed.

Apart from the number and way of setting the scenarios, the simulation results are values of the selected variables. The discussed variables are mostly the smoke density and the temperature which define the different risk levels on the basis of the human endurance in increased values and conditional interacting dependence. On this basis the deterministic risk register is made which is the key element of the dissertation. The register presents a passage between a practical way of using the CFD model and the user who needs clear and fast accessible data of the situation during a fire in a tunnel.

#### **10. References**


Fluid Dynamic Models Application in Risk Assessment 89

Center for Chemical Process Safety, Guidelines for the management of change for process

D.W. Hissong, Keys to modeling LNG spills on water, Journal of Hazardous Materials 140

Det Norske Veritas, MPACT Theory manual, Internal publication, DNV software, June

Det Norske Veritas, Risk Assessment LNG import Koper: nautical and unloading

Fells, I., Rutherford, A. G., Burning velocity of methane-air flames, Combustion and Flame,

Gucma L. "Evaluation of oil spills in the Baltic Sea be means of simulation model and

Handleiding risicoberekeningen Bevi: inleiding, Module A/B/C - Versie 3.0", RIVM,

McGrattan, K., Baum, H., Rehm, R., Hamins, A., Forney, G.P., Floyd, J.E. and Hostikka, S.,

Melhem, G. A., Kalelkar, A. S., Understanding LNG Fire Hazards, ioMosaic Corporation

Michael Hightower, M., Luketa-Hanlin, A., Gritzo, L. A., Covan, J. M., Sandia National

Perkovic, M., Gucma, L., Przywarty, M., Gucma, M., Petelin, S., Vidmar, P., Nautical risk

Petelin S., Vidmar P., Perkovič M., Luin B., Kožuh M., Predlog prometno-varnostnih analiz

Pitblado, R. M., Baik, J., Hughes, G. J., Ferro, C., Shaw, S. J., Consequences of LNG Marine

Raj, P. K. , Lemoff, T., Risk analysis based LNG facility siting standard in NFPA 59A, Journal

Safer System LLC, User's Guide Trace 8-Description of modelling algorithms, Westlake

SFPE Handbook, Fire protection engineering, 2nd edition, National Fire Protection

Vidmar, P., Petelin, S., An analysis of a fire resulting from a traffic accident, Journal of

International Maritime Organisation, IMO Guideline for alternative tanker design, 1995. ioMosaic, Modelling LNG Spreading and Vaporisation, ioMosaic Corporation 2007

statistical data. International Maritime Association of Mediterranean", Kolev and

2001. Fire Dynamics Simulator - Technical reference guide, National Institute of

Laboratories, Review of the independent risk assessment of the proposed Cabrillo

assessment for LNG operations at the Port of Koper, Inernational conference on

safety, John Wiley & Sons, Inc., Hoboken, New Jersey, 2008

operations, Report no/DNV Reg No.: / 124UI0A-4, 2008 Drysdale Dougal, An introduction to Fire dinamics, John Wiley & Sons Ltd - 1998

Volume 13, Issue 2, April 1969, Pages 130-138

Standard and Technology, NISTIR 6783, 2001.

liquefied natural gas deepwater port project, 2006

za plinski terminal (Sovenian only), Portorož 2009

Incidents, Det Norske Veritas (USA) Inc., Houston, 2004.

Trbojevic, V. M., Risk criteria in EU, ESREL'05, Poland, 27-30 June 2005

Mechanical Engineering 49(2003), ISSN 0039-2480, pp. 1-13

of Loss Prevention in the Process Industries 22 (2009) 820-829

Nautical Risk assessment LNG transport Rostock", DNV Energy, December 2007

Soares editors), Balkema 2007.

traffic science, Portorož 2010

Village, California , USA-1996

Association 1995

(2007) 465-477

January 2008 http://www.hse.gov.uk/

2007

2007


#### **References on LNG**


Gann, R.G., Hall, J.R., Fire Conditions for Smoke Toxicity Measurement, Fire and Materials,

Gasser, I., Struckmeier, J., Modelling and simulation of fires in vehicle tunnels,

Haack A., Fire Protection in Traffic Tunnels: General and Results of the EUREKA

Haack, A., Current safety issues in traffic tunnels, Tunnelling and Underground Space

Heskestad, G., SFPE Handbook, chapter Fire Plumes. National Fire Protection Association,

Jang, H.M., Chen, F. A novel approach to the transient ventilation of road tunnels, Journal of

Kirchsteiger, C., On the use of probabilistic and deterministic methods in risk analysis, Journal of Loss Prevention in the Process Industries 12 (1999) 399-419 Kunsch J.P., Critical velocity and range of a pre-gas plume in a ventilated tunnel,

Lesieur, M., Turbulence in Fluids, Kluwer Academic Pudlisher, Netherlands 1997, ISBN: 0-

McGrattan, K., Baum, H., Rehm, R., Hamins, A., Forney, G.P., Floyd, J.E. and Hostikka, S.,

Memorial Trunnel Fire Ventilation Test Program, Massachusetts Highway Department,

Mingchung, L., Vaughan, B., Stoichiometric combustion model with oxygen threshold improved predictionsfor fire simulation using CFD model, pp. 559570, 1999. Persson, M., Quantitative Risk Analysis, Procedure for the Fire Evacuation of a Road

PIARC Technical Committee on Tunnel Operation, Fire and Smoke Control in Road

Sagaut, P., Large Eddy Simulations for Incompressible Flows, second edition, Springer

Vidmar, P., Petelin, S., An analysis of a fire resulting from a traffic accident, Journal of

Woodburn P. J., Britter R. E., CFD Simulations of a Tunnel Fire, Fire Safety Journal 16 (1996)

Bottelberghs, P. H. (2000). Risk analysis and safety policy developments in the Netherlands.

Bubbico, R., Di Cave, S., Mazzarotta, B., Preliminary risk analysis for LNG tankers

Industries, Volume 22, Issue 5, September 2009, Pages 634-638.

approaching a maritime terminal, Journal of Loss Prevention in the Process

Mechanical Engineering 49(2003), ISSN 0039-2480, pp. 1-13

2001. Fire Dynamics Simulator - Technical reference guide, National Institute of

Tunnel, Department of Fire Safety Engineering, Lund University, Sweden, Lund

Wind Engineering and Industrial Aerodynamics 86 (2000) 15-36.

Project, Tunnelling and Underground Space Technology, Vol. 13, No. 4, pp. 377-

Computational fluid dynamics and data analyses, November 2002

vol.18, 193-199, 1994

Technology 17 p.p. 117-127 - 2002.

Quincy, Massachusetts, 2nd edition, 1995.

Atmospheric Environment 33 (1999) pp.: 13-24.

Standard and Technology, NISTIR 6783, 2001.

Central Artery/Tunnel Project, USA-1996

381.1998.

7923-4415-4(HB).

Tunnels, France 2003

Berlin Heidelberg 2002.

Journal of Hazardous Materials, 71

2002

35-62. **References on LNG** 


**5** 

*Japan* 

**Sail Performance Analysis of Sailing Yachts by** 

Sails of a sailing yacht can be considered as multiple soft thin wings (membrane wings) with relative large cambers, and are often used at large attack angles. The shape of sail is determined as an equilibrium state of both aerodynamic force and tension acting on the sail surface. In particular, a spinnaker used for the running condition is a very soft membrane like a parachute, and the shape is simply formed by self-generated aerodynamic forces which are strongly affected by the sail shape itself. These facts lead to new challenges in the present problem, i.e., in the measurements the sail shape must be accurately measured in the flying condition, and in numerical simulation of flow and forces the sail flying shape is correctly given or predicted as a part of solution. The present study concerns the authors' ongoing effort on analyses of sail performance of sailing yachts by numerical calculations and experiments, and in this paper, the focus of discussions is more on the former. Two computational fluid dynamics (CFD) methods are used in the present study, and the results are validated through detailed comparison with experimental data. The data are obtained in onboard full-scale measurements by using a sail dynamometer boat. Our study concerns both the upwind and downwind sailing conditions; however, we focus on the former in the present chapter due to the limitation of space in this book. More detailed background of the

One of the two CFD methods is a Vortex Lattice method (VLM). Although the VLM is a potential flow calculation, it is well known the results agree well with the measured data at the upwind condition of small attack angle. The VLM is used as the sail design and making tool due to the quick convergence ability for the parametric survey of sail shape to obtain the desired sail performance, and also due to good compatibility with the finite element method (FEM) for the strength analysis. In this paper, a method to shed wake vortices stepby-step developed by Fukasawa was adopted in the Vortex Lattice method (Fukasawa, 1993;

Another CFD method is a Multiblock Reynolds-Averaged Navier-Stokes (RANS)-based CFD named "*FLOWPACK*". This code was developed by Tahara specifically for CFD education and research, and design applications for ship hydrodynamics, aerodynamics, and fluid engineering (Tahara, 2008). As part of the developments for application to design problems, a complete multiblock domain decomposition feature was included. The numerical method of *FLOWPACK* solves the unsteady RANS and continuity equations for mean velocity and pressure. Either a zero or a two-equation turbulence model can be used for turbulence flow calculation, and in the present study the former was used. The *FLOWPACK* was included as a

present work is well described in Masuyama et al. (2009).

**1. Introduction**

Fukasawa & Katori, 1993).

**Numerical Calculations and Experiments** 

*National Maritime Research Institute, Kanazawa Institute of Technology* 

Y. Tahara, Y. Masuyama, T. Fukasawa and M. Katori

*Osaka Prefecture University, North Sails* 

Vidmar, P., Petelin, S., Analysis of the effect of an external fire on the safety operation of a power plant, Fire Safety Journal 41 (2006) 486-490

World LNG Carrier Fleet, LNG Journal, Maritime Content Ltd, July/August 2008

## **Sail Performance Analysis of Sailing Yachts by Numerical Calculations and Experiments**

Y. Tahara, Y. Masuyama, T. Fukasawa and M. Katori *National Maritime Research Institute, Kanazawa Institute of Technology Osaka Prefecture University, North Sails Japan* 

#### **1. Introduction**

90 Fluid Dynamics, Computational Modeling and Applications

Vidmar, P., Petelin, S., Analysis of the effect of an external fire on the safety operation of a

World LNG Carrier Fleet, LNG Journal, Maritime Content Ltd, July/August 2008

power plant, Fire Safety Journal 41 (2006) 486-490

Sails of a sailing yacht can be considered as multiple soft thin wings (membrane wings) with relative large cambers, and are often used at large attack angles. The shape of sail is determined as an equilibrium state of both aerodynamic force and tension acting on the sail surface. In particular, a spinnaker used for the running condition is a very soft membrane like a parachute, and the shape is simply formed by self-generated aerodynamic forces which are strongly affected by the sail shape itself. These facts lead to new challenges in the present problem, i.e., in the measurements the sail shape must be accurately measured in the flying condition, and in numerical simulation of flow and forces the sail flying shape is correctly given or predicted as a part of solution. The present study concerns the authors' ongoing effort on analyses of sail performance of sailing yachts by numerical calculations and experiments, and in this paper, the focus of discussions is more on the former. Two computational fluid dynamics (CFD) methods are used in the present study, and the results are validated through detailed comparison with experimental data. The data are obtained in onboard full-scale measurements by using a sail dynamometer boat. Our study concerns both the upwind and downwind sailing conditions; however, we focus on the former in the present chapter due to the limitation of space in this book. More detailed background of the present work is well described in Masuyama et al. (2009).

One of the two CFD methods is a Vortex Lattice method (VLM). Although the VLM is a potential flow calculation, it is well known the results agree well with the measured data at the upwind condition of small attack angle. The VLM is used as the sail design and making tool due to the quick convergence ability for the parametric survey of sail shape to obtain the desired sail performance, and also due to good compatibility with the finite element method (FEM) for the strength analysis. In this paper, a method to shed wake vortices stepby-step developed by Fukasawa was adopted in the Vortex Lattice method (Fukasawa, 1993; Fukasawa & Katori, 1993).

Another CFD method is a Multiblock Reynolds-Averaged Navier-Stokes (RANS)-based CFD named "*FLOWPACK*". This code was developed by Tahara specifically for CFD education and research, and design applications for ship hydrodynamics, aerodynamics, and fluid engineering (Tahara, 2008). As part of the developments for application to design problems, a complete multiblock domain decomposition feature was included. The numerical method of *FLOWPACK* solves the unsteady RANS and continuity equations for mean velocity and pressure. Either a zero or a two-equation turbulence model can be used for turbulence flow calculation, and in the present study the former was used. The *FLOWPACK* was included as a

Sail Performance Analysis of Sailing Yachts by Numerical Calculations and Experiments 93

I, J, P, E of Sail are defined in Fig. 1.

**3. Overview of Vortex Lattice method (VLM)** 

Fig. 2. 2-Dimensional flat plate wing and a vortex filament

**3.1 Basic concept of Vortex Lattice method** 

Table 1. Principal dimensions of *Fujin* and detailed measurements of sails

force acting on the wing is approximated by the force acting on the vortex line.

The Vortex Lattice method is a branch of CFD, and it is often used at the early stage of yacht sail design because of the comparatively less computational time. This method is based on the potential theory, similar to the panel method, and the flow around the sail is expressed by discrete vortices. The Vortex Lattice method has its root in the lifting line theory formulated by Prandtl in 1918. A wing is represented by a single vortex line in the lifting line theory, and the

Firstly, 2-dimensional flow around a flat plate wing is considered. A vortex filament is located at a distance " a " from the leading edge of the wing as shown in Fig.2. Although the onset flow U is constant, the flow over the wing is accelerated, while it is decelerated below the wing, because of the flow induced by the vortex filament. This leads to the pressure decrease on the back surface and the pressure increase on the front surface of the wing accordingly to Bernoulli's theorem. This means that the flow around the wing can be realized by a vortex filament in the flow. The strength of vortex filament, or sometimes called circulation, is determined by a boundary condition on the wing; that is, there is no cross flow through the

HULL SAIL DIMENSIONS Length Over All [m] 10.35 Mainsail 130% Jib Length Water Line [m] 8.80 Peak Height [m] 13.82 10.70 Breadth Maximum [m] 3.37 Luff Length [m] 12.50 11.45 Breadth Water Line [m] 2.64 Foot Length [m] 4.44 4.89 Displacement [ton] 3.86 Sail Area [m2] 33.20 26.10

> SAIL Height [%] Chord Length [m] I [m] 11.00 0 4.44 4.89 J [m] 3.61 10 4.13 4.44 P [m] 12.55 20 3.85 3.94 E [m] 4.51 40 3.23 2.94

> > 60 2.43 1.97 80 1.39 0.98 100 0.15 0.10

solver in a sail performance analyzer named "*Advanced Aero Flow* (*AAF*)" developed by Katori (Katori, 2009). The AAF is a specialized package for the calculation of sail performance of sailing yachts, and composed of both mesh generator and post analyzer.

The sail shapes and performance were measured using a sail dynamometer boat *Fujin* under sailing condition on the sea (Masuyama et al., 1997a, 1997b). *Fujin* is a 34-foot LOA boat, in which load cells and CCD cameras were installed to simultaneously measure the sail forces and shapes. At the same time, the sailing conditions of the boat, e.g., boat speed, heel angle, wind speed, and wind angle, were measured. The shapes and 3D coordinates of the sails were used for the input data of the numerical calculations, and the calculated results were compared with the measured data. The sail coordinates with aerodynamic coefficients are tabulated for some sailing conditions in order to provide benchmark data for the CFD validation.

In this paper, overview of the above-mentioned CFD methods and experiments are described. As the aforementioned, sail flying shapes are considered in the present CFD so that the accurate prediction of flow and aerodynamic forces is possible. Discussion of the results is based on detailed comparison with the measurements. The discussion also includes the current capability of the CFD methods in the present problem, and prognosis for the enhancement of the capability in future work for higher accuracy and/or more complicated flow simulation. It will be noteworthy that the overall trends of the flow and the aerodynamic forces measured in the experiments are fairly well predicted by the present computations; and at the same time, experimental techniques originally implemented and used in the present study are shown very promising and capable to provide very detailed benchmark data for CFD validation.

#### **2. Sail plan for the analysis**

In this study the experiments and numerical calculations were performed for the upwind sailing condition. The sail shapes and performance were measured using a sail dynamometer boat *Fujin*. The sail plan of the *Fujin* and the coordinate system are shown in Fig.1. The principal dimensions of the boat and the detailed measurements of the sails are also shown in Table 1. The measurement system of the boat and testing conditions are described in section 5, and the measured and calculated results are compared and discussed in section 6.

Fig. 1. Schematic showing the sail plan of *Fujin* with 130% jib and the coordinate system


Table 1. Principal dimensions of *Fujin* and detailed measurements of sails

#### **3. Overview of Vortex Lattice method (VLM)**

#### **3.1 Basic concept of Vortex Lattice method**

92 Fluid Dynamics, Computational Modeling and Applications

solver in a sail performance analyzer named "*Advanced Aero Flow* (*AAF*)" developed by Katori (Katori, 2009). The AAF is a specialized package for the calculation of sail performance of

The sail shapes and performance were measured using a sail dynamometer boat *Fujin* under sailing condition on the sea (Masuyama et al., 1997a, 1997b). *Fujin* is a 34-foot LOA boat, in which load cells and CCD cameras were installed to simultaneously measure the sail forces and shapes. At the same time, the sailing conditions of the boat, e.g., boat speed, heel angle, wind speed, and wind angle, were measured. The shapes and 3D coordinates of the sails were used for the input data of the numerical calculations, and the calculated results were compared with the measured data. The sail coordinates with aerodynamic coefficients are tabulated for

In this paper, overview of the above-mentioned CFD methods and experiments are described. As the aforementioned, sail flying shapes are considered in the present CFD so that the accurate prediction of flow and aerodynamic forces is possible. Discussion of the results is based on detailed comparison with the measurements. The discussion also includes the current capability of the CFD methods in the present problem, and prognosis for the enhancement of the capability in future work for higher accuracy and/or more complicated flow simulation. It will be noteworthy that the overall trends of the flow and the aerodynamic forces measured in the experiments are fairly well predicted by the present computations; and at the same time, experimental techniques originally implemented and used in the present study are shown very promising and capable to provide very detailed

In this study the experiments and numerical calculations were performed for the upwind sailing condition. The sail shapes and performance were measured using a sail dynamometer boat *Fujin*. The sail plan of the *Fujin* and the coordinate system are shown in Fig.1. The principal dimensions of the boat and the detailed measurements of the sails are also shown in Table 1. The measurement system of the boat and testing conditions are described in section 5,

and the measured and calculated results are compared and discussed in section 6.

Fig. 1. Schematic showing the sail plan of *Fujin* with 130% jib and the coordinate system

some sailing conditions in order to provide benchmark data for the CFD validation.

sailing yachts, and composed of both mesh generator and post analyzer.

benchmark data for CFD validation.

**2. Sail plan for the analysis** 

The Vortex Lattice method is a branch of CFD, and it is often used at the early stage of yacht sail design because of the comparatively less computational time. This method is based on the potential theory, similar to the panel method, and the flow around the sail is expressed by discrete vortices. The Vortex Lattice method has its root in the lifting line theory formulated by Prandtl in 1918. A wing is represented by a single vortex line in the lifting line theory, and the force acting on the wing is approximated by the force acting on the vortex line.

Fig. 2. 2-Dimensional flat plate wing and a vortex filament

Firstly, 2-dimensional flow around a flat plate wing is considered. A vortex filament is located at a distance " a " from the leading edge of the wing as shown in Fig.2. Although the onset flow U is constant, the flow over the wing is accelerated, while it is decelerated below the wing, because of the flow induced by the vortex filament. This leads to the pressure decrease on the back surface and the pressure increase on the front surface of the wing accordingly to Bernoulli's theorem. This means that the flow around the wing can be realized by a vortex filament in the flow. The strength of vortex filament, or sometimes called circulation, is determined by a boundary condition on the wing; that is, there is no cross flow through the

Sail Performance Analysis of Sailing Yachts by Numerical Calculations and Experiments 95

and the total velocity of inflow into the wing changes to Ue . If the lift acting on the wing is defined as the force perpendicular to the onset flow direction, it is given by the following

In this case, the force in the onset flow direction is generated, which is the apparent drag

The induced drag is a distinctive drag in a 3-dimensional wing, and does not appear in a 2-

In the Vortex Lattice method, the lift, induced drag, and center of pressure are calculated by arranging horseshoe vortices of different strength on the surface of sail. By placing a number of horseshoe vortices, the sail of complex shape with twist, camber, or two or more sails, can be analyzed. Falkner used the name "Vortex Lattice" firstly in his report, in which a wing was covered with a grid of straight horseshoe vortices (Falkner, 1943, 1946). In 1950's, only the analysis where the trailing vortices are placed in the straight line was able to be carried out because of the computer capability, and the accuracy was questionable. It was 1965 when the Vortex Lattice method started to demonstrates its ability along with the development of computer, and the method came to be used for the performance prediction of yacht sail. The yacht sail is one of the most suitable objects for applying the Vortex Lattice method because of its thickness, if the viscous effect of fluid can be disregarded. An application of the Vortex Lattice method to the performance prediction of yacht sail will be explained in the following paragraph with the use of a step-by-step procedure to estimate

Discretized horseshoe vortices are located on the sail plane in the Vortex Lattice method. It is usual to divide the sail plane into quadrilateral panels as shown in Fig.4, and the

<sup>e</sup> i e i ei i <sup>L</sup> <sup>L</sup> D L sin<sup>α</sup> <sup>L</sup> <sup>α</sup> w w

L L cos e ie e α L ρU Γ ρUΓ (3)

U U *<sup>i</sup>* (4)

formula according to Kutta-Joukowski theorem.

called induced drag given by

the trailing vortex deformations.

Fig. 4. Panel discretization of sails

dimensional wing.

wing. This boundary condition is usually satisfied at a certain point called control point. Assuming that the control point is located at a distance " b " from the leading edge of the wing shown in Fig.2, and satisfying the boundary condition at this point, the strength of the vortex filament can be determined. Once the strength of vortex filament is determined, the lift acting on the vortex filament can be calculated according to Kutta-Joukowski theorem, that is,

$$\mathbf{L} = \mathbf{p} \,\mathrm{U}\Gamma\tag{1}$$

where ρ is the density of the fluid and Γ is the strength of vortex filament. If the calculated lift is assumed to equal that generated in a 2-dimentional thin parabolic shape airfoil, the locations of the vortex filament and the control point are determined to be a 4 *c* and b 34 *c* , where " c " is the chord length of the wing. This is called 1/4-3/4 rule, which was shown by Pistolesi (Pistolesi, 1937). This rule is used as the basis of the present Vortex Lattice method.

#### **3.2 Application to sail configuration**

As the yacht sail is a 3-dimension shape body, attention should be paid to the treatment of the end of vortex line. According to the Helmholtz's theorem on vortex, the vortex line should expand from the boundary to the boundary of the flow or shuts oneself and makes vortex ring. Accordingly, in the 3-dimensional body, the vortex line should expand infinity from the edge of the body. In the lifting line theory, or the Vortex Lattice method, the vortex line is assumed to be a horseshoe type shown in Fig.3(a), and the vortex line changes its direction at the edge of the body to extend to infinity as trailing free vortices.

Fig. 3. Horseshoe vortex and downwash effect

In case of the horseshoe vortex shown in Fig.3(a), the flow induced by each vortex line affects the onset flow in the magnitude and the direction. According to Biot-Savart law, the velocity vector induced by a slight part of the vortex line d is given by

$$\mathbf{v} = \frac{\Gamma}{4\,\mathrm{m}} \iint \frac{\mathbf{d}\,\vec{\ell} \times \vec{\mathbf{r}}}{\left|\vec{\mathbf{r}}\right|^3} \tag{2}$$

where r is a position vector from the vortex part to the point concerned. The downward velocity called downwash wi is calculated by using Equation (2), which affects the onset flow. This causes the reduction of the attack angle of total inflow by αi as shown in Fig.3(b),

wing. This boundary condition is usually satisfied at a certain point called control point. Assuming that the control point is located at a distance " b " from the leading edge of the wing shown in Fig.2, and satisfying the boundary condition at this point, the strength of the vortex filament can be determined. Once the strength of vortex filament is determined, the lift acting

where ρ is the density of the fluid and Γ is the strength of vortex filament. If the calculated lift is assumed to equal that generated in a 2-dimentional thin parabolic shape airfoil, the locations of the vortex filament and the control point are determined to be a 4 *c* and b 34 *c* , where " c " is the chord length of the wing. This is called 1/4-3/4 rule, which was shown by Pistolesi (Pistolesi, 1937). This rule is used as the basis of the present Vortex

As the yacht sail is a 3-dimension shape body, attention should be paid to the treatment of the end of vortex line. According to the Helmholtz's theorem on vortex, the vortex line should expand from the boundary to the boundary of the flow or shuts oneself and makes vortex ring. Accordingly, in the 3-dimensional body, the vortex line should expand infinity from the edge of the body. In the lifting line theory, or the Vortex Lattice method, the vortex line is assumed to be a horseshoe type shown in Fig.3(a), and the vortex line changes its

In case of the horseshoe vortex shown in Fig.3(a), the flow induced by each vortex line affects the onset flow in the magnitude and the direction. According to Biot-Savart law, the

> 4π r

velocity called downwash wi is calculated by using Equation (2), which affects the onset flow. This causes the reduction of the attack angle of total inflow by αi as shown in Fig.3(b),

3 Γ d r

is a position vector from the vortex part to the point concerned. The downward

(2)

velocity vector induced by a slight part of the vortex line d is given by

v

direction at the edge of the body to extend to infinity as trailing free vortices.

L ρUΓ (1)

on the vortex filament can be calculated according to Kutta-Joukowski theorem, that is,

Lattice method.

where r

**3.2 Application to sail configuration** 

Fig. 3. Horseshoe vortex and downwash effect

and the total velocity of inflow into the wing changes to Ue . If the lift acting on the wing is defined as the force perpendicular to the onset flow direction, it is given by the following formula according to Kutta-Joukowski theorem.

$$\mathbf{L} = \mathbf{L\_e} \cos \mathbf{u\_i} \approx \mathbf{L\_e} = \mathbf{p} \, \mathbf{U\_e} \, \Gamma \approx \mathbf{p} \, \mathbf{U} \, \Gamma \tag{3}$$

In this case, the force in the onset flow direction is generated, which is the apparent drag called induced drag given by

$$\mathbf{L}\mathbf{D}\_{\mathrm{i}} = \mathbf{L}\_{\mathrm{e}} \sin \mathfrak{a}\_{\mathrm{i}} \approx \mathbf{L}\_{\mathrm{e}} \mathbf{a}\_{\mathrm{i}} = \frac{\mathbf{L}\_{\mathrm{e}}}{\mathbf{U}} \mathbf{w}\_{\mathrm{i}} \approx \frac{\mathbf{L}}{\mathbf{U}} \mathbf{w}\_{\mathrm{i}} \tag{4}$$

The induced drag is a distinctive drag in a 3-dimensional wing, and does not appear in a 2 dimensional wing.

In the Vortex Lattice method, the lift, induced drag, and center of pressure are calculated by arranging horseshoe vortices of different strength on the surface of sail. By placing a number of horseshoe vortices, the sail of complex shape with twist, camber, or two or more sails, can be analyzed. Falkner used the name "Vortex Lattice" firstly in his report, in which a wing was covered with a grid of straight horseshoe vortices (Falkner, 1943, 1946). In 1950's, only the analysis where the trailing vortices are placed in the straight line was able to be carried out because of the computer capability, and the accuracy was questionable. It was 1965 when the Vortex Lattice method started to demonstrates its ability along with the development of computer, and the method came to be used for the performance prediction of yacht sail. The yacht sail is one of the most suitable objects for applying the Vortex Lattice method because of its thickness, if the viscous effect of fluid can be disregarded. An application of the Vortex Lattice method to the performance prediction of yacht sail will be explained in the following paragraph with the use of a step-by-step procedure to estimate the trailing vortex deformations.

Discretized horseshoe vortices are located on the sail plane in the Vortex Lattice method. It is usual to divide the sail plane into quadrilateral panels as shown in Fig.4, and the

Fig. 4. Panel discretization of sails

Sail Performance Analysis of Sailing Yachts by Numerical Calculations and Experiments 97

parallel to the velocity field induced by total vortex system. A step-by-step procedure developed by Fukasawa was adopted in this paper to determine the strength of the bound vortices and the location of wake vortices (Fukasawa, 1993; Fukasawa & Katori, 1993; Masuyama et al., 1997a, 1997b). The wake vortices are shed from the trailing edge in each

<sup>D</sup>Γ<sup>B</sup> <sup>0</sup>

where ΓB and ΓW are the total strength of bound vortices, or the circulation around the

Γ Γ B B <sup>U</sup> t ξ 

Γ Γ B W <sup>U</sup> t ξ 

where ξ is taken to the downstream direction, and U is the local velocity at the wave vortex. Assuming that the wake vortices proceeds Δξ downstream in a time step Δt with the velocity U , the strength of wake vortex shed at time step k, can be given by integrating

> k B k k-1 k-1 k W BB B B

<sup>1</sup> <sup>Γ</sup> Δξ ΔΓ d<sup>ξ</sup> <sup>Γ</sup> - Γ Γ - <sup>Γ</sup>

Equation (12) means that the strength of wake vortex shed at time step k is the increase of the strength of bound vortex from time step k-1 to time step k. Once a vortex filament is shed at time step k, it proceeds downstream with a constant strength according to the local field velocity, i.e., each horseshoe wake vortex moves in the direction of field velocity in each time step. The field velocity is updated in every time step. The calculation is carried forward until the the calculated lift and drag forces converges. The forces vector and the

<sup>F</sup> ρ Γ U ds (13)

with the use of the vortex strengths of the wake vortices and the bound vortices determined

Finally, the overall numerical solution procedure of the present Vortex Lattice method is

Step 1. Divide the sail planes into quadrilateral panels, and allocate horseshoe vortices on

Step 2. Input the mast rake angle, heel angle of the yacht, apparent wind speed and

U U t Δt

moment acting on the sail are calculated accordingly to Kutta-Joukowski theorem.

(12)

sail, and the strength of wake vortices, respectively. From Equation (9), we have

Δt

0

Γ Γ B W 0 (8)

Dt (9)

(10)

(11)

time step according to Helmholtz's theorem; that is,

and substituting Equation (10) into (8), we have

Equation (11), that is,

by solving Equation (7).

summarized as follows:

the sail plane.

apparent wind angle.

horseshoe vortices are placed at 1/4 length of panel from the front end edge of the panel so as to trail the trailing vortices rearwards. The strengths of the horseshoe vortices are determined by satisfying the boundary condition on the sail; that is, the total flow of the onset flow and the induced wake by vortices is parallel to the sail surface at control points. The control point is taken to be the point 3/4 of length of panel from the front end edge of the panel according to the 1/4-3/4 rule.

According to Biot-Savart law, the velocity vector at the control point of i-th panel induced by other vortices are given by

$$\vec{\mathbf{v}}\_{i} = \frac{1}{4\,\mathrm{m}} \int \frac{\Gamma\left(\mathrm{d}\vec{\ell} \times \vec{\mathbf{r}}\right)}{\left|\,\vec{\mathbf{r}}\right|^{3}} = \sum\_{\text{j}} \frac{\Gamma\_{\text{j}}}{4\,\mathrm{m}} \sum\_{\text{k}=1}^{3} \frac{\cos\alpha\_{\text{k}\vec{\text{r}}} + \cos\beta\_{\text{k}\vec{\text{j}}\,\mathrm{i}}}{\mathbf{h}\_{\text{k}\vec{\text{i}}}} \bar{\mathbf{e}}\_{\text{k}\vec{\text{j}}}\tag{5}$$

where αkji βkji , hkji , are the angles and the distance of a k-th filament of a j-th horseshoe vortex and i-th control point shown in Fig.5, which shows a plane containing the vortex filament and the control point. kji e is a unit vector perpendicular to the plane shown in Fig.5. k=1, 2, 3 in Equation (5) denotes each vortex filament of a horseshoe vortex. With the use of Equation (5), the boundary condition on the control point can be given by

$$\sum\_{\mathbf{i}}^{\text{NB}} \vec{\mathbf{v}}\_{\mathbf{i}} \cdot \vec{\mathbf{n}} = \vec{\mathbf{U}} \cdot \vec{\mathbf{n}}\_{\mathbf{i}} - \sum\_{\mathbf{i}}^{\text{NW}} \vec{\mathbf{v}}\_{\mathbf{i}} \cdot \vec{\mathbf{n}} \tag{6}$$

where U is an onset flow velocity vector and n is the unit normal vector at the control points. NB is the number of bound horseshoe vortices on the sail plane and NW is the number of trailing horseshoe vortices in the wake. Equation (6) can be written in the vector matrix form by

$$\begin{bmatrix} \Lambda \end{bmatrix} \begin{Bmatrix} \Gamma \end{Bmatrix} = \begin{Bmatrix} \mathbf{u\_n} \end{Bmatrix} + \begin{Bmatrix} \mathbf{v\_n} \end{Bmatrix} \tag{7}$$

Solving Equation (7), the strength of bound vortices can be obtained.

Fig. 5. K-th filament of j-th horseshoe vortex and i-th control point

#### **3.3 Step-by-step vortex shedding technique**

The important point in the Vortex Lattice approach to yacht sail is the handling of the wake of the sail. The wake vortices proceed downstream from the trailing edge, or leech/foot of sail, in the Vortex Lattice method. The location of wake vortices are determined by the condition that they are free vortices; that is, the stream line of wake vortices should be parallel to the velocity field induced by total vortex system. A step-by-step procedure developed by Fukasawa was adopted in this paper to determine the strength of the bound vortices and the location of wake vortices (Fukasawa, 1993; Fukasawa & Katori, 1993; Masuyama et al., 1997a, 1997b). The wake vortices are shed from the trailing edge in each time step according to Helmholtz's theorem; that is,

$$
\Gamma\_\mathcal{B} + \Gamma\_\mathcal{W} = \mathcal{O} \tag{8}
$$

$$\frac{\mathrm{D}\Gamma\_{\mathrm{B}}}{\mathrm{D}\mathrm{t}} = 0\tag{9}$$

where ΓB and ΓW are the total strength of bound vortices, or the circulation around the sail, and the strength of wake vortices, respectively. From Equation (9), we have

$$\frac{\partial \Gamma\_{\rm B}}{\partial \mathbf{t}} = \tilde{\mathbf{U}} \frac{\partial \Gamma\_{\rm B}}{\partial \xi} \tag{10}$$

and substituting Equation (10) into (8), we have

96 Fluid Dynamics, Computational Modeling and Applications

horseshoe vortices are placed at 1/4 length of panel from the front end edge of the panel so as to trail the trailing vortices rearwards. The strengths of the horseshoe vortices are determined by satisfying the boundary condition on the sail; that is, the total flow of the onset flow and the induced wake by vortices is parallel to the sail surface at control points. The control point is taken to be the point 3/4 of length of panel from the front end edge of

According to Biot-Savart law, the velocity vector at the control point of i-th panel induced

1 Γ d r Γ cosα cosβ v e 4π r 4π h 

where αkji βkji , hkji , are the angles and the distance of a k-th filament of a j-th horseshoe vortex and i-th control point shown in Fig.5, which shows a plane containing the vortex

Fig.5. k=1, 2, 3 in Equation (5) denotes each vortex filament of a horseshoe vortex. With the

points. NB is the number of bound horseshoe vortices on the sail plane and NW is the number of trailing horseshoe vortices in the wake. Equation (6) can be written in the vector

*kji kji*

j

The important point in the Vortex Lattice approach to yacht sail is the handling of the wake of the sail. The wake vortices proceed downstream from the trailing edge, or leech/foot of sail, in the Vortex Lattice method. The location of wake vortices are determined by the condition that they are free vortices; that is, the stream line of wake vortices should be

hk

i

NB NW i ii i i

use of Equation (5), the boundary condition on the control point can be given by

is an onset flow velocity vector and n

Solving Equation (7), the strength of bound vortices can be obtained.

Fig. 5. K-th filament of j-th horseshoe vortex and i-th control point

**3.3 Step-by-step vortex shedding technique** 

 <sup>3</sup> <sup>j</sup> kji kji i k 3 ji

j k1 kji

(5)

is a unit vector perpendicular to the plane shown in

v n Un v n (6)

Λ Γ u v n n (7)

is the unit normal vector at the control

the panel according to the 1/4-3/4 rule.

by other vortices are given by

filament and the control point. kji e

where U 

matrix form by

$$\frac{\partial \Gamma\_{\rm B}}{\partial \mathbf{t}} = -\mathbf{\tilde{U}} \frac{\partial \Gamma\_{\rm W}}{\partial \boldsymbol{\xi}} \tag{11}$$

where ξ is taken to the downstream direction, and U is the local velocity at the wave vortex. Assuming that the wake vortices proceeds Δξ downstream in a time step Δt with the velocity U , the strength of wake vortex shed at time step k, can be given by integrating Equation (11), that is,

$$
\Delta\Gamma\_{\mathrm{W}}^{\mathrm{k}} = -\frac{1}{\bar{\mathbf{U}}} \int\_{0}^{\Delta\mathbf{t}} \frac{\partial}{\partial \text{ t}} \,\mathrm{d}\xi = -\frac{\Delta\xi}{\bar{\mathbf{U}}\,\Delta\mathbf{t}} \left(\Gamma\_{\mathrm{B}}^{\mathrm{k}} \cdot \Gamma\_{\mathrm{B}}^{\mathrm{k}\cdot\mathrm{1}}\right) = \Gamma\_{\mathrm{B}}^{\mathrm{k}\cdot\mathrm{1}} \cdot \Gamma\_{\mathrm{B}}^{\mathrm{k}}\tag{12}
$$

Equation (12) means that the strength of wake vortex shed at time step k is the increase of the strength of bound vortex from time step k-1 to time step k. Once a vortex filament is shed at time step k, it proceeds downstream with a constant strength according to the local field velocity, i.e., each horseshoe wake vortex moves in the direction of field velocity in each time step. The field velocity is updated in every time step. The calculation is carried forward until the the calculated lift and drag forces converges. The forces vector and the moment acting on the sail are calculated accordingly to Kutta-Joukowski theorem.

$$
\vec{\mathbf{F}} = \mathfrak{p} \int \Gamma \, \vec{\mathbf{U}} \times d\vec{s} \tag{13}
$$

with the use of the vortex strengths of the wake vortices and the bound vortices determined by solving Equation (7).

Finally, the overall numerical solution procedure of the present Vortex Lattice method is summarized as follows:


Sail Performance Analysis of Sailing Yachts by Numerical Calculations and Experiments 99

decomposition feature was included. At present, *FLOWPACK* has a good interface with the authors' inhouse automatic grid generator as well as with commercial grid generation software. For a complete documentation of the method is available in Tahara (2008). In the

Let us consider a sail system fixed in the uniform onset flow (see Fig.1 for the basic coordinate system). The non-dimensional RANS equations for unsteady, three-dimensional

> *j i j ji U U u u <sup>p</sup> U U*

0 *<sup>i</sup> i U x*

where *Ui* (*i*=1,2,3) *=(U,V,W)* and *ui* (*i*=1,2,3) *=*(*u,v,w*) are the Cartesian components of mean

=(*X,Y,Z*) is the dimensionless coordinates normalized by a characteristic length *L*, Re=*U0L*/*ν* is the Reynolds number, *ν* is the kinematic viscosity, the barred quantities *u ui <sup>j</sup>* are the

related to the corresponding mean rate of strain through an isotropic eddy viscosity *νt*, i.e.,

*j i i j t ij j i U U u u k x x*

when a suitable turbulence model is employed to calculate the eddy-viscosity distribution. Either a zero or a two-equation turbulence model can be used for turbulent flow calculation, and a model used for the present study is the former, i.e., Baldwin-Lomax model (Baldwin & Lomax, 1978), which is an algebraic scheme that makes use of a two-layer isotropic eddyviscosity formulation. Detailed validation study of this model for boundary layer flows around three-dimensional bodies was done by the author (Tahara, 1995; Tahara & Stern,

> () ()

*t inner c t outer c*

where *y* is the distance normal to the wall surface and *yc* is the minimum value of *y* where both the inner and outer viscosities match. The inner viscosity follows the Prandtl-Van

*y y y y*

and fluctuating velocities, respectively, normalized by the reference velocity *U0*, *xi*

Reynolds stresses normalized by <sup>2</sup> *U*<sup>0</sup> , and *p* is the pressure normalized by <sup>2</sup>

where *k*=( *uu vv ww* )/2 is the turbulent kinetic energy, Equation (15) becomes

 

*j i ti t j jj j i*

*U U U*

1996). In this model, the eddy viscosity is evaluated as follows:

 

*t xx x x*

*U*

where 1/*R*

=1/Re+*t*, and  <sup>1</sup> <sup>2</sup> <sup>0</sup>

(16)

(17)

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

(18)

(19)

3 *<sup>i</sup> <sup>i</sup> <sup>p</sup> pk U x R*

=*Ui* (*i*=1,2,3). Equations (16) and (18) can be solved for *Ui* and *p*

(*i*=1,2,3)

*U*<sup>0</sup> . If *u ui <sup>j</sup>* are

(15)

Re

2 3

following, an overview of the numerical method is given.

incompressible flow can be written in Cartesian tensor notation as

*i j i i*

*t x x x* 

**4.1 Governing equations** 


Fig.6 shows the example calculation results. In the present study, the mast and rigging were not considered for the series calculations, and the mirror image was taken into account about the deck plane of the boat. Since the vortex lattice methods do not predict viscous drag, the viscous drag acting on the sails and rigging was calculated empirically using a drag coefficient *CDp*. The value of *CDp* was obtained from the measured data in the previous papers and formulated for the upwind condition as follows:

$$\mathbf{C\_{Dp}} = 0.0026 \,\,\mathbf{y\_A} + 0.005 \,\tag{14}$$

where *γA* is apparent wind angle in degrees.

Fig. 6. Calculated wake by Vortex Lattice method

#### **4. Overview of Reynolds-averaged Navier-Stokes equation method**

The RANS-based CFD method used in the present study was *FLOWPACK*. The code was developed by Tahara specifically for CFD education and research and for design applications for ship hydrodynamics, aerodynamics, and fluid engineering. As part of the developments for application to design problems, a complete multiblock domain decomposition feature was included. At present, *FLOWPACK* has a good interface with the authors' inhouse automatic grid generator as well as with commercial grid generation software. For a complete documentation of the method is available in Tahara (2008). In the following, an overview of the numerical method is given.

#### **4.1 Governing equations**

98 Fluid Dynamics, Computational Modeling and Applications

Step 3. Solve the strength of bound horseshoe vortices on the sail plane with non-

Step 5. Compute the increment of the total circulation and shed free horseshoe vortices

Step 8. Solve the strength of bound horseshoe vortices on the sail plane with trailing free

Fig.6 shows the example calculation results. In the present study, the mast and rigging were not considered for the series calculations, and the mirror image was taken into account about the deck plane of the boat. Since the vortex lattice methods do not predict viscous drag, the viscous drag acting on the sails and rigging was calculated empirically using a drag coefficient *CDp*. The value of *CDp* was obtained from the measured data in the previous

C 0.0026 Dp γ<sup>A</sup> 0.005 (14)

Step 4. Compute the total circulation around the sail caused by horseshoe vortices.

Step 6. Calculate the local velocities along the wake, and deform the trailing vortices.

Step 7. Compute the force vector and moment acting on the sails.

Step 9. Repeat Step 4 through Step 8 until the force is converged.

papers and formulated for the upwind condition as follows:

deformed trailing vortices.

according to Equation (12).

where *γA* is apparent wind angle in degrees.

Fig. 6. Calculated wake by Vortex Lattice method

**4. Overview of Reynolds-averaged Navier-Stokes equation method** 

The RANS-based CFD method used in the present study was *FLOWPACK*. The code was developed by Tahara specifically for CFD education and research and for design applications for ship hydrodynamics, aerodynamics, and fluid engineering. As part of the developments for application to design problems, a complete multiblock domain

vortices.

Let us consider a sail system fixed in the uniform onset flow (see Fig.1 for the basic coordinate system). The non-dimensional RANS equations for unsteady, three-dimensional incompressible flow can be written in Cartesian tensor notation as

$$\frac{\partial \mathcal{U}\_i}{\partial t} + \mathcal{U}\_j \frac{\partial \mathcal{U}\_i}{\partial \mathbf{x}^j} + \frac{\partial u\_i u\_j}{\partial \mathbf{x}^j} + \frac{\partial p}{\partial \mathbf{x}^i} - \frac{\mathbf{1}}{\mathbf{Re}} \nabla^2 \mathcal{U}\_i = 0 \tag{15}$$

$$\frac{\partial \mathcal{U}\_i}{\partial \mathbf{x}^i} = \mathbf{0} \tag{16}$$

where *Ui* (*i*=1,2,3) *=(U,V,W)* and *ui* (*i*=1,2,3) *=*(*u,v,w*) are the Cartesian components of mean and fluctuating velocities, respectively, normalized by the reference velocity *U0*, *xi* (*i*=1,2,3) =(*X,Y,Z*) is the dimensionless coordinates normalized by a characteristic length *L*, Re=*U0L*/*ν* is the Reynolds number, *ν* is the kinematic viscosity, the barred quantities *u ui <sup>j</sup>* are the Reynolds stresses normalized by <sup>2</sup> *U*<sup>0</sup> , and *p* is the pressure normalized by <sup>2</sup> *U*<sup>0</sup> . If *u ui <sup>j</sup>* are related to the corresponding mean rate of strain through an isotropic eddy viscosity *νt*, i.e.,

$$-\overline{\boldsymbol{\mu}\_{i}\boldsymbol{\mu}\_{j}} = \nu\_{t} \left( \frac{\partial \boldsymbol{U}\_{i}}{\partial \mathbf{x}^{j}} + \frac{\partial \boldsymbol{U}\_{j}}{\partial \mathbf{x}^{i}} \right) - \frac{2}{3} \boldsymbol{\delta}\_{ij} k \tag{17}$$

where *k*=( *uu vv ww* )/2 is the turbulent kinetic energy, Equation (15) becomes

$$\frac{\partial \mathcal{U}\_i}{\partial t} + \left( \mathcal{U}\_j - \frac{\partial \mathcal{V}\_t}{\partial \mathbf{x}^j} \right) \frac{\partial \mathcal{U}\_i}{\partial \mathbf{x}^j} - \frac{\partial \mathcal{V}\_t}{\partial \mathbf{x}^j} \frac{\partial \mathcal{U}\_j}{\partial \mathbf{x}^i} + \frac{\partial p}{\partial \mathbf{x}^j} \left( p + \frac{2}{3}k \right) - \frac{1}{R\_\phi} \nabla^2 \mathcal{U}\_i = 0 \tag{18}$$

where 1/*R*=1/Re+*t*, and =*Ui* (*i*=1,2,3). Equations (16) and (18) can be solved for *Ui* and *p* when a suitable turbulence model is employed to calculate the eddy-viscosity distribution. Either a zero or a two-equation turbulence model can be used for turbulent flow calculation, and a model used for the present study is the former, i.e., Baldwin-Lomax model (Baldwin & Lomax, 1978), which is an algebraic scheme that makes use of a two-layer isotropic eddyviscosity formulation. Detailed validation study of this model for boundary layer flows around three-dimensional bodies was done by the author (Tahara, 1995; Tahara & Stern, 1996). In this model, the eddy viscosity is evaluated as follows:

$$\begin{cases} (\nu\_t)\_{inner} & \mathcal{Y} \le \mathcal{Y}\_c \\ (\nu\_t)\_{outer} & \mathcal{Y} > \mathcal{Y}\_c \end{cases} \tag{19}$$

where *y* is the distance normal to the wall surface and *yc* is the minimum value of *y* where both the inner and outer viscosities match. The inner viscosity follows the Prandtl-Van

Sail Performance Analysis of Sailing Yachts by Numerical Calculations and Experiments 101

12 13 23 12 13 23 *Ss g g g* 2

Thompson et al. (1985). The transport equations (22) on a computational cell (shown in Fig.7(a)) are linearized and evaluating coefficients and source term at the center node P of

<sup>P</sup> <sup>P</sup> P P

12 3 P P P PP

(b) (c)

*Pg* . The above equation is discretized by the finite-analytic scheme. Solution dependent coefficients are analytically derived by solving the above linearized transport equation using a hybrid method which combines a two-dimensional analytic solution in

conditions on the faces of the cell as a combination of exponential and linear functions, which are the natural solutions for the linearized transport equation, Equation (24) can be solved by the method of separation of variables. When the solution is evaluated at the center node P of

PU D <sup>1</sup>

*<sup>R</sup> CC C <sup>R</sup>*

*n n n nn P EC WC NC SC*

 

NE NW SE SW EC WC NC SC

*C C CC CC C C*

 

*nn n U D*

*CC C C S* 

PU D P P <sup>P</sup>

*n n nn NE NW SE SW*

 

Fig. 7. Definition sketch of a computational cell (a), and nodes in regular grid (b) and

The dimensions of the computational cell are 2*l*×2*k*×2*h*, where *l*=1/ <sup>11</sup>

the element, the following twelve-point finite analytic formula is obtained:

  

 

*g a RS* 

 

 

2 *jj <sup>j</sup> j j <sup>j</sup> <sup>j</sup>*

 

3 2

 

ξ3

1

1 1 2 2 3 3 11 22 33 PP P *ggg*

ξ2

P NC NE EC

SW SE SC

h h

plane with one dimensional analytic solution in the

1

1

 

D

ξ3

WC NW k k

The geometric coefficients *<sup>j</sup>*

the element yields

ξ1

ξ2

U

continuity cell (c).

*h*=1/ <sup>33</sup>

where

or

(a)

222

 

 

 

(23)

 

> 

 *t* (24)

*<sup>i</sup> <sup>b</sup>* , *ij <sup>g</sup>* , and *<sup>j</sup> <sup>f</sup>* appearing in the above equations are defined by

222 *C B A RS*

  *Pg* , *k*=1/ <sup>22</sup>

direction. By specifying boundary

*Pg* , and

 -

(25)

Driest formula, i.e., (t)inner=2|*ω*|, where =*y*[1-exp(-*y*+/*A*+)] is the turbulent length scale for the inner region, and *A*+ are model constants, |*ω*| is the vorticity magnitude, and *y*+ is the dimensionless distance to the wall. In the outer region, eddy viscosity is given by (t)outer=*KC*cp*F*wake*F*Kleb, where *K* and *C*cp are model constants, *F*wake=min(*y*max*F*max, *C*wk*y*max*U*dif2/*F*max), and *F*Kleb= [1+5.5(*C*Kleb*y*/*y*max)6]-1. The *F*max and *y*max are determined by the value and corresponding location, respectively, of the maximum of *F*=*y*|*ω*|[1-exp( *y*+/*A*+)]. The quantity *U*dif is the difference between maximum and minimum velocity magnitudes in the profile and is expressed as *U*dif= 2 2 2 1/2 max ( ) *UVW* - 2 2 2 1/2 min ( ) *UVW* . *C*Kleb and *C*wk are additional model constants. Numerical values for the model constants are *A*+=26, =0.4, *K*=0.0168, *C*cp=1.6, *C*wk=1.0, and *C*Kleb=0.3.

#### **4.2 Discretization and velocity-pressure coupling**

In the following, discretization and velocity-pressure coupling of the present RANS method are described. First, it is convenient to rewrite the transport equations for momentum (*Ui*) in the following general form:

$$\nabla^2 \phi = R\_{\phi} \left[ \sum\_{j=1}^3 \left( \mathcal{U}\_j - \frac{1}{\sigma\_{\phi}} \frac{\partial \mathcal{V}\_t}{\partial \mathbf{x}^j} \right) \frac{\partial \phi}{\partial \mathbf{x}^j} + \frac{\partial \phi}{\partial t} \right] + s\_{\phi} \tag{20}$$

where again represents any one of the convective transport quantities (*Ui*), and *s* is the source function for the corresponding quantity. We transform the physical space (*xi* ,*t*) into a rectangular region in the computational space (*i* ,) using the following coordinate transformations:

$$\begin{aligned} t &= \tau, \quad \mathbf{x}^i = \mathbf{x}^i \left(\boldsymbol{\xi}^j\right), & \mathbf{e}\_i \cdot \nabla \phi &= \frac{1}{J} \sum\_{j=1}^3 b\_i^j \frac{\partial \phi}{\partial \boldsymbol{\xi}^j} \\\\ \nabla^2 \phi &= \sum\_{i=1}^3 \sum\_{j=1}^3 \mathbf{g}^{ij} \frac{\hat{\mathcal{O}}^2 \phi}{\hat{\mathcal{O}} \boldsymbol{\xi}^i \boldsymbol{\xi}^j} + \sum\_{j=1}^3 f^j \frac{\partial \phi}{\hat{\mathcal{O}} \boldsymbol{\xi}^j} \end{aligned} \qquad \begin{aligned} \boldsymbol{e}\_i \cdot \nabla \phi &= \frac{1}{J} \sum\_{j=1}^3 b\_i^j \frac{\partial \phi}{\partial \boldsymbol{\xi}^j} \\ \frac{\partial \phi}{\partial t} &= \frac{\partial \phi}{\partial \boldsymbol{\tau}} - \frac{1}{J} \sum\_{i=1}^3 \sum\_{j=1}^3 b\_i^j \frac{\partial \mathbf{x}^i}{\partial \boldsymbol{\tau}} \frac{\partial \phi}{\partial \boldsymbol{\xi}^j} \end{aligned}$$

Then the continuity equation (16) and the transport equations (20) for momentum parameters can be written as

$$\frac{1}{J} \sum\_{i=1}^{3} \sum\_{j=1}^{3} \frac{\mathcal{C}}{\mathcal{O}\xi^{i}} \left( b\_{j}^{\ \ i} \mathcal{U}\_{j} \right) = 0 \tag{21}$$

$$\sum\_{j=1}^{3} \left( \mathbf{g}^{\circ j} \frac{\hat{\sigma}^{2} \phi}{\hat{\sigma} \xi^{j} \hat{\sigma} \xi^{j}} - \mathbf{2} a\_{\phi}^{j} \frac{\hat{\sigma} \phi}{\hat{\sigma} \xi^{j}} \right) = R\_{\phi} \frac{\hat{\sigma} \phi}{\hat{\sigma} \tau} + S\_{\phi} \tag{22}$$

where

$$\mathfrak{L}a\_{\phi}^{j} = \frac{R\_{\phi}}{J} \sum\_{n=1}^{3} b\_{n}^{j} \left( \mathcal{U}\_{n} - \frac{\partial \mathfrak{x}^{n}}{\partial \, \boldsymbol{\tau}} - \frac{1}{J \sigma\_{\phi}} \sum\_{m=1}^{3} b\_{n}^{m} \frac{\partial \nu\_{t}}{\partial \, \xi^{m}} \right) - f^{j}$$

$$S\_{\phi} = s\_{\phi} - 2\left(g^{12}\frac{\partial^2 \phi}{\partial \xi^1 \partial \xi^2} + g^{13}\frac{\partial^2 \phi}{\partial \xi^1 \partial \xi^3} + g^{23}\frac{\partial^2 \phi}{\partial \xi^2 \partial \xi^3}\right)$$

The geometric coefficients *<sup>j</sup> <sup>i</sup> <sup>b</sup>* , *ij <sup>g</sup>* , and *<sup>j</sup> <sup>f</sup>* appearing in the above equations are defined by Thompson et al. (1985). The transport equations (22) on a computational cell (shown in Fig.7(a)) are linearized and evaluating coefficients and source term at the center node P of the element yields

$$\sum\_{j=1}^{3} \left( \operatorname{g}\_{\text{P}}^{\text{ij}} \frac{\widehat{\operatorname{\mathcal{O}}}^{2} \phi}{\widehat{\operatorname{\mathcal{O}}} \xi^{j} \widehat{\operatorname{\mathcal{O}}}^{\text{j}}} - \operatorname{2} \left( a\_{\phi}^{j} \right)\_{\text{P}} \frac{\widehat{\operatorname{\mathcal{O}} \phi}}{\widehat{\operatorname{\mathcal{O}}} \xi^{j}} \right) = \left( R\_{\phi} \right)\_{\text{P}} \frac{\widehat{\operatorname{\mathcal{O}} \phi}}{\widehat{\operatorname{\mathcal{O}} \pi}} + \left( S\_{\phi} \right)\_{\text{P}} \tag{23}$$

or

100 Fluid Dynamics, Computational Modeling and Applications

the dimensionless distance to the wall. In the outer region, eddy viscosity is given by

*C*Kleb and *C*wk are additional model constants. Numerical values for the model constants are

In the following, discretization and velocity-pressure coupling of the present RANS method are described. First, it is convenient to rewrite the transport equations for momentum (*Ui*) in

> <sup>1</sup> *<sup>t</sup> <sup>j</sup> j j <sup>j</sup> R U s*

Then the continuity equation (16) and the transport equations (20) for momentum

3 3

 

3 3

1 1

 

*<sup>n</sup> <sup>j</sup> j j <sup>m</sup> <sup>t</sup> n n <sup>n</sup> <sup>m</sup> n m <sup>R</sup> <sup>x</sup> a bU b f J J*

<sup>1</sup> <sup>0</sup> *<sup>i</sup> j j i*

*g a RS* 

 

*b U*

1 1

2 *jj <sup>j</sup> jj j <sup>j</sup>*

<sup>1</sup> <sup>2</sup>

*i j*

*J*

 

source function for the corresponding quantity. We transform the physical space (*xi*

again represents any one of the convective transport quantities (*Ui*), and *s*

*x x t*

*i* ,

 

3

*t J*

(22)

1 1 *<sup>j</sup> i i <sup>j</sup> <sup>j</sup> e b J*

3 3

(21)

1 1 1 *<sup>i</sup> j <sup>i</sup> <sup>j</sup> i j <sup>x</sup> <sup>b</sup>*

(20)

t)outer=*KC*cp*F*wake*F*Kleb, where *K* and *C*cp are model constants, *F*wake=min(*y*max*F*max, *C*wk*y*max*U*dif2/*F*max), and *F*Kleb= [1+5.5(*C*Kleb*y*/*y*max)6]-1. The *F*max and *y*max are determined by the value and corresponding location, respectively, of the maximum of *F*=*y*|*ω*|[1-exp( *y*+/*A*+)]. The quantity *U*dif is the difference between maximum and minimum velocity

and *A*+ are model constants, |*ω*| is the vorticity magnitude, and *y*+ is

*y*[1-exp(-*y*+/*A*+)] is the turbulent length scale

max ( ) *UVW* - 2 2 2 1/2

) using the following coordinate

 

 

min ( ) *UVW* .

is the

,*t*) into a

Driest formula, i.e., (

for the inner region,

the following general form:

(

*A*+=26, 

where 

where

transformations:

t)inner=2|*ω*|, where =

magnitudes in the profile and is expressed as *U*dif= 2 2 2 1/2

3

> ,

*ij j <sup>i</sup> j j i j <sup>j</sup> g f* 

,

3 2

1

1

=0.4, *K*=0.0168, *C*cp=1.6, *C*wk=1.0, and *C*Kleb=0.3.

**4.2 Discretization and velocity-pressure coupling** 

2

rectangular region in the computational space (

2

parameters can be written as

, *i i <sup>j</sup> t xx*

3 3 2 3

1 1 1

 

$$\mathcal{g}\_{\rm P}^{11} \phi\_{\hat{\boldsymbol{\xi}}^{1} \hat{\boldsymbol{\xi}}^{1}} + \mathcal{g}\_{\rm P}^{22} \phi\_{\hat{\boldsymbol{\xi}}^{2} \hat{\boldsymbol{\xi}}^{2}} + \mathcal{g}\_{\rm P}^{33} \phi\_{\hat{\boldsymbol{\xi}}^{3} \hat{\boldsymbol{\xi}}^{3}} = \mathcal{2} \left( \mathbf{C}\_{\phi} \right)\_{\rm P} \phi\_{\hat{\boldsymbol{\xi}}^{1}} + \mathcal{2} \left( \mathbf{B}\_{\phi} \right)\_{\rm P} \phi\_{\hat{\boldsymbol{\xi}}^{2}} + \mathcal{2} \left( \mathbf{A}\_{\phi} \right)\_{\rm P} \phi\_{\hat{\boldsymbol{\xi}}^{3}} + \left( \mathbf{R}\_{\phi} \right)\_{\rm P} \phi\_{\boldsymbol{t}} + \left( \mathbf{S}\_{\phi} \right)\_{\rm P} \quad \text{(24)}$$

Fig. 7. Definition sketch of a computational cell (a), and nodes in regular grid (b) and continuity cell (c).

The dimensions of the computational cell are 2*l*×2*k*×2*h*, where *l*=1/ <sup>11</sup> *Pg* , *k*=1/ <sup>22</sup> *Pg* , and *h*=1/ <sup>33</sup> *Pg* . The above equation is discretized by the finite-analytic scheme. Solution dependent coefficients are analytically derived by solving the above linearized transport equation using a hybrid method which combines a two-dimensional analytic solution in plane with one dimensional analytic solution in the direction. By specifying boundary conditions on the faces of the cell as a combination of exponential and linear functions, which are the natural solutions for the linearized transport equation, Equation (24) can be solved by the method of separation of variables. When the solution is evaluated at the center node P of the element, the following twelve-point finite analytic formula is obtained:

$$\boldsymbol{\phi}\_{\text{P}}^{n} = \frac{1}{1 + \mathsf{C}\_{\text{P}} \left( \mathsf{C}\_{\text{U}} + \mathsf{C}\_{\text{D}} + \frac{\mathsf{R}\_{\text{NN}}}{\Delta \mathsf{r}} \right)} + \begin{bmatrix} \mathsf{C}\_{\text{NE}} \boldsymbol{\phi}\_{\text{NE}}^{n} + \mathsf{C}\_{\text{NN}} \boldsymbol{\phi}\_{\text{NN}}^{n} + \mathsf{C}\_{\text{SE}} \boldsymbol{\phi}\_{\text{SE}}^{n} + \mathsf{C}\_{\text{SN}} \boldsymbol{\phi}\_{\text{SN}}^{n} \\\\ + \mathsf{C}\_{\text{EC}} \boldsymbol{\phi}\_{\text{EC}}^{n} + \mathsf{C}\_{\text{WC}} \boldsymbol{\phi}\_{\text{NC}}^{n} + \mathsf{C}\_{\text{NC}} \boldsymbol{\phi}\_{\text{NC}}^{n} + \mathsf{C}\_{\text{SC}} \boldsymbol{\phi}\_{\text{NC}}^{n} \\\\ + \mathsf{C}\_{\text{P}} \left( \mathsf{C}\_{\text{U}} \boldsymbol{\phi}\_{\text{U}}^{n} + \mathsf{C}\_{\text{D}} \boldsymbol{\phi}\_{\text{D}}^{n} + \frac{\mathsf{R}\_{\phi}}{\Delta \mathsf{r}} \boldsymbol{\phi}\_{\text{P}}^{n-1} \right) - \mathsf{C}\_{\text{P}} \left( \boldsymbol{\S}\_{\phi} \right)\_{\text{P}} \end{bmatrix} \tag{25}$$

where

Sail Performance Analysis of Sailing Yachts by Numerical Calculations and Experiments 103

PU D

The coefficients and the modified pseudovelocities in the above equations are defined at the staggered node, and obtained from those at the regular node by the one-dimensional linear interpolation. The solution of the complete flow equations involves a global iteration process, in which the velocity-pressure coupling is effected by predictor-corrector steps. In the predictor step, the pressure field at the previous time step is used in the solution of the implicit equations (25) to obtain the corresponding velocity field. Since the velocity field generally does not satisfy mass conservation, a corrector step is needed. In the corrector step, the explicit momentum equations (28) and the implicit pressure equation (26) are

**4.3 Multiblock (domain decomposition) capability, and overall numerical solution** 

al., 2009, for the results for which mast influences in computation are considered).

Fig. 8. Overview of the present multiblock computational grid.

As mentioned earlier, the multiblock (domain decomposition) capability is facilitated in the present RANS method. This capability is essential for simulation of flow around complex geometry, e.g., multiple sail system for sailing yacht as focused in the present study. Fig.8 shows overview of the present multiblock computational grid, while the grid is generated by using an automatic gridding scheme developed by the present author (Masuyama et al., 2009) Note that the gridding engine together with the present RANS method was recently implemented into a comprehensive sail performance prediction software "*Advanced Aero Flow*" (Katori, 2009). See the reference for more details of the scheme. Total number of grids is around a half million, and the number of multiblock is 48. Free-stream, symmetry, and wall-surface (no slip) boundary conditions are imposed on outer and top boundaries, bottom boundary, and sail surface boundary, respectively. For the results shown in this paper, the mast and rigging are not considered in the series calculations, and the bottom boundary is located at the same height as that of deck plane of the boat (see Masuyama et

*i i <sup>i</sup> <sup>j</sup> <sup>j</sup> R C <sup>p</sup> U U <sup>b</sup> <sup>R</sup> J CC C*

ˆ

1

solved iteratively to ensure the satisfaction of the continuity equation.

**procedure** 

P

3

1

*j*

(28)

$$\mathbf{C}\_{\rm U} = \frac{\mathbf{C}e^{\mathbf{C}l}}{I \sinh \mathbf{C}l}, \quad \mathbf{C}\_{\rm D} = \frac{\mathbf{C}e^{-\mathbf{C}l}}{I \sinh \mathbf{C}l}, \; \mathbf{C}\_{\rm SC} = \frac{e^{\mathbf{R}k}}{2 \cosh \mathbf{B}k} P\_{A}, \; \mathbf{C}\_{\rm WC} = e^{-2\mathbf{R}l} \mathbf{C}\_{\rm SC}$$

$$\mathbf{C}\_{\rm WC} = \frac{e^{4h}}{2 \cosh \mathbf{A}h} P\_{\mathbf{B}}, \; \mathbf{C}\_{\rm EC} = e^{-2Ah} \mathbf{C}\_{\rm WC}, \; \mathbf{C}\_{\rm SW} = \frac{e^{4h + 2k}}{4 \cosh \mathbf{A}h \cosh \mathbf{B}k} (1 - P\_A - P\_B)$$

$$\mathbf{C}\_{\rm SE} = e^{-2Ah} \mathbf{C}\_{\rm SW}, \; \mathbf{C}\_{\rm SW} = e^{-2Bk} \mathbf{C}\_{\rm SW}, \; \mathbf{C}\_{\rm NE} = e^{-2Ah - 2Bk} \mathbf{C}\_{\rm SW}$$

$$\mathbf{C}\_{\rm P} = \frac{h \tanh \mathbf{A}h}{2A} \left( 1 - P\_A \right) = \frac{k \tanh \mathbf{B}k}{2B} \left( 1 - P\_B \right), \; P\_A = 4E\_2 Ah \cosh \mathbf{A}h \cosh \mathbf{B}k \coth \mathbf{A}h$$

$$P\_B = 1 + \frac{Bh \coth \mathbf{B}k}{4A \cosh \mathbf{A}h} (P\_A - 1), \; E\_B = \sum\_{i=1}^{\infty} \frac{-(-1)^{\text{int}} \left( \lambda\_m^i h \right)}{\lambda\_m^{-2} \lambda\_i^2}$$

coth *B A Ak Ah* , <sup>2</sup> <sup>2</sup> <sup>2</sup> <sup>2</sup> 22 2 <sup>1</sup> cos *<sup>m</sup> m m Ah h A B k* 1 2 *mh m* 

The subscripts P, U and D denote the center, upstream and downstream nodes, respectively, and NC, NW, WC, etc. denote the nodes in the -plane in terms of compass directions. The superscripts (n) and (n-1) refer to the current and previous time levels, and is the time step. The solution of the complete flow equations involves a global iteration process, in which the velocity-pressure coupling is effected by PISO-type predictor-corrector steps. The pressure equation is derived by introducing pseudo-velocities at staggered locations while maintaining the regular grid arrangement for all the transport equations. Fig.7(b) and (c) show the locations of nodes in the regular grid in the -plane. All transport quantities and pressure are evaluated at the regular nodes. In deriving the pressure equation, a control volume is employed as a continuity cell, to establish the coupling between the velocity and pressure fields. The pressure equation used in this study is written as

$$\begin{aligned} \left(E\_{\rm d}^{11} + E\_{\rm u}^{11} + E\_{\rm n}^{22} + E\_{\rm s}^{22} + E\_{\rm e}^{33} + E\_{\rm w}^{33}\right) p\_{\rm P} \\ = E\_{\rm d}^{11} p\_{\rm D} + E\_{\rm u}^{11} p\_{\rm U} + E\_{\rm n}^{22} p\_{\rm NC} + E\_{\rm s}^{22} p\_{\rm SC} + E\_{\rm e}^{33} p\_{\rm EC} + E\_{\rm w}^{33} p\_{\rm WC} - \hat{D} \end{aligned} \tag{26}$$

with

$$
\hat{D} = \hat{\mathcal{U}}\_{\rm d}^{1} - \hat{\mathcal{U}}\_{\rm u}^{1} + \hat{\mathcal{U}}\_{\rm n}^{2} - \hat{\mathcal{U}}\_{\rm s}^{2} + \hat{\mathcal{U}}\_{\rm e}^{3} - \hat{\mathcal{U}}\_{\rm w}^{3}
$$

Here *ij E* and a modified pseudovelocity ˆ *<sup>i</sup> U* at the regular node are

$$E^{ij} = \frac{R \cdot \mathbb{C}\_{\text{P}}}{J \left[1 + \mathbb{C}\_{\text{P}} \left(\mathbb{C}\_{\text{U}} + \mathbb{C}\_{\text{D}} + \frac{R}{\Delta \tau} \right) \right]} \sum\_{m=1}^{3} b\_{m}^{i} b\_{m}^{j} \; , \; \hat{\mathcal{U}}^{i} = \sum\_{n=1}^{3} b\_{n}^{i} \hat{\mathcal{U}}\_{n} - E^{ij} \frac{\partial p}{\partial \xi^{j}} - E^{ik} \frac{\partial p}{\partial \xi^{k}} \tag{27}$$

where ˆ*Ui* is a pseudovelocity given by the decomposition of Equation (25) for *Ui* into ˆ*Ui* plus the pressure gradient terms, such that

*Bk*

*Ah C eC* , SW <sup>1</sup>

*<sup>e</sup> C P*

*Bk C eC* , 2 2

3

1

ˆ ˆ *i i ij ik n n <sup>j</sup> <sup>k</sup> <sup>n</sup> <sup>p</sup> <sup>p</sup> U bU E E*

(27)

1 2

 

 

The subscripts P, U and D denote the center, upstream and downstream nodes, respectively,

time step. The solution of the complete flow equations involves a global iteration process, in which the velocity-pressure coupling is effected by PISO-type predictor-corrector steps. The pressure equation is derived by introducing pseudo-velocities at staggered locations while maintaining the regular grid arrangement for all the transport equations. Fig.7(b) and (c)

pressure are evaluated at the regular nodes. In deriving the pressure equation, a control volume is employed as a continuity cell, to establish the coupling between the velocity and

> 11 11 22 22 33 33 d D u U n NC s SC e EC w WC <sup>ˆ</sup>

> > 3

*m*

*Ep Ep Ep Ep Ep Ep D*

112233 duns ew *DU U U U U U* ˆˆ ˆ ˆ ˆ ˆ ˆ

*m m*

,

where ˆ*Ui* is a pseudovelocity given by the decomposition of Equation (25) for *Ui* into ˆ*Ui*

*A*

NC SC *Bk C eC*

*A B*

 

> is the

(26)

*Bk* , <sup>2</sup>

4 cosh cosh *Ah Bk*

<sup>2</sup> <sup>2</sup> <sup>2</sup> <sup>2</sup> 22 2 <sup>1</sup> 1

*<sup>m</sup> m m*

cos

*Ah h A B k*

*h*



*m m*

*<sup>e</sup> <sup>C</sup> P P Ah Bk*

> NE SW *Ah Bk Ce C*

 

, SC 2 cosh

NW SW

*A B* , 2 *P E Ah Ah Bk Ah <sup>A</sup>* 4 cosh cosh coth

*Ak Ah* ,

*E*

*mh m*

The superscripts (n) and (n-1) refer to the current and previous time levels, and

<sup>U</sup> sinh *Cl Ce <sup>C</sup>*

*Ah*

<sup>P</sup> tanh tanh 1 1 2 2 *<sup>A</sup> <sup>B</sup>*

> coth 1 1 coth *B A Bh Bk P P*

and NC, NW, WC, etc. denote the nodes in the

show the locations of nodes in the regular grid in the

with

pressure fields. The pressure equation used in this study is written as

 11 11 22 22 33 33 d u n s e wP

*E E E E E Ep*

Here *ij E* and a modified pseudovelocity ˆ *<sup>i</sup> U* at the regular node are

P <sup>1</sup> 1 PU D *ij j i*

*R C <sup>E</sup> b b <sup>R</sup> J CC C*

plus the pressure gradient terms, such that

*h Ah k Bk CPP*

WC 2 cosh

*<sup>e</sup> C P*

*l Cl* , D sinh

*B*

2 SE SW

*Ah C eC* , <sup>2</sup>

*Ah* , <sup>2</sup>

*Cl Ce <sup>C</sup> l Cl* 

EC WC

$$\hat{\mathcal{U}}I\_i = \hat{\mathcal{U}}\_i - \frac{R \, \mathcal{C}\_\text{P}}{J \left[ 1 + \mathcal{C}\_\text{P} \left( \mathcal{C}\_\text{U} + \mathcal{C}\_\text{D} + \frac{R}{\Delta \tau} \right) \right]} \sum\_{j=1}^3 b\_i^j \frac{\partial p}{\partial \xi^j} \tag{28}$$

The coefficients and the modified pseudovelocities in the above equations are defined at the staggered node, and obtained from those at the regular node by the one-dimensional linear interpolation. The solution of the complete flow equations involves a global iteration process, in which the velocity-pressure coupling is effected by predictor-corrector steps. In the predictor step, the pressure field at the previous time step is used in the solution of the implicit equations (25) to obtain the corresponding velocity field. Since the velocity field generally does not satisfy mass conservation, a corrector step is needed. In the corrector step, the explicit momentum equations (28) and the implicit pressure equation (26) are solved iteratively to ensure the satisfaction of the continuity equation.

#### **4.3 Multiblock (domain decomposition) capability, and overall numerical solution procedure**

As mentioned earlier, the multiblock (domain decomposition) capability is facilitated in the present RANS method. This capability is essential for simulation of flow around complex geometry, e.g., multiple sail system for sailing yacht as focused in the present study. Fig.8 shows overview of the present multiblock computational grid, while the grid is generated by using an automatic gridding scheme developed by the present author (Masuyama et al., 2009) Note that the gridding engine together with the present RANS method was recently implemented into a comprehensive sail performance prediction software "*Advanced Aero Flow*" (Katori, 2009). See the reference for more details of the scheme. Total number of grids is around a half million, and the number of multiblock is 48. Free-stream, symmetry, and wall-surface (no slip) boundary conditions are imposed on outer and top boundaries, bottom boundary, and sail surface boundary, respectively. For the results shown in this paper, the mast and rigging are not considered in the series calculations, and the bottom boundary is located at the same height as that of deck plane of the boat (see Masuyama et al., 2009, for the results for which mast influences in computation are considered).

Fig. 8. Overview of the present multiblock computational grid.

Sail Performance Analysis of Sailing Yachts by Numerical Calculations and Experiments 105

Step 12. Return to Step 4 for the next time step, until the time step reaches the given

More details of the present RANS method are described in Tahara et al. (2006a, 2006b) in

**5. Measurements of upwind sail performance in full-scale condition using sail** 

Full-scale onboard measurements are free from scale-effect problem by wind tunnel tests and appear more promising, but the challenge becomes how to accurately measure forces acting on the sail. Such studies on sail force measurements were performed by Milgram et al. (1993), Masuyama et al. (1997a, 1997b), and Hochkirch et al. (1999), who built full-scale

Milgram (1993) showed in his pioneering work that the sail dynamometer boat, *Amphetrete*, is quite capable. This measurement system consists of a 35-foot boat with an internal frame connected to the hull by six load cells, which were configured to measure all forces and moments acting on the sails. In his work, the sail shapes were also measured and used for CFD analyses; however unfortunately, details of the sail shape and performance data were not presented. Hochkirch et al. (1999) also built a 33-foot dynamometer boat *DYNA*. The aerodynamic forces acting on the sail were measured and compared with the results from wind tunnel tests (Hansen et al. 2003). The measured data were also used as input to the CFD calculation and a parametric survey was carried out (Krebber et al. 2006). Masuyama and Fukasawa were encouraged by Milgram's work, and built a sail dynamometer boat, *Fujin*. The measurement system installed in the *Fujin* and the results of calibration test and

The *Fujin* was originally built for conducting tests on sails for the Japanese America's Cup entry in 1994. *Fujin* is a 10.3m-long ocean cruiser with a sail dynamometer system in the hull which can directly measure sail forces and moments. Fig. 1 shows the general arrangement of the *Fujin*. The test sails were made to correspond to a typical sail plan for an International Measurement System (IMS) class boat. The rigging of the *Fujin* was originally designed for testing sails for the International America's Cup Class (IACC) boat. The jib of IACC boat is relatively small. Therefore, the longitudinal position of the jib rail track of the *Fujin* was located further forward than that of the typical IMS boat. For this reason, the tests were performed using a fully batten mainsail and a 130% jib instead of a 150% jib. The sails were made by North Sails Japan. The axes system is also shown in Fig. 1. The origin is located on the vessel's centerline at the aft face of the mast (*x*-direction), and the height of deck level at the base of the forestay (*z*-direction). Table 1 shows the principal dimensions of the boat and the detailed measurements of the sails, where "I", "J", "P" and "E" are the measurement

The aerodynamic coefficients and the coordinates of the center of effort of the sails are

Step 11. Repeat Step 8 through Step 10 for the specified number of times.

maximum value.

**dynamometer boat** *Fujin* **5.1 Full-scale measurements** 

addition to the above-cited references.

boats with onboard sail dynamometer systems.

sailing test were reported by Masuyama et al. (1997a and 1997b).

lengths of sail dimensions for the IMS rule as defined in Fig.1.

defined as follows:

**5.2 Measurements by sail dynamometer boat** *Fujin* 

The basic strategy to handle the multiblock follows domain decomposition technique to solve the elliptic PDE by using several subdomains. After adequate discretization is applied and a simple preconditioner is introduced, the discrete alternating Schwarz's method to solve the PDE is used for boundary matching. Finally, Fig.9 shows the code structure of the present RANS method, and the overall numerical solution procedure of the present RANS method is summarized as follows:

Fig. 9. PISO type solution algorithm for the present multiblock RANS method.


Step 11. Repeat Step 8 through Step 10 for the specified number of times.

Step 12. Return to Step 4 for the next time step, until the time step reaches the given maximum value.

More details of the present RANS method are described in Tahara et al. (2006a, 2006b) in addition to the above-cited references.

#### **5. Measurements of upwind sail performance in full-scale condition using sail dynamometer boat** *Fujin*

#### **5.1 Full-scale measurements**

104 Fluid Dynamics, Computational Modeling and Applications

The basic strategy to handle the multiblock follows domain decomposition technique to solve the elliptic PDE by using several subdomains. After adequate discretization is applied and a simple preconditioner is introduced, the discrete alternating Schwarz's method to solve the PDE is used for boundary matching. Finally, Fig.9 shows the code structure of the present RANS method, and the overall numerical solution procedure of the present RANS

Receive Grid

,

Boundary Matching

Predictor Compute U,V,W

Corrector Compute P Correct U,V,W

> Boundary Matching

> > Inter Block Communications

 

 

*m n m n m m m m*

*m n m n m m m m*

*A u*

I

*A u*

I

 

*n n m m*

 

*n n m m*

 <sup>1</sup> <sup>1</sup> \ <sup>1</sup> <sup>1</sup> *<sup>n</sup>*

*<sup>f</sup> <sup>A</sup> <sup>g</sup> <sup>A</sup> <sup>u</sup> <sup>u</sup> <sup>u</sup> <sup>A</sup>*

*m m m*

*m m m*

 <sup>1</sup> <sup>1</sup> \ <sup>1</sup> <sup>1</sup> *<sup>n</sup>*

*<sup>f</sup> <sup>A</sup> <sup>g</sup> <sup>A</sup> <sup>u</sup> <sup>u</sup> <sup>u</sup> <sup>A</sup>*

 

 

*n*

*n*

method is summarized as follows:

Receive Grid

, Compute <sup>t</sup>

Boundary Matching

Predictor Compute U,V,W

Corrector Compute P Correct U,V,W

> Boundary Matching

Conv.? Conv.? Conv.? No No

Yes Yes Yes

Send sols. Send sols.

Fig. 9. PISO type solution algorithm for the present multiblock RANS method.

Step 4. Compute the finite-analytic coefficients for the transport equation.

Step 1. Input the computational grid, setup parameters, and boundary condition

Step 6. Solve transport equation for velocities (*U*,*V*,*W*) using the previous pressure field

Step 9. Using the newly obtained pressure, calculate the new velocity field explicitly

Step 10. Update the finite-analytic coefficients for the transport equation for velocities

Step 2. Specify the initial conditions for the velocity, pressure and turbulence fields.

(Block 2) .... (Block m)

, Compute <sup>t</sup>

START (Block 1)

Read Grid Send to Block. m

Predictor Compute U,V,W Boundary Matching

Compute <sup>t</sup>

Corrector Compute P Correct U,V,W

> Boundary Matching

Receive sols. from Bk. m Write sols.

STOP

information.

(*U*,*V*,*W*).

Step 3. Compute the geometric coefficients.

Step 5. Compute eddy viscosity distribution.

(predictor stage for the velocity field).

(corrector stage for the velocity field).

Step 8. Solve pressure equation.

Step 7. Compute the coefficients of pressure equation.

No

Full-scale onboard measurements are free from scale-effect problem by wind tunnel tests and appear more promising, but the challenge becomes how to accurately measure forces acting on the sail. Such studies on sail force measurements were performed by Milgram et al. (1993), Masuyama et al. (1997a, 1997b), and Hochkirch et al. (1999), who built full-scale boats with onboard sail dynamometer systems.

Milgram (1993) showed in his pioneering work that the sail dynamometer boat, *Amphetrete*, is quite capable. This measurement system consists of a 35-foot boat with an internal frame connected to the hull by six load cells, which were configured to measure all forces and moments acting on the sails. In his work, the sail shapes were also measured and used for CFD analyses; however unfortunately, details of the sail shape and performance data were not presented. Hochkirch et al. (1999) also built a 33-foot dynamometer boat *DYNA*. The aerodynamic forces acting on the sail were measured and compared with the results from wind tunnel tests (Hansen et al. 2003). The measured data were also used as input to the CFD calculation and a parametric survey was carried out (Krebber et al. 2006). Masuyama and Fukasawa were encouraged by Milgram's work, and built a sail dynamometer boat, *Fujin*. The measurement system installed in the *Fujin* and the results of calibration test and sailing test were reported by Masuyama et al. (1997a and 1997b).

#### **5.2 Measurements by sail dynamometer boat** *Fujin*

The *Fujin* was originally built for conducting tests on sails for the Japanese America's Cup entry in 1994. *Fujin* is a 10.3m-long ocean cruiser with a sail dynamometer system in the hull which can directly measure sail forces and moments. Fig. 1 shows the general arrangement of the *Fujin*. The test sails were made to correspond to a typical sail plan for an International Measurement System (IMS) class boat. The rigging of the *Fujin* was originally designed for testing sails for the International America's Cup Class (IACC) boat. The jib of IACC boat is relatively small. Therefore, the longitudinal position of the jib rail track of the *Fujin* was located further forward than that of the typical IMS boat. For this reason, the tests were performed using a fully batten mainsail and a 130% jib instead of a 150% jib. The sails were made by North Sails Japan. The axes system is also shown in Fig. 1. The origin is located on the vessel's centerline at the aft face of the mast (*x*-direction), and the height of deck level at the base of the forestay (*z*-direction). Table 1 shows the principal dimensions of the boat and the detailed measurements of the sails, where "I", "J", "P" and "E" are the measurement lengths of sail dimensions for the IMS rule as defined in Fig.1.

The aerodynamic coefficients and the coordinates of the center of effort of the sails are defined as follows:

Sail Performance Analysis of Sailing Yachts by Numerical Calculations and Experiments 107

as an Optical Fiber Gyroscope (roll and pitch angles), a Flux Gate Compass (heading angle), a Differential type GPS receiver, a speedometer (velocity in the x direction) and a potentiometer for rudder angle. These data are recorded by an onboard computer

simultaneously with the data from the load cells.

Fig. 10. General arrangement of dynamometer frame in *Fujin*.

Fig. 11. Sea trial condition in light wind with 130% jib.

$$\mathbf{C}\_{X} = \frac{X\_{S}}{\frac{1}{2}\rho\_{a}\mathcal{U}\_{A}\mathcal{^2S}\_{A}}, \qquad \mathbf{C}\_{Y} = \frac{Y\_{S}}{\frac{1}{2}\rho\_{a}\mathcal{U}\_{A}\mathcal{^2S}\_{A}}, \qquad \mathbf{x}\_{\rm CE} = \frac{N\_{S}}{Y\_{S}}, \qquad \mathbf{z}\_{\rm CE} = \frac{K\_{S}}{Y\_{S}} \tag{29}$$

where *XS* and *YS* are the force components along the *x* and *y* axes of the boat respectively, and *KS* and *NS* are the moments around the *x* and *z* axes. *xCE* and *zCE* are the *x* and *z* coordinates of the center of effort of the sails (CE). The thrust force coefficient *CX* is expressed as positive for the forward direction and the side force coefficient *CY* is positive for both port and starboard directions. It should be noted that the coordinates are given in the body axes system. Therefore, when the boat heels the *YS* force component is not in the horizontal plane but is normal to the mast. The aerodynamic forces acting on the mast and rigging are included in the measured sail forces.

#### **5.2.1 Measurement system of aerodynamic performance and sail shape**

The sail dynamometer system is composed of a rigid aluminum frame and four load cells. The frame is separated structurally from the hull and connected to it by the load cells. The general arrangement of the dynamometer frame is given in Fig.10. The load cells are numbered in the figure. Two of these are 1-component load cells and the others are 2 component ones. Hence, these load cells form a 6-component dynamometer system, and their outputs can be transformed to the forces and moments about the boat axes using a calibration matrix. All rig components such as mast, chain plates, winches, lead blocks, etc. are attached to the aluminum frame. The under deck portion of the mast is held by the frame, and the other rig components are attached to the frame through the deck holes. The data acquisition system and calibration method for the *Fujin* were described by Masuyama et al. (1997a and 1997b).

The sail shape was recorded using pairs of CCD cameras. The lower part of the mainsail was photographed using the CCD camera pair designated A in Fig.11. These were located at the mast top, 50 cm transversely from each side of the mast. The upper part of the mainsail was photographed using a portable video camera from below the boom. The lower part of the jib was photographed using the camera pair designated B in Fig.11, which were located at the intersection point of the forestay and the mast, 10 cm transversely from each side of the mast. The upper part of the jib was photographed using a portable video camera from inside the bow hatch. For measuring convenience, horizontal stripes were drawn on the mainsail and jib at heights of 10, 20, 40, 60 and 80% of each sail. The sail shape images were analyzed using the sail shape analyzing software, SSA-2D, developed by Armonicos Co. Fig.12(a) shows an example of processed image of the mainsail using the SSA-2D. This software calculates the curvature of the sail section by marking several points of the sail stripe and the reference line on the PC display, and indicates the parameters such as chord length, maximum draft, maximum draft position, entry angle at the luff, i.e., leading edge, and exit angle at the leech, i.e., trailing edge, as shown in Fig.12(b). The apparent wind speed (AWS) and apparent wind angle (AWA) are measured by an anemometer attached on the "Bow unit" as shown in Fig.11. This unit post can rotate freely to maintain its vertical attitude when the boat heels in order to measure the wind data in the horizontal plane. The height of the anemometer coincides with the geometric center of effort (GCE) of the sail plan. The wind speed and wind angle sensors were calibrated by wind tunnel tests in advance and the calibration equations were obtained. The *Fujin* also has motion measuring instruments such

*X Y NK S S SS <sup>C</sup> , C ,x , z <sup>X</sup> <sup>Y</sup> CE Y Y CE*

where *XS* and *YS* are the force components along the *x* and *y* axes of the boat respectively, and *KS* and *NS* are the moments around the *x* and *z* axes. *xCE* and *zCE* are the *x* and *z* coordinates of the center of effort of the sails (CE). The thrust force coefficient *CX* is expressed as positive for the forward direction and the side force coefficient *CY* is positive for both port and starboard directions. It should be noted that the coordinates are given in the body axes system. Therefore, when the boat heels the *YS* force component is not in the horizontal plane but is normal to the mast. The aerodynamic forces acting on the mast and

The sail dynamometer system is composed of a rigid aluminum frame and four load cells. The frame is separated structurally from the hull and connected to it by the load cells. The general arrangement of the dynamometer frame is given in Fig.10. The load cells are numbered in the figure. Two of these are 1-component load cells and the others are 2 component ones. Hence, these load cells form a 6-component dynamometer system, and their outputs can be transformed to the forces and moments about the boat axes using a calibration matrix. All rig components such as mast, chain plates, winches, lead blocks, etc. are attached to the aluminum frame. The under deck portion of the mast is held by the frame, and the other rig components are attached to the frame through the deck holes. The data acquisition system and calibration method for the *Fujin* were described by Masuyama

The sail shape was recorded using pairs of CCD cameras. The lower part of the mainsail was photographed using the CCD camera pair designated A in Fig.11. These were located at the mast top, 50 cm transversely from each side of the mast. The upper part of the mainsail was photographed using a portable video camera from below the boom. The lower part of the jib was photographed using the camera pair designated B in Fig.11, which were located at the intersection point of the forestay and the mast, 10 cm transversely from each side of the mast. The upper part of the jib was photographed using a portable video camera from inside the bow hatch. For measuring convenience, horizontal stripes were drawn on the mainsail and jib at heights of 10, 20, 40, 60 and 80% of each sail. The sail shape images were analyzed using the sail shape analyzing software, SSA-2D, developed by Armonicos Co. Fig.12(a) shows an example of processed image of the mainsail using the SSA-2D. This software calculates the curvature of the sail section by marking several points of the sail stripe and the reference line on the PC display, and indicates the parameters such as chord length, maximum draft, maximum draft position, entry angle at the luff, i.e., leading edge, and exit angle at the leech, i.e., trailing edge, as shown in Fig.12(b). The apparent wind speed (AWS) and apparent wind angle (AWA) are measured by an anemometer attached on the "Bow unit" as shown in Fig.11. This unit post can rotate freely to maintain its vertical attitude when the boat heels in order to measure the wind data in the horizontal plane. The height of the anemometer coincides with the geometric center of effort (GCE) of the sail plan. The wind speed and wind angle sensors were calibrated by wind tunnel tests in advance and the calibration equations were obtained. The *Fujin* also has motion measuring instruments such

**5.2.1 Measurement system of aerodynamic performance and sail shape** 

*<sup>ρ</sup> U S <sup>ρ</sup> U S S S aA A aA A* (29)

1 1 2 2

2 2

rigging are included in the measured sail forces.

et al. (1997a and 1997b).

as an Optical Fiber Gyroscope (roll and pitch angles), a Flux Gate Compass (heading angle), a Differential type GPS receiver, a speedometer (velocity in the x direction) and a potentiometer for rudder angle. These data are recorded by an onboard computer simultaneously with the data from the load cells.

Fig. 10. General arrangement of dynamometer frame in *Fujin*.

Fig. 11. Sea trial condition in light wind with 130% jib.

Sail Performance Analysis of Sailing Yachts by Numerical Calculations and Experiments 109

Fig.13 shows the performance variation for the mainsail and 130% jib configuration as a function of AWA. In the figure the solid symbols indicate the experimental results and the open symbols indicate the calculated results using the VLM and the RANS-based CFD. For the experimental results, both data from the starboard (Stbd) and port tack (Port) are shown. All of the measured coefficients are plotted with error bars indicating the range of deviation over the averaging period. There are some discrepancies between the data from each tack. During the experiments, efforts were made to remove this asymmetrical performance. However, the boat speed actually differed on each tack. It can be concluded that there was a slight asymmetry in the combination of the hull, keel, rudder and dynamometer frame. The numerical calculations were performed using the measured shape data. In order to avoid confusion when interpreting the figure, the calculated results are indicated only for the port tack. Therefore, the calculated and experimental

Fig. 13. Performance variation as a function of apparent wind angle (AWA) for mainsail

**5.3.1 Variation with apparent wind angle** 

points for the port tack correspond to each other.

and 130% jib.

Fig. 12. (a) Example of a processed image of the mainsail using SSA-2D. (b) Measured sail shape parameters.

#### **5.2.2 Test condition and error analysis**

The sea tests were performed on Nanao Bay off the Noto Peninsula. The bay is approximately eight nautical miles from east to west and five from north to south. The bay is surrounded by low hills, and the mouth connecting it to the Japan Sea is narrow. Therefore, there is little tidal current in the bay, and the wave heights are relatively low even though the wind can be strong. The close-hauled tests were conducted over the apparent wind angle (AWA) range of 20 to 40 degrees, and the apparent wind speed (AWS) range of 5 to 11m/s. The effect of the AWA, and the draft and twist of the mainsail on the sail performance were measured. Data sampling was started when the sailing condition was considered to be in steady state. The sampling rate for the data acquisition system was set at 10Hz. Data sampling was continued for 90 seconds, and during this time the sail shapes were recorded using the CCD cameras. The boat was steered carefully during this time. However, the measured data contained some variation due to wind fluctuation and wave reflection on the hull. Therefore the steady state values for the aerodynamic coefficients were obtained by averaging the data over a 30 to 60 seconds period, in which the AWA was closer to the target value than during the whole 90 second period. For these tests if the range of deviation of AWA exceeded ±5 degrees, the results were discarded. All of the measured coefficients are plotted with error bars indicating the range of deviation over the averaging period.

#### **5.3 Comparison between experimental and calculated results**

In this chapter, the experimental results and the calculated results for the following cases will be compared:


For each series, first the sail coefficients: *CL*, *CD*, *CX* and *CY*, and the coordinates of *xCE* and *zCE* are given. Then, the calculated the sail surface pressure and streamlines using the RANSbased CFD are presented for two typical cases in each series. Finally, the shapes and threedimensional coordinates of the sails are tabulated for each case corresponding to those where the RANS-based CFD results are given.

#### **5.3.1 Variation with apparent wind angle**

108 Fluid Dynamics, Computational Modeling and Applications

Maximum draft position

(a) (b)

The sea tests were performed on Nanao Bay off the Noto Peninsula. The bay is approximately eight nautical miles from east to west and five from north to south. The bay is surrounded by low hills, and the mouth connecting it to the Japan Sea is narrow. Therefore, there is little tidal current in the bay, and the wave heights are relatively low even though the wind can be strong. The close-hauled tests were conducted over the apparent wind angle (AWA) range of 20 to 40 degrees, and the apparent wind speed (AWS) range of 5 to 11m/s. The effect of the AWA, and the draft and twist of the mainsail on the sail performance were measured. Data sampling was started when the sailing condition was considered to be in steady state. The sampling rate for the data acquisition system was set at 10Hz. Data sampling was continued for 90 seconds, and during this time the sail shapes were recorded using the CCD cameras. The boat was steered carefully during this time. However, the measured data contained some variation due to wind fluctuation and wave reflection on the hull. Therefore the steady state values for the aerodynamic coefficients were obtained by averaging the data over a 30 to 60 seconds period, in which the AWA was closer to the target value than during the whole 90 second period. For these tests if the range of deviation of AWA exceeded ±5 degrees, the results were discarded. All of the measured coefficients are plotted with error bars indicating the

In this chapter, the experimental results and the calculated results for the following cases

For each series, first the sail coefficients: *CL*, *CD*, *CX* and *CY*, and the coordinates of *xCE* and *zCE* are given. Then, the calculated the sail surface pressure and streamlines using the RANSbased CFD are presented for two typical cases in each series. Finally, the shapes and threedimensional coordinates of the sails are tabulated for each case corresponding to those

Fig. 12. (a) Example of a processed image of the mainsail using SSA-2D. (b) Measured sail

shape parameters.

will be compared:

**5.2.2 Test condition and error analysis** 

range of deviation over the averaging period.

a. Variation with apparent wind angle b. Variation with mainsail twist angle

where the RANS-based CFD results are given.

**5.3 Comparison between experimental and calculated results** 

Entry angle

Maximum draft

Chord length

> Exit angle

Sail section

Leech

Twist angle Center line Luff Fig.13 shows the performance variation for the mainsail and 130% jib configuration as a function of AWA. In the figure the solid symbols indicate the experimental results and the open symbols indicate the calculated results using the VLM and the RANS-based CFD. For the experimental results, both data from the starboard (Stbd) and port tack (Port) are shown. All of the measured coefficients are plotted with error bars indicating the range of deviation over the averaging period. There are some discrepancies between the data from each tack. During the experiments, efforts were made to remove this asymmetrical performance. However, the boat speed actually differed on each tack. It can be concluded that there was a slight asymmetry in the combination of the hull, keel, rudder and dynamometer frame. The numerical calculations were performed using the measured shape data. In order to avoid confusion when interpreting the figure, the calculated results are indicated only for the port tack. Therefore, the calculated and experimental points for the port tack correspond to each other.

Fig. 13. Performance variation as a function of apparent wind angle (AWA) for mainsail and 130% jib.

Sail Performance Analysis of Sailing Yachts by Numerical Calculations and Experiments 111

The calculated results for *CL* using the VLM show good agreement with the experiments at AWA angles less than about 35 degrees. Over about 35 degrees, the calculated results are lower than the measured ones. This shows that the calculated results strongly indicate the effect of incorrect sail trimming. The results for *CL* using the RANS-based CFD show the same trends with the experiments, but are slight higher than those from the experiments for AWA between 20 degrees to 30 degrees and lower for AWA greater than 30 degrees. In particular, the decrease in *CL* for AWA values greater than 30 degrees is considerably large. This will be discussed later with the calculated sail surface pressure and streamlines. The calculated results for *CD* slightly over predict those from the experiments. Fig.13(c) shows the coordinates of the center of effort of the sails. The *x* and *z* coordinates of the geometric center of effort (*xGCE* and *zGCE*) are 0.63m aft and 4.80m above the origin, which are indicated by alternate long and short dashed lines in the figure. It is seen that both the experimental and the calculated coordinates of *xCE* are near *xGCE* and move slightly forward with increasing AWA. Unfortunately, there is a wide scatter in the experimental values of *zCE*. This is thought to be because the measured *Ks* moment contains a large component from the mass of the dynamometer frame and rigging (659kg). This moment was subtracted from the measurement, taking into account the measured heel angle. If there is a slight error in the position of center of gravity of the dynamometer frame, or in the measured heel angle, the error in the calculated moment will be large. However, though there is a scatter in the measured data, it can be seen that *zCE* is decreasing as AWA increases. The trends in the movement of both *xCE* and *zCE* as functions of AWA might be caused by the decrement of force acting on the aft and upper parts of the sails due to the loosening of main and jib sheets with increasing AWA. The calculated results for *zCE* obtained using theRANS-based CFD show the same trend as the experiments. On the other hand, the calculated results using VLM are considerably higher than the experimental ones. This might be caused by over estimation of the force acting on the upper portion of the mainsail. In this area, since the jib is not overlapping, flow separation may occur easily. However, the VLM does not

Figures 14(1) and 14(2) show the calculated results of the sail surface pressure and streamlines using RANS- based CFD. Fig.14(1) indicates the case of experiment ID 96092335 (AWA= 30.7deg.), and 14(2) indicates ID 96080248 (AWA= 37.9deg.). These data correspond to the plotted points on the vertical dotted lines (1) and (2) in Fig.13. In Fig.14, the left and right diagrams correspond to the port and starboard sides, i.e., pressure and suction sides, respectively. In 14(1), although slight flow separation on the suction side of mainsail is seen, the streamlines of both sides run smoothly. On the other hand, in 14(2), considerable flow separation is occurring, in particular, on the suction side of jib. This is the main reason for the reduction of *CL* value in the RANS-based CFD calculation at (2) in Fig.13(a). The shapes and three-dimensional coordinates of the sails are given in Table 2. The numbered (1) and (2) tables correspond to the cases of experiment ID 96092335 and ID 96080248, respectively. These also correspond to the calculated results shown in Fig.14. The figures described above the tables show the sail section profiles at 0, 20, 40, 60 and 80% of the sail height. The dimensions of these three-dimensional coordinates are given in the tables including 100% height section data. The positive direction of the *x* coordinate is aft. The four lines at the top of the tables are the measured values for the wind and sail trim conditions, the boat attitude

take flow separation into account.

and the sail performance coefficients.

In this figure, AWA ranges from 20.3 degrees to 37.9 degrees for the port tack. The former is the closest angle to the wind that was achieved, and the latter is typical of a close reaching condition, where the sail is eased for the power down mode. There is some scatter in the experimental data because this is made up from measurements taken with the sails trimmed in slightly different ways. The experimental value of *CL* in Fig.13(a) varies with AWA from 0.91 to 1.58. For the close reaching condition, the sails were not well trimmed to satisfy the power down mode. A sample of measured sail sections at this condition is shown in a figure attached to Table 2(2). From the figure, it can be seen that both the mainsail and the jib are not eased sufficiently to correspond to the large AWA. This is the reason for the decrement in the measured lift curve slope of *CL* at the range of AWA angles over about 35 degrees.


Table 2. Sail shapes, measured experimental data and three-dimensional coordinates of the sails for the cases of (1) 96092335 and (2) 96080248.

In this figure, AWA ranges from 20.3 degrees to 37.9 degrees for the port tack. The former is the closest angle to the wind that was achieved, and the latter is typical of a close reaching condition, where the sail is eased for the power down mode. There is some scatter in the experimental data because this is made up from measurements taken with the sails trimmed in slightly different ways. The experimental value of *CL* in Fig.13(a) varies with AWA from 0.91 to 1.58. For the close reaching condition, the sails were not well trimmed to satisfy the power down mode. A sample of measured sail sections at this condition is shown in a figure attached to Table 2(2). From the figure, it can be seen that both the mainsail and the jib are not eased sufficiently to correspond to the large AWA. This is the reason for the decrement in the measured lift curve slope of *CL* at the range of AWA angles over about 35 degrees.

Table 2. Sail shapes, measured experimental data and three-dimensional coordinates of the

sails for the cases of (1) 96092335 and (2) 96080248.

The calculated results for *CL* using the VLM show good agreement with the experiments at AWA angles less than about 35 degrees. Over about 35 degrees, the calculated results are lower than the measured ones. This shows that the calculated results strongly indicate the effect of incorrect sail trimming. The results for *CL* using the RANS-based CFD show the same trends with the experiments, but are slight higher than those from the experiments for AWA between 20 degrees to 30 degrees and lower for AWA greater than 30 degrees. In particular, the decrease in *CL* for AWA values greater than 30 degrees is considerably large. This will be discussed later with the calculated sail surface pressure and streamlines. The calculated results for *CD* slightly over predict those from the experiments. Fig.13(c) shows the coordinates of the center of effort of the sails. The *x* and *z* coordinates of the geometric center of effort (*xGCE* and *zGCE*) are 0.63m aft and 4.80m above the origin, which are indicated by alternate long and short dashed lines in the figure. It is seen that both the experimental and the calculated coordinates of *xCE* are near *xGCE* and move slightly forward with increasing AWA. Unfortunately, there is a wide scatter in the experimental values of *zCE*. This is thought to be because the measured *Ks* moment contains a large component from the mass of the dynamometer frame and rigging (659kg). This moment was subtracted from the measurement, taking into account the measured heel angle. If there is a slight error in the position of center of gravity of the dynamometer frame, or in the measured heel angle, the error in the calculated moment will be large. However, though there is a scatter in the measured data, it can be seen that *zCE* is decreasing as AWA increases. The trends in the movement of both *xCE* and *zCE* as functions of AWA might be caused by the decrement of force acting on the aft and upper parts of the sails due to the loosening of main and jib sheets with increasing AWA. The calculated results for *zCE* obtained using theRANS-based CFD show the same trend as the experiments. On the other hand, the calculated results using VLM are considerably higher than the experimental ones. This might be caused by over estimation of the force acting on the upper portion of the mainsail. In this area, since the jib is not overlapping, flow separation may occur easily. However, the VLM does not take flow separation into account.

Figures 14(1) and 14(2) show the calculated results of the sail surface pressure and streamlines using RANS- based CFD. Fig.14(1) indicates the case of experiment ID 96092335 (AWA= 30.7deg.), and 14(2) indicates ID 96080248 (AWA= 37.9deg.). These data correspond to the plotted points on the vertical dotted lines (1) and (2) in Fig.13. In Fig.14, the left and right diagrams correspond to the port and starboard sides, i.e., pressure and suction sides, respectively. In 14(1), although slight flow separation on the suction side of mainsail is seen, the streamlines of both sides run smoothly. On the other hand, in 14(2), considerable flow separation is occurring, in particular, on the suction side of jib. This is the main reason for the reduction of *CL* value in the RANS-based CFD calculation at (2) in Fig.13(a). The shapes and three-dimensional coordinates of the sails are given in Table 2. The numbered (1) and (2) tables correspond to the cases of experiment ID 96092335 and ID 96080248, respectively. These also correspond to the calculated results shown in Fig.14. The figures described above the tables show the sail section profiles at 0, 20, 40, 60 and 80% of the sail height. The dimensions of these three-dimensional coordinates are given in the tables including 100% height section data. The positive direction of the *x* coordinate is aft. The four lines at the top of the tables are the measured values for the wind and sail trim conditions, the boat attitude and the sail performance coefficients.

Sail Performance Analysis of Sailing Yachts by Numerical Calculations and Experiments 113

twist angle ranges from 4.5 degrees to 24.9 degrees for the port tack. Varying the twist angle by 20.4 degrees, results in the value of *CX* in Fig. 15(b) changing from 0.33 to 0.39 (18%), and the value of *CY* changing from 1.16 to 1.39 (20%). It can be seen that the maximum *CX* occurs at a twist angle of around 15 degrees. The considerable decrease in *CY* with increasing twist angle is also worth noticing. In this case, the calculated results of both VLM and RANSbased CFD for *CX* and *CY*, and *CL* and *CD* correspond to the measured values very well.

Fig. 15. Performance variation as a function of mainsail twist angle for mainsail and 130% jib. Figures 16(1) and 16(2) show the calculated results using RANS-based CFD. Fig.16(1) corresponds to ID 97072213 (twist angle = 8.2 deg.), and 16(2) corresponds to ID 97072218 (twist angle = 24.1 deg.). It can be seen in Fig.16(1) that the streamlines on the upper part of the suction side of the mainsail for the smaller twist angle, show considerable flow

Fig. 14. (1) Surface pressure and streamlines obtained by RANS-based CFD at experimental ID 96092335 (AWA=30.7 deg.) and (2) ID 96080248 (AWA=37.9 deg.).

#### **5.3.2 Variation with mainsail twist angle**

Fig.15 shows the performance variation for the mainsail and 130% jib configuration as a function of mainsail twist angle. The mainsail twist was changed by varying the main sheet tension. The boom angle was kept parallel with the boat centerline by moving the main sheet traveler. The experiment was performed for an average value of AWA of 30 ± 2 degrees and mean draft at around 10%. The jib shape was fixed. The twist angle is defined as the angle between the boom line and section chord line at 80% height. In the figure, the

Fig. 14. (1) Surface pressure and streamlines obtained by RANS-based CFD at experimental

Fig.15 shows the performance variation for the mainsail and 130% jib configuration as a function of mainsail twist angle. The mainsail twist was changed by varying the main sheet tension. The boom angle was kept parallel with the boat centerline by moving the main sheet traveler. The experiment was performed for an average value of AWA of 30 ± 2 degrees and mean draft at around 10%. The jib shape was fixed. The twist angle is defined as the angle between the boom line and section chord line at 80% height. In the figure, the

ID 96092335 (AWA=30.7 deg.) and (2) ID 96080248 (AWA=37.9 deg.).

**5.3.2 Variation with mainsail twist angle** 

(1)

(2)

twist angle ranges from 4.5 degrees to 24.9 degrees for the port tack. Varying the twist angle by 20.4 degrees, results in the value of *CX* in Fig. 15(b) changing from 0.33 to 0.39 (18%), and the value of *CY* changing from 1.16 to 1.39 (20%). It can be seen that the maximum *CX* occurs at a twist angle of around 15 degrees. The considerable decrease in *CY* with increasing twist angle is also worth noticing. In this case, the calculated results of both VLM and RANSbased CFD for *CX* and *CY*, and *CL* and *CD* correspond to the measured values very well.

Fig. 15. Performance variation as a function of mainsail twist angle for mainsail and 130% jib.

Figures 16(1) and 16(2) show the calculated results using RANS-based CFD. Fig.16(1) corresponds to ID 97072213 (twist angle = 8.2 deg.), and 16(2) corresponds to ID 97072218 (twist angle = 24.1 deg.). It can be seen in Fig.16(1) that the streamlines on the upper part of the suction side of the mainsail for the smaller twist angle, show considerable flow

Sail Performance Analysis of Sailing Yachts by Numerical Calculations and Experiments 115

Fig. 16. (1) Surface pressure and streamlines obtained by RANS-based CFD at experimental

The flow is dominated by multiple-lifting-surface aerodynamic interactions. For larger AWA values, in particular, a large-scale flow separation exists on the leeward side of the sails. In general, there is complex vortex generation in the wake, especially near the top and bottom of the sails, i.e., tip vortices are generated and are influenced by the boundary layer flows on the sails. The resultant aerodynamic forces are mostly dominated by the pressure component, whereas the contribution of the frictional component is generally small. The accurate prediction of the boundary layer flows on the sails and the three-dimensional flow separation, associated with the abovementioned vortex generation, are big challenges for

ID 97072213 (twist angle= 8.2 deg.) and (2) ID 97072218 (twist angle= 24.1 deg.).

**6. Discussion of numerical calculation methods** 

(1)

(2)

separation. This is caused by the large angle of attack at the upper part of the sail due to the small twist angle. On the other hand, for the higher twist angle shown in Fig.16(2), there is a low negative pressure area at the luff on the suction side of mainsail due to the small angle of attack. This is what causes the considerable reduction in the calculated value for *CX* at (2) in Fig.15(b). Table 3 shows the shapes and three-dimensional coordinates of the sails for cases (1) and (2), which correspond to the calculated results shown in Fig.16. Further measured data and comparison with the numerical calculations are described by Masuyama et al. (2007 and 2009)


Table 3. Sail shapes, measured experimental data and three-dimensional coordinates of the sails for the cases of (1) 97072213 and (2) 97072218.

separation. This is caused by the large angle of attack at the upper part of the sail due to the small twist angle. On the other hand, for the higher twist angle shown in Fig.16(2), there is a low negative pressure area at the luff on the suction side of mainsail due to the small angle of attack. This is what causes the considerable reduction in the calculated value for *CX* at (2) in Fig.15(b). Table 3 shows the shapes and three-dimensional coordinates of the sails for cases (1) and (2), which correspond to the calculated results shown in Fig.16. Further measured data and comparison with the numerical calculations are described by Masuyama

Table 3. Sail shapes, measured experimental data and three-dimensional coordinates of the

sails for the cases of (1) 97072213 and (2) 97072218.

et al. (2007 and 2009)

Fig. 16. (1) Surface pressure and streamlines obtained by RANS-based CFD at experimental ID 97072213 (twist angle= 8.2 deg.) and (2) ID 97072218 (twist angle= 24.1 deg.).

#### **6. Discussion of numerical calculation methods**

The flow is dominated by multiple-lifting-surface aerodynamic interactions. For larger AWA values, in particular, a large-scale flow separation exists on the leeward side of the sails. In general, there is complex vortex generation in the wake, especially near the top and bottom of the sails, i.e., tip vortices are generated and are influenced by the boundary layer flows on the sails. The resultant aerodynamic forces are mostly dominated by the pressure component, whereas the contribution of the frictional component is generally small. The accurate prediction of the boundary layer flows on the sails and the three-dimensional flow separation, associated with the abovementioned vortex generation, are big challenges for

Sail Performance Analysis of Sailing Yachts by Numerical Calculations and Experiments 117

moment acting on the model are measured simultaneously. The numerical simulation by using a RANS-based CFD method is also in progress. Along with integration with the aforementioned sail design and performance prediction software *AAF*, more details of our

The sail dynamometer boat *Fujin* was built for sail tests for the Japanese America's Cup entry by a Grant-in-Aid from the Nippon Foundation and the authors would like to thank the Nippon Foundation for providing them this invaluable tool. The authors wish to express their thanks to Yamaha Motor Co. Ltd. for building *Fujin* and to North Sails Japan Co. for making the sails. We would like to thank Dr. Martin Renilson for his valuable discussions and comments on this article. We would also like to thank Mr. H. Mitsui, the harbormaster of the Anamizu Bay Seminar House of the Kanazawa Institute of Technology, for his assistance with the sea trials. Help with the sea trials given by graduate and undergraduate students of the Kanazawa Institute of Technology is also gratefully acknowledged. The graduate students were Masaya Miyagawa, Takashi Hasegawa, and Munehiko Ogihara. Finally, we would like to acknowledge the effort by Mr. Naotoshi Maeda who carried out most of the RANS-CFD simulations as a part of his work on M.S. thesis at Osaka Prefecture

Baldwin, B.S. & Lomax, H. (1978). Thin Layer Approximation and Algebraic Model for

Falkner, V.M. (1943). The Calculation of Aerodynamic Loading on Surface of Any Shape,

Falkner, V.M. (1946). The Accuracy of Calculations Based on Vortex Lattice Theory, *British* 

Fukasawa, T. & Katori, M. (1993). Numerical Approach to Aeroelastic Responses of Three-

Fukasawa, T. (1993). Aeroelastic Transient Response of 3-Dimensional Flexible Sail,

Hansen, H., Jackson, P. & Hochkirch, K. (2003). Comparison of Wind Tunnel and Full-scale

Hochkirch, K. & Brandt, H. (1999). Fullscale Hydrodynamic Force Measurement on the

Katori, M. (2009). Advanced Aero Flow Software Manual, *North Sails Japan*, December 2009. Krebber, B. & Hochkirch, K. (2006). Numerical Investigation on the Effects of Trim for a

*Symposium*, pp.33-44, Annapolis, Maryland, USA, January, 1999

Separated Turbulent Flows, *AIAA Paper, 78-257, AIAA 16th Aerospace Science* 

*Aeronautical Research Committee Reports and Memoranda No.1910,* Ministry of Aircraft

Dimensional Flexible Sails, *Proceedings of 11th Chesapeake Sailing Yacht Symposium*,

*Proceedings of International Conference on Aero-Hydroelasticity*, pp.57-62, Beijing,

Aerodynamic Sail Force, *International Journal of Small Craft Technology* (IJSCT), Vol.

Berlin Sailing Dynamometer, *Proceedings of SNAME 14th Chesapeake Sailing Yacht* 

Yacht Rig, *Proceedings of 2nd High Performance Yacht Design Conference*, pp.13-21,

work will be reported in our future publications.

*Meeting*, Reno, NV., USA, 1978

*Aeronautical Research Council,* No.9621

pp.87-105, Annapolis, USA, January 29-30, 1993

Production, London, UK

China, October 18-21, 1993,

145 Part B1, (2003), pp. 23-31

Auckland, New Zealand, February, 2006

**8. Acknowledgments** 

University.

**9. References** 

RANS-based CFD. The geometrical complexity is also another significant challenge to RANS-based CFD. The accuracy in the prediction of the CE is of great interest, in association with the correct prediction of the above-mentioned three-dimensional flow separation.

Through the analyses of the multiblock RANS-based CFD, it appears that the overall trends of the flow and the aerodynamic forces measured in the experiments are fairly well predicted by the present computations. It is also seen that the multiblock domain decomposition considered here is very effective for the present mainsail and jib configurations. The automatic gridding scheme used successfully generates high-quality structured grids for the various sail geometries, AWA, and heel angles considered in the present study. Although there are advantages to a structured grid system for highresolution in the boundary layer flow, building a grid in this fashion is difficult to apply to complex geometries. This problem appears to be resolved by the present scheme.

The Vortex Lattice method is, on the other hand, a convenient tool to predict the lift and induced drag acting on the sail accurately at the apparent wind angles less than about 35 degrees. The computational time of the method is about a few minutes for one calculation condition. The longitudinal coordinates, or x coordinate, of the center of effort of the sail can also be calculated with accuracy by the Vortex Lattice method, however, the estimated vertical coordinate, or z coordinate, of the center of effort by the Vortex Lattice method is considerably higher than the experimental ones. This may be caused by the fact that the flow at the upper portion of mailsail is easily separated because of the absence of jib overlapping, while the flow separation cannot be taken into account in the Vortex Lattice method

#### **7. Conclusion**

The sail performance analysis of sailing yacht was carried out by using numerical calculations and experiments. Focus in the present manuscript is especially on the upwind sailing condition. The sails considered here are IMS type, and the shapes and performance were measured by using the sail dynamometer boat *Fujin*. The measured sail flying shapes were used by numerical analysis, where two CFD methods developed by the authors were used, i.e. a multiblock RANS-based CFD method by Tahara and a VLM-based CFD method by Fukasawa. It appears that the overall trends of the flow and the aerodynamic forces measured in the experiments are fairly well predicted by the present computations; and at the same time, the present sail performance database based on the full scale onboard measurements are very useful for validation study of numerical methods. As compared to maturity of VLM, that for RANS-based CFD is still in underway but the future prospect is shown promising, especially for capability in predicting separation flow filed where viscous effects of fluid are significant. The authors believe that our sail performance database associated with accurate sail flying shape measurements will be able to contribute to the further development of more advanced CFD methods.

Although details are not described in the present manuscript, our current effort is directed toward the more challenging problem, i.e., extension of the present work for the downwind sailing condition. Since the onboard sail shape measurement system of *Fujin* is incompetent for the spinnaker measurement due to its balloon shape, the sail shapes and performance are measured using wind tunnel equipment, and such activities are already in progress. The sail shapes are recorded using digital cameras and processed to obtain 3D coordinates using solid shape analyzer software, which provides 3D coordinates from digital photographs taken from several different directions. Importantly, the sail shapes and the sail forces and moment acting on the model are measured simultaneously. The numerical simulation by using a RANS-based CFD method is also in progress. Along with integration with the aforementioned sail design and performance prediction software *AAF*, more details of our work will be reported in our future publications.

#### **8. Acknowledgments**

116 Fluid Dynamics, Computational Modeling and Applications

RANS-based CFD. The geometrical complexity is also another significant challenge to RANS-based CFD. The accuracy in the prediction of the CE is of great interest, in association with the correct prediction of the above-mentioned three-dimensional flow separation. Through the analyses of the multiblock RANS-based CFD, it appears that the overall trends of the flow and the aerodynamic forces measured in the experiments are fairly well predicted by the present computations. It is also seen that the multiblock domain decomposition considered here is very effective for the present mainsail and jib configurations. The automatic gridding scheme used successfully generates high-quality structured grids for the various sail geometries, AWA, and heel angles considered in the present study. Although there are advantages to a structured grid system for highresolution in the boundary layer flow, building a grid in this fashion is difficult to apply to

complex geometries. This problem appears to be resolved by the present scheme.

separation cannot be taken into account in the Vortex Lattice method

further development of more advanced CFD methods.

**7. Conclusion** 

The Vortex Lattice method is, on the other hand, a convenient tool to predict the lift and induced drag acting on the sail accurately at the apparent wind angles less than about 35 degrees. The computational time of the method is about a few minutes for one calculation condition. The longitudinal coordinates, or x coordinate, of the center of effort of the sail can also be calculated with accuracy by the Vortex Lattice method, however, the estimated vertical coordinate, or z coordinate, of the center of effort by the Vortex Lattice method is considerably higher than the experimental ones. This may be caused by the fact that the flow at the upper portion of mailsail is easily separated because of the absence of jib overlapping, while the flow

The sail performance analysis of sailing yacht was carried out by using numerical calculations and experiments. Focus in the present manuscript is especially on the upwind sailing condition. The sails considered here are IMS type, and the shapes and performance were measured by using the sail dynamometer boat *Fujin*. The measured sail flying shapes were used by numerical analysis, where two CFD methods developed by the authors were used, i.e. a multiblock RANS-based CFD method by Tahara and a VLM-based CFD method by Fukasawa. It appears that the overall trends of the flow and the aerodynamic forces measured in the experiments are fairly well predicted by the present computations; and at the same time, the present sail performance database based on the full scale onboard measurements are very useful for validation study of numerical methods. As compared to maturity of VLM, that for RANS-based CFD is still in underway but the future prospect is shown promising, especially for capability in predicting separation flow filed where viscous effects of fluid are significant. The authors believe that our sail performance database associated with accurate sail flying shape measurements will be able to contribute to the

Although details are not described in the present manuscript, our current effort is directed toward the more challenging problem, i.e., extension of the present work for the downwind sailing condition. Since the onboard sail shape measurement system of *Fujin* is incompetent for the spinnaker measurement due to its balloon shape, the sail shapes and performance are measured using wind tunnel equipment, and such activities are already in progress. The sail shapes are recorded using digital cameras and processed to obtain 3D coordinates using solid shape analyzer software, which provides 3D coordinates from digital photographs taken from several different directions. Importantly, the sail shapes and the sail forces and The sail dynamometer boat *Fujin* was built for sail tests for the Japanese America's Cup entry by a Grant-in-Aid from the Nippon Foundation and the authors would like to thank the Nippon Foundation for providing them this invaluable tool. The authors wish to express their thanks to Yamaha Motor Co. Ltd. for building *Fujin* and to North Sails Japan Co. for making the sails. We would like to thank Dr. Martin Renilson for his valuable discussions and comments on this article. We would also like to thank Mr. H. Mitsui, the harbormaster of the Anamizu Bay Seminar House of the Kanazawa Institute of Technology, for his assistance with the sea trials. Help with the sea trials given by graduate and undergraduate students of the Kanazawa Institute of Technology is also gratefully acknowledged. The graduate students were Masaya Miyagawa, Takashi Hasegawa, and Munehiko Ogihara. Finally, we would like to acknowledge the effort by Mr. Naotoshi Maeda who carried out most of the RANS-CFD simulations as a part of his work on M.S. thesis at Osaka Prefecture University.

#### **9. References**


**Part 2** 

**Multiphase Flow, Structures, and Gases** 

