# On the Analytical Approach to Present Engineering Problems: Photovoltaic Systems Behavior, Wind Speed Sensors Performance, and High-Speed Train Pressure Wave Effects in Tunnels

Authors: Santiago Pindado; Javier Cubas; Félix Sorribes-Palmer
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

At present, engineering problems required quite a sophisticated calculation means. However, analytical models still can prove to be a useful tool for engineers and scientists when dealing with complex physical phenomena. The mathematical models developed to analyze three different engineering problems: photovoltaic devices analysis; cup anemometer performance; and high-speed train pressure wave effects in tunnels are described. In all cases, the results are quite accurate when compared to testing measurements.


---

## Body

## 1. Introduction

The evolution of engineering has provided engineers with better tools, that is, better calculation possibilities (more sophisticated numerical approximations and computers) and better testing equipment. However, the classical approach, based on mathematical models directly linked to the physics of the analyzed phenomena, should not be left aside. Even in the rush to obtain a suitable solution to an engineering problem, the analytical approximations can state the limits to accept the solution from more complex calculations. Besides, a very good explanation for the need of analytical methods is offered by Simmonds and Mann [1] in relation to the use of perturbation methods, as they “produce analytic approximations that often reveal the essential dependence of the exact solution on the parameters in a more satisfactory way than does a numerical solution.” It cannot be more properly expressed.At the IDR/UPM Institute different analytical models have been developed along the last decades to study quite complicated problems. As an example, the stability analysis of liquid bridges has been thoroughly studied since the late 70s, when the effect of microgravity on the liquid bridges was brought to Professor Da Riva’s and Professor Meseguer’s attention [2–9]. Another good example is the analysis of train-induced pressure effects on pedestrian and sign panels carried out by Professor Sanz-Andrés [10–12].At present, several problems are analyzed at the IDR/UPM Institute under the guidance of Professor Meseguer. In the following sections the analytical models developed to study three different phenomena are described. In Section2, the analytical models used to study solar cells/panels are described. These models have been developed to analyze the solar panels of the UPMSat-2 satellite and model its power subsystem. This 50 kg microsatellite is planned to be launched by the end of 2015/beginning of 2016. In the present work, a new and a bit more complex analytical approximation to solar panels behavior has been developed. This approximation is based on a 2-diode equivalent electrical circuit, which is obviously less simple than the one based on the 1-diode equivalent electrical circuit already used at the IDR/UPM to estimate photovoltaic power supply [13–16].Additionally, a model used to study the rotor aerodynamics of the cup anemometer is described in Section3. The cup anemometer, invented by J. T. R. Robinson in the 19th century, is today the standardized instrument regarding the wind power performance measurements [17]. It has been extensively studied and optimized during the 20th century (see in [18] a thorough review of the literature) and, at present, it is the reference instrument in relation to wind production forecasting in the field and wind turbine performance control. Since 2009 the performance of this instrument is studied at the IDR/UPM Institute [18–26]. As a result, a quite accurate analytical method has been developed to analyze the anemometer performance based on the cups’ aerodynamics, including the effect of perturbations. In the present work, the analytical model is used to study the equilibrium positions (hereinafter referred as “flag” positions) of a damaged anemometer with one cup missing. This is not a rare case, bearing in mind that 30% of mast-mounted anemometers return for recalibration in bad shape [27].Finally, a mathematical model developed to study the effect of pressure waves in high-speed railway tunnels is included in Section4. This model can be used in preliminary design of tunnel sections and to evaluate tunnel costs and also to find the optimal location for damper devices depending on the trains and the travel speeds during the pass through the tunnel. Conclusions are summarized in Section 5.


## 2. Solar Cell/Panel Modeling

The behavior of a solar cell is normally modeled by a 1-diode and 2-resistor electric circuit (see Figure1); this simple approximation captures the nature of the more important physical effects related to the photovoltaic conversion of the light to current. However, this model that fits quite correctly the behavior of monocrystalline and polycrystalline silicon cells has certain deviation when the solar cell is not operating at standard test conditions [28, 29]. It should also be fair to say that today this model can be considered a simplification from the more accurate 2-diode and 2-resistor electric circuit model (see Figure 1), which was proposed by Wolf and Rauschenbach in 1961 to introduce the effect of recombination of electrons and holes in the depletion region [30]. These two models have been successfully applied to model solar panels made of cells grouped in parallel series [31, 32]. In the following subsection, the 1-diode/2-resistor equivalent circuit model is firstly described as an introduction to the second subsection, in which the 2-diode/2-resistor equivalent circuit model developed in the present work is presented.Figure 1
Sketches of the 1-diode/2-resistor electric circuit model (a) and 2-diode/2-resistor electric circuit model (b).
(a)
(b)The use of these models is conditioned by the identification of the parameters of the equivalent circuits. Two main approaches are used to fulfill this task, analytical and numerical. Although many more numerical solutions can be found in the literature [33–44] than the ones based on analytical reasoning [13, 32, 45–51], this last approach is still used, especially when high degree of reliability is needed within the model recalculations (e.g., calculating the power supply in operational spacecraft photovoltaic subsystems [52, 53]).


### 2.1.  The 1-Diode/2-Resistor Electric Circuit Model

The equation that describes the behavior of photovoltaic devices following this approximation is(1)
I
=
I
pv
-
I
0
exp
⁡
V
+
I
R
s
a
V
T
-
1
-
V
+
I
R
s
R
sh
,where I
pv is the photovoltaic current, R
s is the series resistance, R
sh is the shunt resistance, I
0 is the saturation current of the diode, a is the ideality factor of the diode, and V
T is the thermal voltage of the cell, defined as(2)
V
T
=
k
T
q
,where T is the temperature, k is the Boltzmann constant, and q is the charge of the electron.According to (1), five parameters (I
pv, I
0, R
s, R
sh, and a) must be identified before using the model to calculate the performance of the studied photovoltaic device. The analytical approach is normally based on data from the most representative points (short circuit: V
=
0, I
=
I
sc; open circuit: V
=
V
oc, I
=
0; and maximum power: V
=
V
mp, I
=
I
mp points) of the measured I-V curve of the solar cell/panel (as an example, see in Figure 2 the current-voltage and power-voltage curves of a triple-junction solar cell (Emcore ZTJ), developed in recent years for space applications [54–58]. The characteristic points of the I-V curve are I
sc
=
0.462 A; V
oc
=
2.73 V; and V
mp
=
2.41 V and I
mp
=
0.441 A [56]). The aforementioned three points help to derive four equations which allow identifying the value of four of the aforementioned parameters in relation to the fifth (normally, this fifth parameter is the ideality factor, a).Figure 2
Emcore ZTJ triple-junction (InGaP/InGaAs/Ge) solar cell current-voltage (I-V) and power-voltage curves (solid and dashed lines, resp.), measured at AM0 (135.3 mW/cm²) and 28°C.With data (current and voltage levels) from the short circuit and the open circuit points of theI-V curve and taking into account the magnitude of each one of the terms, it is possible to derive the following equations [13]:(3)
I
pv
=
R
sh
+
R
s
R
sh
I
sc
,
(4)
I
0
=
R
sh
+
R
s
I
sc
-
V
oc
R
sh
exp
⁡
V
oc
/
a
V
T
,whereas from the maximum power an additional equation is obtained(5)
I
mp
=
I
sc
-
I
sc
-
V
oc
-
R
s
I
sc
R
sh
exp
⁡
V
mp
+
I
mp
R
s
-
V
oc
a
V
T
-
V
mp
+
I
mp
R
s
-
I
sc
R
s
R
sh
.Finally, taking into account the condition of the maximum power at the maximum power point(6)
∂
P
∂
V
=
V
∂
I
∂
V
+
I
=
0an implicit expression of the series resistor, R
s, as a function of the initial parameters can be obtained:(7)
a
V
T
V
mp
2
I
mp
-
I
sc
·
V
mp
I
sc
+
V
oc
I
mp
-
I
sc
h
h
h
h
h
h
·
V
mp
-
I
mp
R
s
-
a
V
T
V
mp
I
sc
-
V
oc
I
mp
-
1
=
exp
⁡
V
mp
+
I
mp
R
s
-
V
oc
a
V
Ttogether with another for the shunt resistor, R
sh, as a function of R
s and the initial parameters:(8)
R
sh
=
V
mp
-
I
mp
R
s
V
mp
-
R
s
I
sc
-
I
mp
-
a
V
T
V
mp
-
I
mp
R
s
I
sc
-
I
mp
-
a
V
T
I
mp
.As aforementioned, parametera needs to be estimated in order to solve the above equations, and one of them is unfortunately implicit for R
s. One solution to this difficulty (leaving aside any numerical help) consists in using the Lambert W-function, W
(
z
), which is defined as(9)
z
=
W
z
e
W
(
z
)
,where z is any complex number [59]. This function is not injective; in the real variable x the relation W is defined only for x
≥
-
1
/
e and is double valued in the bracket [
-
1
/
e
,
0
], the two branches of the function being expressed as W
0
(
x
) for W
(
x
)
≥
-
1 and W
-
1
(
x
) for W
(
x
)
≤
-
1 in the aforementioned bracket. The general strategy to apply the Lambert W-function in solving exponential equations is to use the following equivalence:(10)
X
=
Y
e
Y
⟺
Y
=
W
X
.Equation (7) can be rewritten as follows:(11)
-
V
mp
2
I
mp
-
I
sc
V
mp
I
sc
+
V
oc
I
mp
-
I
sc
·
exp
⁡
⁡
-
2
V
mp
-
V
oc
a
V
T
+
V
mp
I
sc
-
V
oc
I
mp
V
mp
I
sc
+
V
oc
I
mp
-
I
sc
=
I
mp
R
s
-
V
mp
a
V
T
+
V
mp
I
sc
-
V
oc
I
mp
V
mp
I
sc
+
V
oc
I
mp
-
I
sc
·
exp
⁡
⁡
I
mp
R
s
-
V
mp
a
V
T
+
V
mp
I
sc
-
V
oc
I
mp
V
mp
I
sc
+
V
oc
I
mp
-
I
scand then, using the mathematical equivalence from (10), the above equation turns into(12)
I
mp
R
s
-
V
mp
a
V
T
+
V
mp
I
sc
-
V
oc
I
mp
V
mp
I
sc
+
V
oc
I
mp
-
I
sc
=
W
-
1
-
V
mp
2
I
mp
-
I
sc
V
mp
I
sc
+
V
oc
I
mp
-
I
sc
h
h
h
h
h
h
h
h
·
exp
⁡
V
mp
I
sc
-
V
oc
I
mp
V
mp
I
sc
+
V
oc
I
mp
-
I
sc
-
2
V
mp
-
V
oc
a
V
T
h
h
h
h
h
h
h
h
h
h
h
h
.
h
+
V
mp
I
sc
-
V
oc
I
mp
V
mp
I
sc
+
V
oc
I
mp
-
I
sc
,where W
-
1 is the negative branch of the Lambert W-function (as the left part of (14) is lower than −1 for typical cells and solar panels). And then an explicit expression for R
s is obtained:(13)
R
s
=
A
W
-
1
B
exp
⁡
⁡
C
-
D
+
C
,where (14)
A
=
a
V
T
I
mp
,
B
=
-
V
mp
2
I
mp
-
I
sc
V
mp
I
sc
+
V
oc
I
mp
-
I
sc
,
C
=
-
2
V
mp
-
V
oc
a
V
T
+
V
mp
I
sc
-
V
oc
I
mp
V
mp
I
sc
+
V
oc
I
mp
-
I
sc
,
D
=
V
mp
-
V
oc
a
V
T
.An alternative to the estimation of the ideality factor,a, consists in using other specific characteristics of the I-V curve, such as the slope at short circuit and open circuit points, R
sh
0 and R
s
0, respectively,(15)
R
sh
0
=
-
∂
V
∂
I
I
=
I
sc
,
(16)
R
s
0
=
-
∂
V
∂
I
V
=
V
octo calculate the ideality factor, a. It is preferable to use (15) rather than (16), as the final derived expressions are less complicated. Once R
sh
0 has been obtained, (1) at the short circuit point of the I-V curve can be expressed as [49](17)
1
R
sh
-
1
R
sh
0
-
R
s
+
I
0
a
V
T
exp
⁡
R
s
I
sc
a
V
T
=
0which, taking into account a magnitude analysis of the three terms (in [13] the order of magnitude of terms from the above equation was compared for different technologies of solar cells; R
sh
-
1 and R
sh
0
-
R
s
-
1 had very similar values, whereas the term I
0
/
a
V
T
exp
⁡
R
s
I
sc
/
a
V
T was, in the best case, around 60 times lower), can be simplified as(18)
1
R
sh
-
1
R
sh
0
-
R
s
=
0
;that is,(19)
R
sh
0
=
R
sh
+
R
s
.Once the above expression has been reached, it is possible to derive an expression that gives the ideality factor as a function of constants already known,R
s and the slope R
sh
0:(20)
a
V
T
=
V
mp
-
I
mp
R
s
V
mp
+
I
mp
-
I
sc
R
sh
0
V
mp
-
I
mp
R
sh
0
.Finally, from (7) and (20) the value of the series resistor, R
s, is obtained:(21)
R
s
=
E
-
F
E
+
F
V
mp
I
mp
+
F
E
+
F
V
oc
I
mp
,where(22)
E
=
V
mp
+
I
mp
-
I
sc
R
sh
0
log
⁡
V
mp
+
I
mp
-
I
sc
R
sh
0
V
oc
-
I
sc
R
sh
0
,
F
=
V
mp
-
R
sh
0
I
mp
.As said in Section1, this 1-diode/2-resistor electric circuit model has been successfully correlated to solar cells behavior [13–16].


### 2.2. The 2-Diode/2-Resistor Electric Circuit Model

As this equivalent circuit has one more diode, the model equation contains two exponential terms:(23)
I
=
I
pv
-
I
01
exp
⁡
V
+
I
R
s
a
1
V
T
-
1
-
I
02
exp
⁡
V
+
I
R
s
a
2
V
T
-
1
-
V
+
I
R
s
R
sh
;therefore, four new parameters, two corresponding to different saturation currents, I
01 and I
02, together with other two corresponding to the ideality factors of each diode, a
1 and a
2, are introduced in the calculations instead of the single saturation current and ideality factor from the 1-diode/2-resistor model equation. To identify the parameters of the above equation, the three characteristic points of the I-V curve are considered as boundary conditions, together with the slope of the curve at two of them. As a consequence, six equations are derived to identify seven parameters and, therefore, one of these parameters needs to be estimated.Operating as in the previous subsection, the short circuit point of theI-V curve leads to the following expression: (24)
I
sc
=
I
pv
-
I
01
exp
⁡
I
sc
R
s
a
1
V
T
-
1
-
I
02
exp
⁡
I
sc
R
s
a
2
V
T
-
1
-
I
sc
R
s
R
shwhich, taking into account the much lower magnitude of the terms involving the exponentials, can be simplified as the same relationship (3) obtained with the 1-diode/2-resistor model:(25)
I
pv
=
R
sh
+
R
s
R
sh
I
sc
.Besides, the derivative of (23) at the short circuit point can be rewritten as(26)
1
+
R
sh
+
R
s
(
d
I
/
d
V
)
sc
R
sh
1
+
R
s
(
d
I
/
d
V
)
sc
=
-
1
a
1
V
T
I
01
exp
⁡
I
sc
R
s
a
1
V
T
-
1
a
2
V
T
I
02
exp
⁡
I
sc
R
s
a
2
V
Twhich, bearing in mind the slope of the I-V curve (i.e., (15)), leads to (27)
R
sh
0
-
R
sh
+
R
s
R
sh
R
sh
0
-
R
s
=
-
1
a
1
V
T
I
01
exp
⁡
I
sc
R
s
a
1
V
T
-
1
a
2
V
T
I
02
exp
⁡
I
sc
R
s
a
2
V
T
.Again, comparing the magnitudes of the different terms in the above equation, the already known equation (19) can be obtained:(28)
R
sh
0
=
R
sh
+
R
s
.From the open circuit point, the following expression can be obtained:(29)
R
sh
+
R
s
I
sc
-
V
oc
R
sh
=
I
01
exp
⁡
V
oc
a
1
V
T
+
I
02
exp
⁡
V
oc
a
2
V
T
.If the derivative of (23) at that point(30)
1
+
R
sh
+
R
s
d
I
/
d
V
oc
R
sh
1
+
R
s
d
I
/
d
V
oc
=
-
1
a
1
V
T
I
01
exp
⁡
V
oc
a
1
V
T
-
1
a
2
V
T
I
02
exp
⁡
V
oc
a
2
V
Tis considered, together with the slope of the I-V (i.e., (16)), it is possible to derive the following equation:(31)
R
sh
+
R
s
-
R
s
0
R
sh
R
s
0
-
R
s
=
1
a
1
V
T
I
01
exp
⁡
V
oc
a
1
V
T
+
1
a
2
V
T
I
02
exp
⁡
V
oc
a
2
V
T
.From the maximum power point data the equation(32)
V
mp
+
I
mp
-
I
sc
R
sh
+
R
s
R
sh
=
-
I
01
exp
⁡
V
mp
+
I
mp
R
s
a
1
V
T
-
I
02
exp
⁡
V
mp
+
I
mp
R
s
a
2
V
Tis obtained. Furthermore, bearing in mind the maximum power condition (6)(33)
∂
I
∂
V
I
mp
,
V
mp
=
-
I
mp
V
mpit is possible to derive the last expression needed:(34)
1
+
R
sh
+
R
s
d
I
/
d
V
mp
R
sh
1
+
R
s
d
I
/
d
V
mp
=
-
1
a
1
V
T
I
01
exp
⁡
V
mp
+
I
mp
R
s
a
1
V
T
-
1
a
2
V
T
I
02
exp
⁡
V
mp
+
I
mp
R
s
a
2
V
T
,
(35)
R
sh
+
R
s
I
mp
-
V
mp
R
sh
V
mp
-
R
s
I
mp
=
1
a
1
V
T
I
01
exp
⁡
V
mp
+
I
mp
R
s
a
1
V
T
+
1
a
2
V
T
I
02
exp
⁡
V
mp
+
I
mp
R
s
a
2
V
T
.Now, using (25) and (27) it is possible to derive from expressions for both saturation currents from (29) and (31)(36)
I
01
=
a
1
a
2
-
a
1
exp
⁡
-
V
oc
a
1
V
T
·
a
2
V
T
R
sh
0
-
R
s
0
-
R
s
0
-
R
s
R
sh
0
I
sc
-
V
oc
R
sh
0
-
R
s
R
s
0
-
R
s
,
(37)
I
02
=
a
2
a
1
-
a
2
exp
⁡
-
V
oc
a
2
V
T
·
a
1
V
T
R
sh
0
-
R
s
0
-
R
s
0
-
R
s
R
sh
0
I
sc
-
V
oc
R
sh
0
-
R
s
R
s
0
-
R
sand, from (32) and (35),(38)
a
1
V
T
=
R
sh
0
I
sc
-
I
mp
-
V
mp
V
mp
+
I
mp
R
s
-
V
oc
a
2
V
T
hhhhh
-
R
sh
0
I
sc
-
V
oc
exp
⁡
V
mp
+
I
mp
R
s
-
V
oc
a
2
V
T
·
R
sh
0
-
V
mp
/
I
mp
V
mp
/
I
mp
-
R
s
h
h
h
h
h
h
h
h
-
R
sh
0
-
R
s
0
R
s
0
-
R
s
exp
⁡
V
mp
+
I
mp
R
s
-
V
oc
a
2
V
T
R
sh
0
-
V
mp
/
I
mp
V
mp
/
I
mp
-
R
s
-
1
,
(39)
a
2
V
T
=
R
sh
0
I
sc
-
I
mp
-
V
mp
V
mp
+
I
mp
R
s
-
V
oc
a
1
V
T
h
h
h
h
h
-
R
sh
0
I
sc
-
V
oc
exp
⁡
V
mp
+
I
mp
R
s
-
V
oc
a
1
V
T
·
R
sh
0
-
V
mp
/
I
mp
V
mp
/
I
mp
-
R
s
h
h
h
h
h
h
h
h
-
R
sh
0
-
R
s
0
R
s
0
-
R
s
exp
⁡
V
mp
+
I
mp
R
s
-
V
oc
a
1
V
T
R
sh
0
-
V
mp
/
I
mp
V
mp
/
I
mp
-
R
s
-
1
.Finally, an implicit expression for the series resistor,R
s, can be obtained from the above equations:(40)
exp
⁡
V
mp
+
I
mp
R
s
-
V
oc
a
1
V
T
=
R
sh
0
I
sc
-
I
mp
-
V
mp
R
sh
0
-
V
mp
/
I
mp
V
mp
/
I
mp
-
R
s
h
h
h
h
h
h
-
a
2
V
T
R
sh
0
-
V
mp
/
I
mp
V
mp
/
I
mp
-
R
s
·
R
sh
0
I
sc
-
V
oc
-
a
2
V
T
R
sh
0
-
R
s
0
R
s
0
-
R
s
-
1which, if a
2 is estimated (a
2
~
2 can be considered as starting point [13]), leads to(41)
log
⁡
R
sh
0
I
sc
-
I
mp
-
V
mp
-
a
2
V
T
R
sh
0
-
V
mp
/
I
mp
/
V
mp
/
I
mp
-
R
s
R
sh
0
I
sc
-
V
oc
-
a
2
V
T
R
sh
0
-
R
s
0
/
R
s
0
-
R
s
-
V
mp
+
I
mp
R
s
-
V
oc
R
sh
0
-
V
mp
/
I
mp
/
V
mp
/
I
mp
-
R
s
-
R
sh
0
-
R
s
0
/
R
s
0
-
R
s
exp
⁡
V
mp
+
I
mp
R
s
-
V
oc
/
a
2
V
T
R
sh
0
I
sc
-
I
mp
-
V
mp
-
R
sh
0
I
sc
-
V
oc
exp
⁡
V
mp
+
I
mp
R
s
-
V
oc
/
a
2
V
T
=
0
.Therefore, the process to extract the seven parameters of (23) can be summarized as follows.(1)
Estimate the parametera
2.
(2)
ObtainR
s from (41).
(3)
Obtain the parametera
1 from (38).
(4)
ObtainI
01 from (36).
(5)
ObtainI
02 from (37).
(6)
ObtainR
sh from (28).
(7)
ObtainI
pv from (25).


### 2.3. Results and Discussion

In Table1 the results of the parameter extraction in relation to the 1-diode/2-resistor and 2-diode/2-resistor equivalent circuit models of the Emcore ZTJ solar cell (see Figure 2) are included. The results, that is, the I-V curves of both models, are compared to the experimental data in Figure 3. With regard to the 1-diode/2-resistor equivalent circuit model, results from the two methods described in Section 2.1 to obtain the series resistor, R
s, have been included in the figure. The first one (labeled “A-1” in Table 1 and Figure 3) is based on a previous estimation of the ideality factor a = 3.6 (i.e., three times a
=
1.2, as this is a three-layer solar cell, and this value of the ideality factor is widely accepted for single-junction cells [60]), and the use of the Lambert W-function. The second one (labeled “A-2” in Table 1 and Figure 3) is based on the slopes of the I-V curve at short circuit and open circuit points. As it can be appreciated, the results fit quite well the testing data, with a maximum deviation of 0.022 A in the case of the 1-diode/2-resistor model and 0.010 A in the case of the 2-diode/2-resistor model.Table 1
Parameters of the 1-diode/2-resistor (A-1 and A-2 methods) and 2-diode/2-resistor equivalent circuit models to the Emcore ZTJ triple-junction (InGaP/InGaAs/Ge) solar cell.
1-diode/2-resistor (A-1)
1-diode/2-resistor (A-2)
2-diode/2-resistor
a
3.6
a
3.6307
a
1
0.9980
I
0 [A]
9.6742
·
10
-
1
4
I
0 [A]
1.2385
·
10
-
13
I
01 [A]
2.4409
·
10
-
16
a
2
2
I
02 [A]
8.4461
·
10
-
10
R
s [Ω]
0.0178
R
s [Ω]
0.0134
R
s [Ω]
0.0802
R
sh [Ω]
337.1056
R
sh [Ω]
343.6025
R
sh [Ω]
343.5357
I
pv [A]
0.4629
I
pv [A]
0.4629
I
pv [A]
0.4629
RMSE
4.43258
·
10
-
3
RMSE
4.4814
·
10
-
3
RMSE
2.8206
·
10
-
3
ξ
9.3469
·
10
-
3
ξ
9.6833
·
10
-
3
ξ
6.0946
·
10
-
3Figure 3
Current-voltage (I-V) curves from the 1-diode/2-resistor (A-1 and A-2 methods) and 2-diode/2-resistor equivalent circuit models to the Emcore ZTJ triple-junction (InGaP/InGaAs/Ge) solar cell (left). The experimental data corresponding to this cell, measured at AM0 (135.3 mW/cm²) and 28°C, has been also included in the graph. Differences between the currents from the studied models and the testing data are included in the right graph.In order to quantify the differences between both models the RMSE may be a better option:(42)
RMSE
=
1
m
∑
i
=
1
m
I
cal
,
i
-
I
i
2as it has been used extensively in the literature [28, 37, 43, 44, 49, 61–72]. In the above expression I
cal
,
i stands for the values calculated using the proposed method, whereas I
i stands for the values from the testing data. Nevertheless, using RMSE does not allow the direct comparison among different solar cells technologies, as the levels of current are different between them. The nondimensional parameter ξ has been proposed to carry out such comparisons [13]:(43)
ξ
=
RMSE
I
sc
=
1
I
sc
1
m
∑
i
=
1
m
I
cal
,
i
-
I
i
2
.Values of RMSE andξ are included in Table 1 for the analyzed methods. The 2-diode/2-resistor model represents a more accurate approach to the data. Regarding the A-1 and A-2 methods followed in relation with the 1-diode/2-resistor model, no big differences are shown with regard to the comparison parameters RMSE and ξ. Slight differences, however, exist when the parameters of the circuit are examined. Taking into account the small difference between both approximations, the first one was selected in the calculations and models carried out to study the power subsystem of the UPMSat-2 satellite [73].Finally, it can be said that, based on the values of the nondimensional parameterξ, the results (i.e., the accuracy) showed by the analytical approximations proposed are similar to the ones from different numerical fittings proposed in the available literature (in [13] such a comparison was carried out based on well-known experimental data from a silicon cell [68], and numerical approximations were characterized by values from ξ = 2.2 · 10−3 to ξ = 6.16 · 10−3).


## 2.1.  The 1-Diode/2-Resistor Electric Circuit Model

The equation that describes the behavior of photovoltaic devices following this approximation is(1)
I
=
I
pv
-
I
0
exp
⁡
V
+
I
R
s
a
V
T
-
1
-
V
+
I
R
s
R
sh
,where I
pv is the photovoltaic current, R
s is the series resistance, R
sh is the shunt resistance, I
0 is the saturation current of the diode, a is the ideality factor of the diode, and V
T is the thermal voltage of the cell, defined as(2)
V
T
=
k
T
q
,where T is the temperature, k is the Boltzmann constant, and q is the charge of the electron.According to (1), five parameters (I
pv, I
0, R
s, R
sh, and a) must be identified before using the model to calculate the performance of the studied photovoltaic device. The analytical approach is normally based on data from the most representative points (short circuit: V
=
0, I
=
I
sc; open circuit: V
=
V
oc, I
=
0; and maximum power: V
=
V
mp, I
=
I
mp points) of the measured I-V curve of the solar cell/panel (as an example, see in Figure 2 the current-voltage and power-voltage curves of a triple-junction solar cell (Emcore ZTJ), developed in recent years for space applications [54–58]. The characteristic points of the I-V curve are I
sc
=
0.462 A; V
oc
=
2.73 V; and V
mp
=
2.41 V and I
mp
=
0.441 A [56]). The aforementioned three points help to derive four equations which allow identifying the value of four of the aforementioned parameters in relation to the fifth (normally, this fifth parameter is the ideality factor, a).Figure 2
Emcore ZTJ triple-junction (InGaP/InGaAs/Ge) solar cell current-voltage (I-V) and power-voltage curves (solid and dashed lines, resp.), measured at AM0 (135.3 mW/cm²) and 28°C.With data (current and voltage levels) from the short circuit and the open circuit points of theI-V curve and taking into account the magnitude of each one of the terms, it is possible to derive the following equations [13]:(3)
I
pv
=
R
sh
+
R
s
R
sh
I
sc
,
(4)
I
0
=
R
sh
+
R
s
I
sc
-
V
oc
R
sh
exp
⁡
V
oc
/
a
V
T
,whereas from the maximum power an additional equation is obtained(5)
I
mp
=
I
sc
-
I
sc
-
V
oc
-
R
s
I
sc
R
sh
exp
⁡
V
mp
+
I
mp
R
s
-
V
oc
a
V
T
-
V
mp
+
I
mp
R
s
-
I
sc
R
s
R
sh
.Finally, taking into account the condition of the maximum power at the maximum power point(6)
∂
P
∂
V
=
V
∂
I
∂
V
+
I
=
0an implicit expression of the series resistor, R
s, as a function of the initial parameters can be obtained:(7)
a
V
T
V
mp
2
I
mp
-
I
sc
·
V
mp
I
sc
+
V
oc
I
mp
-
I
sc
h
h
h
h
h
h
·
V
mp
-
I
mp
R
s
-
a
V
T
V
mp
I
sc
-
V
oc
I
mp
-
1
=
exp
⁡
V
mp
+
I
mp
R
s
-
V
oc
a
V
Ttogether with another for the shunt resistor, R
sh, as a function of R
s and the initial parameters:(8)
R
sh
=
V
mp
-
I
mp
R
s
V
mp
-
R
s
I
sc
-
I
mp
-
a
V
T
V
mp
-
I
mp
R
s
I
sc
-
I
mp
-
a
V
T
I
mp
.As aforementioned, parametera needs to be estimated in order to solve the above equations, and one of them is unfortunately implicit for R
s. One solution to this difficulty (leaving aside any numerical help) consists in using the Lambert W-function, W
(
z
), which is defined as(9)
z
=
W
z
e
W
(
z
)
,where z is any complex number [59]. This function is not injective; in the real variable x the relation W is defined only for x
≥
-
1
/
e and is double valued in the bracket [
-
1
/
e
,
0
], the two branches of the function being expressed as W
0
(
x
) for W
(
x
)
≥
-
1 and W
-
1
(
x
) for W
(
x
)
≤
-
1 in the aforementioned bracket. The general strategy to apply the Lambert W-function in solving exponential equations is to use the following equivalence:(10)
X
=
Y
e
Y
⟺
Y
=
W
X
.Equation (7) can be rewritten as follows:(11)
-
V
mp
2
I
mp
-
I
sc
V
mp
I
sc
+
V
oc
I
mp
-
I
sc
·
exp
⁡
⁡
-
2
V
mp
-
V
oc
a
V
T
+
V
mp
I
sc
-
V
oc
I
mp
V
mp
I
sc
+
V
oc
I
mp
-
I
sc
=
I
mp
R
s
-
V
mp
a
V
T
+
V
mp
I
sc
-
V
oc
I
mp
V
mp
I
sc
+
V
oc
I
mp
-
I
sc
·
exp
⁡
⁡
I
mp
R
s
-
V
mp
a
V
T
+
V
mp
I
sc
-
V
oc
I
mp
V
mp
I
sc
+
V
oc
I
mp
-
I
scand then, using the mathematical equivalence from (10), the above equation turns into(12)
I
mp
R
s
-
V
mp
a
V
T
+
V
mp
I
sc
-
V
oc
I
mp
V
mp
I
sc
+
V
oc
I
mp
-
I
sc
=
W
-
1
-
V
mp
2
I
mp
-
I
sc
V
mp
I
sc
+
V
oc
I
mp
-
I
sc
h
h
h
h
h
h
h
h
·
exp
⁡
V
mp
I
sc
-
V
oc
I
mp
V
mp
I
sc
+
V
oc
I
mp
-
I
sc
-
2
V
mp
-
V
oc
a
V
T
h
h
h
h
h
h
h
h
h
h
h
h
.
h
+
V
mp
I
sc
-
V
oc
I
mp
V
mp
I
sc
+
V
oc
I
mp
-
I
sc
,where W
-
1 is the negative branch of the Lambert W-function (as the left part of (14) is lower than −1 for typical cells and solar panels). And then an explicit expression for R
s is obtained:(13)
R
s
=
A
W
-
1
B
exp
⁡
⁡
C
-
D
+
C
,where (14)
A
=
a
V
T
I
mp
,
B
=
-
V
mp
2
I
mp
-
I
sc
V
mp
I
sc
+
V
oc
I
mp
-
I
sc
,
C
=
-
2
V
mp
-
V
oc
a
V
T
+
V
mp
I
sc
-
V
oc
I
mp
V
mp
I
sc
+
V
oc
I
mp
-
I
sc
,
D
=
V
mp
-
V
oc
a
V
T
.An alternative to the estimation of the ideality factor,a, consists in using other specific characteristics of the I-V curve, such as the slope at short circuit and open circuit points, R
sh
0 and R
s
0, respectively,(15)
R
sh
0
=
-
∂
V
∂
I
I
=
I
sc
,
(16)
R
s
0
=
-
∂
V
∂
I
V
=
V
octo calculate the ideality factor, a. It is preferable to use (15) rather than (16), as the final derived expressions are less complicated. Once R
sh
0 has been obtained, (1) at the short circuit point of the I-V curve can be expressed as [49](17)
1
R
sh
-
1
R
sh
0
-
R
s
+
I
0
a
V
T
exp
⁡
R
s
I
sc
a
V
T
=
0which, taking into account a magnitude analysis of the three terms (in [13] the order of magnitude of terms from the above equation was compared for different technologies of solar cells; R
sh
-
1 and R
sh
0
-
R
s
-
1 had very similar values, whereas the term I
0
/
a
V
T
exp
⁡
R
s
I
sc
/
a
V
T was, in the best case, around 60 times lower), can be simplified as(18)
1
R
sh
-
1
R
sh
0
-
R
s
=
0
;that is,(19)
R
sh
0
=
R
sh
+
R
s
.Once the above expression has been reached, it is possible to derive an expression that gives the ideality factor as a function of constants already known,R
s and the slope R
sh
0:(20)
a
V
T
=
V
mp
-
I
mp
R
s
V
mp
+
I
mp
-
I
sc
R
sh
0
V
mp
-
I
mp
R
sh
0
.Finally, from (7) and (20) the value of the series resistor, R
s, is obtained:(21)
R
s
=
E
-
F
E
+
F
V
mp
I
mp
+
F
E
+
F
V
oc
I
mp
,where(22)
E
=
V
mp
+
I
mp
-
I
sc
R
sh
0
log
⁡
V
mp
+
I
mp
-
I
sc
R
sh
0
V
oc
-
I
sc
R
sh
0
,
F
=
V
mp
-
R
sh
0
I
mp
.As said in Section1, this 1-diode/2-resistor electric circuit model has been successfully correlated to solar cells behavior [13–16].


## 2.2. The 2-Diode/2-Resistor Electric Circuit Model

As this equivalent circuit has one more diode, the model equation contains two exponential terms:(23)
I
=
I
pv
-
I
01
exp
⁡
V
+
I
R
s
a
1
V
T
-
1
-
I
02
exp
⁡
V
+
I
R
s
a
2
V
T
-
1
-
V
+
I
R
s
R
sh
;therefore, four new parameters, two corresponding to different saturation currents, I
01 and I
02, together with other two corresponding to the ideality factors of each diode, a
1 and a
2, are introduced in the calculations instead of the single saturation current and ideality factor from the 1-diode/2-resistor model equation. To identify the parameters of the above equation, the three characteristic points of the I-V curve are considered as boundary conditions, together with the slope of the curve at two of them. As a consequence, six equations are derived to identify seven parameters and, therefore, one of these parameters needs to be estimated.Operating as in the previous subsection, the short circuit point of theI-V curve leads to the following expression: (24)
I
sc
=
I
pv
-
I
01
exp
⁡
I
sc
R
s
a
1
V
T
-
1
-
I
02
exp
⁡
I
sc
R
s
a
2
V
T
-
1
-
I
sc
R
s
R
shwhich, taking into account the much lower magnitude of the terms involving the exponentials, can be simplified as the same relationship (3) obtained with the 1-diode/2-resistor model:(25)
I
pv
=
R
sh
+
R
s
R
sh
I
sc
.Besides, the derivative of (23) at the short circuit point can be rewritten as(26)
1
+
R
sh
+
R
s
(
d
I
/
d
V
)
sc
R
sh
1
+
R
s
(
d
I
/
d
V
)
sc
=
-
1
a
1
V
T
I
01
exp
⁡
I
sc
R
s
a
1
V
T
-
1
a
2
V
T
I
02
exp
⁡
I
sc
R
s
a
2
V
Twhich, bearing in mind the slope of the I-V curve (i.e., (15)), leads to (27)
R
sh
0
-
R
sh
+
R
s
R
sh
R
sh
0
-
R
s
=
-
1
a
1
V
T
I
01
exp
⁡
I
sc
R
s
a
1
V
T
-
1
a
2
V
T
I
02
exp
⁡
I
sc
R
s
a
2
V
T
.Again, comparing the magnitudes of the different terms in the above equation, the already known equation (19) can be obtained:(28)
R
sh
0
=
R
sh
+
R
s
.From the open circuit point, the following expression can be obtained:(29)
R
sh
+
R
s
I
sc
-
V
oc
R
sh
=
I
01
exp
⁡
V
oc
a
1
V
T
+
I
02
exp
⁡
V
oc
a
2
V
T
.If the derivative of (23) at that point(30)
1
+
R
sh
+
R
s
d
I
/
d
V
oc
R
sh
1
+
R
s
d
I
/
d
V
oc
=
-
1
a
1
V
T
I
01
exp
⁡
V
oc
a
1
V
T
-
1
a
2
V
T
I
02
exp
⁡
V
oc
a
2
V
Tis considered, together with the slope of the I-V (i.e., (16)), it is possible to derive the following equation:(31)
R
sh
+
R
s
-
R
s
0
R
sh
R
s
0
-
R
s
=
1
a
1
V
T
I
01
exp
⁡
V
oc
a
1
V
T
+
1
a
2
V
T
I
02
exp
⁡
V
oc
a
2
V
T
.From the maximum power point data the equation(32)
V
mp
+
I
mp
-
I
sc
R
sh
+
R
s
R
sh
=
-
I
01
exp
⁡
V
mp
+
I
mp
R
s
a
1
V
T
-
I
02
exp
⁡
V
mp
+
I
mp
R
s
a
2
V
Tis obtained. Furthermore, bearing in mind the maximum power condition (6)(33)
∂
I
∂
V
I
mp
,
V
mp
=
-
I
mp
V
mpit is possible to derive the last expression needed:(34)
1
+
R
sh
+
R
s
d
I
/
d
V
mp
R
sh
1
+
R
s
d
I
/
d
V
mp
=
-
1
a
1
V
T
I
01
exp
⁡
V
mp
+
I
mp
R
s
a
1
V
T
-
1
a
2
V
T
I
02
exp
⁡
V
mp
+
I
mp
R
s
a
2
V
T
,
(35)
R
sh
+
R
s
I
mp
-
V
mp
R
sh
V
mp
-
R
s
I
mp
=
1
a
1
V
T
I
01
exp
⁡
V
mp
+
I
mp
R
s
a
1
V
T
+
1
a
2
V
T
I
02
exp
⁡
V
mp
+
I
mp
R
s
a
2
V
T
.Now, using (25) and (27) it is possible to derive from expressions for both saturation currents from (29) and (31)(36)
I
01
=
a
1
a
2
-
a
1
exp
⁡
-
V
oc
a
1
V
T
·
a
2
V
T
R
sh
0
-
R
s
0
-
R
s
0
-
R
s
R
sh
0
I
sc
-
V
oc
R
sh
0
-
R
s
R
s
0
-
R
s
,
(37)
I
02
=
a
2
a
1
-
a
2
exp
⁡
-
V
oc
a
2
V
T
·
a
1
V
T
R
sh
0
-
R
s
0
-
R
s
0
-
R
s
R
sh
0
I
sc
-
V
oc
R
sh
0
-
R
s
R
s
0
-
R
sand, from (32) and (35),(38)
a
1
V
T
=
R
sh
0
I
sc
-
I
mp
-
V
mp
V
mp
+
I
mp
R
s
-
V
oc
a
2
V
T
hhhhh
-
R
sh
0
I
sc
-
V
oc
exp
⁡
V
mp
+
I
mp
R
s
-
V
oc
a
2
V
T
·
R
sh
0
-
V
mp
/
I
mp
V
mp
/
I
mp
-
R
s
h
h
h
h
h
h
h
h
-
R
sh
0
-
R
s
0
R
s
0
-
R
s
exp
⁡
V
mp
+
I
mp
R
s
-
V
oc
a
2
V
T
R
sh
0
-
V
mp
/
I
mp
V
mp
/
I
mp
-
R
s
-
1
,
(39)
a
2
V
T
=
R
sh
0
I
sc
-
I
mp
-
V
mp
V
mp
+
I
mp
R
s
-
V
oc
a
1
V
T
h
h
h
h
h
-
R
sh
0
I
sc
-
V
oc
exp
⁡
V
mp
+
I
mp
R
s
-
V
oc
a
1
V
T
·
R
sh
0
-
V
mp
/
I
mp
V
mp
/
I
mp
-
R
s
h
h
h
h
h
h
h
h
-
R
sh
0
-
R
s
0
R
s
0
-
R
s
exp
⁡
V
mp
+
I
mp
R
s
-
V
oc
a
1
V
T
R
sh
0
-
V
mp
/
I
mp
V
mp
/
I
mp
-
R
s
-
1
.Finally, an implicit expression for the series resistor,R
s, can be obtained from the above equations:(40)
exp
⁡
V
mp
+
I
mp
R
s
-
V
oc
a
1
V
T
=
R
sh
0
I
sc
-
I
mp
-
V
mp
R
sh
0
-
V
mp
/
I
mp
V
mp
/
I
mp
-
R
s
h
h
h
h
h
h
-
a
2
V
T
R
sh
0
-
V
mp
/
I
mp
V
mp
/
I
mp
-
R
s
·
R
sh
0
I
sc
-
V
oc
-
a
2
V
T
R
sh
0
-
R
s
0
R
s
0
-
R
s
-
1which, if a
2 is estimated (a
2
~
2 can be considered as starting point [13]), leads to(41)
log
⁡
R
sh
0
I
sc
-
I
mp
-
V
mp
-
a
2
V
T
R
sh
0
-
V
mp
/
I
mp
/
V
mp
/
I
mp
-
R
s
R
sh
0
I
sc
-
V
oc
-
a
2
V
T
R
sh
0
-
R
s
0
/
R
s
0
-
R
s
-
V
mp
+
I
mp
R
s
-
V
oc
R
sh
0
-
V
mp
/
I
mp
/
V
mp
/
I
mp
-
R
s
-
R
sh
0
-
R
s
0
/
R
s
0
-
R
s
exp
⁡
V
mp
+
I
mp
R
s
-
V
oc
/
a
2
V
T
R
sh
0
I
sc
-
I
mp
-
V
mp
-
R
sh
0
I
sc
-
V
oc
exp
⁡
V
mp
+
I
mp
R
s
-
V
oc
/
a
2
V
T
=
0
.Therefore, the process to extract the seven parameters of (23) can be summarized as follows.(1)
Estimate the parametera
2.
(2)
ObtainR
s from (41).
(3)
Obtain the parametera
1 from (38).
(4)
ObtainI
01 from (36).
(5)
ObtainI
02 from (37).
(6)
ObtainR
sh from (28).
(7)
ObtainI
pv from (25).


## 2.3. Results and Discussion

In Table1 the results of the parameter extraction in relation to the 1-diode/2-resistor and 2-diode/2-resistor equivalent circuit models of the Emcore ZTJ solar cell (see Figure 2) are included. The results, that is, the I-V curves of both models, are compared to the experimental data in Figure 3. With regard to the 1-diode/2-resistor equivalent circuit model, results from the two methods described in Section 2.1 to obtain the series resistor, R
s, have been included in the figure. The first one (labeled “A-1” in Table 1 and Figure 3) is based on a previous estimation of the ideality factor a = 3.6 (i.e., three times a
=
1.2, as this is a three-layer solar cell, and this value of the ideality factor is widely accepted for single-junction cells [60]), and the use of the Lambert W-function. The second one (labeled “A-2” in Table 1 and Figure 3) is based on the slopes of the I-V curve at short circuit and open circuit points. As it can be appreciated, the results fit quite well the testing data, with a maximum deviation of 0.022 A in the case of the 1-diode/2-resistor model and 0.010 A in the case of the 2-diode/2-resistor model.Table 1
Parameters of the 1-diode/2-resistor (A-1 and A-2 methods) and 2-diode/2-resistor equivalent circuit models to the Emcore ZTJ triple-junction (InGaP/InGaAs/Ge) solar cell.
1-diode/2-resistor (A-1)
1-diode/2-resistor (A-2)
2-diode/2-resistor
a
3.6
a
3.6307
a
1
0.9980
I
0 [A]
9.6742
·
10
-
1
4
I
0 [A]
1.2385
·
10
-
13
I
01 [A]
2.4409
·
10
-
16
a
2
2
I
02 [A]
8.4461
·
10
-
10
R
s [Ω]
0.0178
R
s [Ω]
0.0134
R
s [Ω]
0.0802
R
sh [Ω]
337.1056
R
sh [Ω]
343.6025
R
sh [Ω]
343.5357
I
pv [A]
0.4629
I
pv [A]
0.4629
I
pv [A]
0.4629
RMSE
4.43258
·
10
-
3
RMSE
4.4814
·
10
-
3
RMSE
2.8206
·
10
-
3
ξ
9.3469
·
10
-
3
ξ
9.6833
·
10
-
3
ξ
6.0946
·
10
-
3Figure 3
Current-voltage (I-V) curves from the 1-diode/2-resistor (A-1 and A-2 methods) and 2-diode/2-resistor equivalent circuit models to the Emcore ZTJ triple-junction (InGaP/InGaAs/Ge) solar cell (left). The experimental data corresponding to this cell, measured at AM0 (135.3 mW/cm²) and 28°C, has been also included in the graph. Differences between the currents from the studied models and the testing data are included in the right graph.In order to quantify the differences between both models the RMSE may be a better option:(42)
RMSE
=
1
m
∑
i
=
1
m
I
cal
,
i
-
I
i
2as it has been used extensively in the literature [28, 37, 43, 44, 49, 61–72]. In the above expression I
cal
,
i stands for the values calculated using the proposed method, whereas I
i stands for the values from the testing data. Nevertheless, using RMSE does not allow the direct comparison among different solar cells technologies, as the levels of current are different between them. The nondimensional parameter ξ has been proposed to carry out such comparisons [13]:(43)
ξ
=
RMSE
I
sc
=
1
I
sc
1
m
∑
i
=
1
m
I
cal
,
i
-
I
i
2
.Values of RMSE andξ are included in Table 1 for the analyzed methods. The 2-diode/2-resistor model represents a more accurate approach to the data. Regarding the A-1 and A-2 methods followed in relation with the 1-diode/2-resistor model, no big differences are shown with regard to the comparison parameters RMSE and ξ. Slight differences, however, exist when the parameters of the circuit are examined. Taking into account the small difference between both approximations, the first one was selected in the calculations and models carried out to study the power subsystem of the UPMSat-2 satellite [73].Finally, it can be said that, based on the values of the nondimensional parameterξ, the results (i.e., the accuracy) showed by the analytical approximations proposed are similar to the ones from different numerical fittings proposed in the available literature (in [13] such a comparison was carried out based on well-known experimental data from a silicon cell [68], and numerical approximations were characterized by values from ξ = 2.2 · 10−3 to ξ = 6.16 · 10−3).


## 3. Cup Anemometer Analytical Model Used to Study the “Flag” Position of Damaged Anemometers

The performance of the cup anemometer rotor can be analyzed with the following expression:(44)
I
z
d
ω
d
t
=
Q
A
+
Q
f
,where I
z is the moment of inertia of the rotor, Q
A is the aerodynamic torque, and Q
f is the frictional torque, which depends on the air temperature, T, and the rotation speed, ω. If the effect of the friction is left out from the above equation (as its effect is much lower than the aerodynamic torque in standard conditions) and the three cups are considered, the following equation can be derived:(45)
I
z
d
ω
d
t
=
1
2
ρ
S
c
R
rc
V
r
2
θ
c
N
α
θ
+
1
2
ρ
S
c
R
rc
V
r
2
θ
+
120
°
c
N
α
θ
+
120
°
+
1
2
ρ
S
c
R
rc
V
r
2
θ
+
240
°
c
N
α
θ
+
240
°
,where ρ is the air density, V
r is the wind speed relative to the cups, c
N is the aerodynamic normal force coefficient, α is the local wind direction with respect to the cups, θ is the angle of the rotor with respect to a reference line (see Figure 4), S
c is the front area of the cups (S
c
=
π
R
c
2, R
c being the cup radius), and R
rc is the cups’ center rotation radius. Wind speed V
r, relative to the cup at rotor angle θ with respect to the reference line, is expressed as(46)
V
r
θ
=
V
2
+
ω
R
rc
2
-
2
V
ω
R
rc
cos
⁡
θ
,the wind direction with respect to the cup, α, being derived, as a function of the rotor’s position angle, θ, from the following expression [22, 23]:(47)
tan
α
=
K
sin
θ
K
cos
⁡
θ
-
1
,where K is the anemometer factor, defined by expression (48)
K
=
V
ω
R
rc
.Figure 4
Sketch of a cup anemometer.These variables are sketched in Figure5, together with the wind tunnel measurement of the c
N coefficient regarding a conical cup, which are plotted in relation to angle α.Figure 5
Wind tunnel measurements of the normal aerodynamic force coefficient,c
N, of a conical cup, plotted as a function of the wind angle, α, regarding the cup. These measurements were carried out with the cup in affixed position (i.e., not considering any rotation).From (45) it is possible to calculate the anemometer performance on steady wind and the effect of perturbations (i.e., the time constant) [23, 26]. It also can be used to analyze the equilibrium or “flag” positions that sometimes damaged cup anemometers can present. Furthermore, a damaged anemometer staying in one of these particular “flag” positions can give false measurements of the wind speed thanks to the tiny oscillations around this position, which activate the optoelectronic output signal systems used in most of class 1 anemometers.


### 3.1. Results and Discussion

In Figure6 a damaged anemometer (one cup missing) is shown, together with its calibration curves. It can be observed that the damaged rotor presents a very linear behavior and even better aerodynamic characteristics, that is, more rotational frequency when compared to the nondamaged situation. Leaving aside the anomaly detection problem that this effect represents [19], the frequency measured by the cup anemometer can be observed in the figure when it stays at a “flag” position. The behavior (i.e., the rotational frequency) can be considered linearly dependent on the wind speed, which makes this state difficult to be identified if there are no means of comparison with other anemometer measurements.Figure 6
Cup anemometer with damaged rotor, that is, with one cup missing (a). Calibration curves of the anemometer equipped with damaged and nondamaged rotors (b). The rotation frequency measured by the anemometer at “flag” position is also included in the graph.
(a)
(b)If the “flag” position of a damaged anemometer with one cup missing is considered, (45) should be rewritten as(49)
I
z
d
ω
d
t
=
1
2
ρ
S
c
R
rc
V
r
2
θ
c
N
θ
+
1
2
ρ
S
c
R
rc
V
r
2
θ
+
120
°
c
N
θ
+
120
°
.Obviously, at that equilibrium positionα
=
θ, as the cups do not rotate. Then, if small fluctuations around the “flag” position are taken into account:(50)
θ
=
θ
0
+
γ
γ
≪
1together with the approximation [23](51)
c
N
θ
=
c
0
+
c
1
cos
⁡
θ
,the following relationships can be derived, taking only into account first-order terms:(52)
c
N
θ
=
c
0
+
c
1
cos
⁡
θ
0
-
c
1
γ
sin
θ
0
,
V
r
2
=
V
2
-
2
V
d
γ
d
t
R
rc
cos
⁡
θ
0
.Therefore, (49) can be rewritten as(53)
I
z
(
1
/
2
)
ρ
S
c
R
rc
d
2
γ
d
t
2
=
V
2
-
2
V
d
γ
d
t
R
rc
cos
⁡
θ
0
·
c
0
+
c
1
cos
⁡
θ
0
-
c
1
γ
sin
θ
0
+
V
2
-
2
V
d
γ
d
t
R
rc
cos
⁡
θ
0
+
120
°
·
c
0
+
c
1
cos
⁡
θ
0
+
120
°
-
c
1
γ
sin
θ
0
+
120
°
.Now it is possible to classify the different terms into two groups, one non-time-dependent, which gives the positions of the equilibrium state (i.e., the values ofθ
0)(54)
0
=
2
c
0
+
c
1
cos
⁡
θ
0
+
c
1
cos
⁡
θ
0
+
120
°
,and a second one, which allows analyzing the nonsteady state:(55)
I
z
(
1
/
2
)
ρ
S
c
R
rc
V
2
d
2
γ
d
t
2
=
-
c
1
γ
sin
θ
0
+
sin
θ
0
+
120
°
-
2
R
rc
V
d
γ
d
t
cos
⁡
2
⁡
c
0
cos
⁡
θ
0
+
cos
⁡
θ
0
+
120
°
hhhhhhhhhhhhh
+
c
1
cos
⁡
2
⁡
θ
0
+
cos
⁡
2
⁡
θ
0
+
120
°
.Going back to (54), it can be simplified as(56)
sin
θ
0
-
30
°
=
2
c
0
c
1which gives two solutions. Conical cups have a coefficient ratio around c
0
/
c
1
=
0.302; therefore, θ
0
=
67
° and θ
0
=
172
°. See in Figure 7 the cup anemometer with one cup missing at “flag” position, with one of the cups around the first solution (θ
0
=
67
°). It was checked that the equilibrium state was reached at a position around the first analytical solution obtained. However, the second solution was never obtained during the testing campaign. Even more, it was found to be unstable in all attempts to force the anemometer to remain in that position.Figure 7
Damaged cup anemometer with one cup missing at “flag” position during the testing campaign (wind tunnel wind speed is aroundV = 6 m/s; wind-flow direction is indicated with an arrow in the picture). View from the top-window of the calibration tunnel at the IDR/UPM Institute.In order to check the stability of the solutions obtained it is necessary to analyze (55), which can be simplified taking into account (56) and using the reduced time variable τ
=
t
V
/
R
rc, as(57)
I
z
(
1
/
2
)
ρ
S
c
R
rc
3
c
1
d
2
γ
d
τ
2
+
3
1
-
sin
θ
0
-
30
°
d
γ
d
τ
+
cos
⁡
θ
0
-
30
°
γ
=
0which, taking into account that the last term is negative for θ
0
=
172
°, demonstrates that this solution is not stable (as it has a root with positive real part).


## 3.1. Results and Discussion

In Figure6 a damaged anemometer (one cup missing) is shown, together with its calibration curves. It can be observed that the damaged rotor presents a very linear behavior and even better aerodynamic characteristics, that is, more rotational frequency when compared to the nondamaged situation. Leaving aside the anomaly detection problem that this effect represents [19], the frequency measured by the cup anemometer can be observed in the figure when it stays at a “flag” position. The behavior (i.e., the rotational frequency) can be considered linearly dependent on the wind speed, which makes this state difficult to be identified if there are no means of comparison with other anemometer measurements.Figure 6
Cup anemometer with damaged rotor, that is, with one cup missing (a). Calibration curves of the anemometer equipped with damaged and nondamaged rotors (b). The rotation frequency measured by the anemometer at “flag” position is also included in the graph.
(a)
(b)If the “flag” position of a damaged anemometer with one cup missing is considered, (45) should be rewritten as(49)
I
z
d
ω
d
t
=
1
2
ρ
S
c
R
rc
V
r
2
θ
c
N
θ
+
1
2
ρ
S
c
R
rc
V
r
2
θ
+
120
°
c
N
θ
+
120
°
.Obviously, at that equilibrium positionα
=
θ, as the cups do not rotate. Then, if small fluctuations around the “flag” position are taken into account:(50)
θ
=
θ
0
+
γ
γ
≪
1together with the approximation [23](51)
c
N
θ
=
c
0
+
c
1
cos
⁡
θ
,the following relationships can be derived, taking only into account first-order terms:(52)
c
N
θ
=
c
0
+
c
1
cos
⁡
θ
0
-
c
1
γ
sin
θ
0
,
V
r
2
=
V
2
-
2
V
d
γ
d
t
R
rc
cos
⁡
θ
0
.Therefore, (49) can be rewritten as(53)
I
z
(
1
/
2
)
ρ
S
c
R
rc
d
2
γ
d
t
2
=
V
2
-
2
V
d
γ
d
t
R
rc
cos
⁡
θ
0
·
c
0
+
c
1
cos
⁡
θ
0
-
c
1
γ
sin
θ
0
+
V
2
-
2
V
d
γ
d
t
R
rc
cos
⁡
θ
0
+
120
°
·
c
0
+
c
1
cos
⁡
θ
0
+
120
°
-
c
1
γ
sin
θ
0
+
120
°
.Now it is possible to classify the different terms into two groups, one non-time-dependent, which gives the positions of the equilibrium state (i.e., the values ofθ
0)(54)
0
=
2
c
0
+
c
1
cos
⁡
θ
0
+
c
1
cos
⁡
θ
0
+
120
°
,and a second one, which allows analyzing the nonsteady state:(55)
I
z
(
1
/
2
)
ρ
S
c
R
rc
V
2
d
2
γ
d
t
2
=
-
c
1
γ
sin
θ
0
+
sin
θ
0
+
120
°
-
2
R
rc
V
d
γ
d
t
cos
⁡
2
⁡
c
0
cos
⁡
θ
0
+
cos
⁡
θ
0
+
120
°
hhhhhhhhhhhhh
+
c
1
cos
⁡
2
⁡
θ
0
+
cos
⁡
2
⁡
θ
0
+
120
°
.Going back to (54), it can be simplified as(56)
sin
θ
0
-
30
°
=
2
c
0
c
1which gives two solutions. Conical cups have a coefficient ratio around c
0
/
c
1
=
0.302; therefore, θ
0
=
67
° and θ
0
=
172
°. See in Figure 7 the cup anemometer with one cup missing at “flag” position, with one of the cups around the first solution (θ
0
=
67
°). It was checked that the equilibrium state was reached at a position around the first analytical solution obtained. However, the second solution was never obtained during the testing campaign. Even more, it was found to be unstable in all attempts to force the anemometer to remain in that position.Figure 7
Damaged cup anemometer with one cup missing at “flag” position during the testing campaign (wind tunnel wind speed is aroundV = 6 m/s; wind-flow direction is indicated with an arrow in the picture). View from the top-window of the calibration tunnel at the IDR/UPM Institute.In order to check the stability of the solutions obtained it is necessary to analyze (55), which can be simplified taking into account (56) and using the reduced time variable τ
=
t
V
/
R
rc, as(57)
I
z
(
1
/
2
)
ρ
S
c
R
rc
3
c
1
d
2
γ
d
τ
2
+
3
1
-
sin
θ
0
-
30
°
d
γ
d
τ
+
cos
⁡
θ
0
-
30
°
γ
=
0which, taking into account that the last term is negative for θ
0
=
172
°, demonstrates that this solution is not stable (as it has a root with positive real part).


## 4. Modeling the Effect of Pressure Waves in High-Speed Railway Tunnels

A well-known aerodynamic phenomenon in high-speed trains is the propagation of the pressure waves generated when a train passes through a tunnel, which can give rise to several problems. The reflections of these waves at the discontinuities of the tunnel can meet the train several times inside the tunnel and can cause discomfort to the train passengers and also to the pedestrians near the tunnel entrance/exit due to the micropressure radiated outside from the tunnel portals. When a train approaches a tunnel entrance portal, the air inside the tunnel begins to be compressed due to the open-air local field preceding the train nose. Once the train head reaches the tunnel entrance the pressure increases faster until the nose is completely inside the tunnel. This interaction between the tunnel entrance hood and the train nose determines the shape of the compression wave generated, but the total pressure gaugeΔ
p
N remains constant. The pressure continues to increase slowly due to the viscous effects on the tunnel wall and the train surface Δ
p
f
r, until the end of the train reaches the tunnel entrance portal, and a rarefaction wave is then generated with a pressure decrease Δ
p
T. The pressure decreases suddenly Δ
p
H
P when the train nose passes the measurement point. These pressure variations are known as the train reference signature and are defined in EN14067-5:2006:E standard [74] (see Figure 8).Figure 8
Train reference signature defined in EN14067-5:2006:E standard [74].The design of railway tunnels for high-speed train operation requires accurate models to study the wave propagation process and its interaction with the vehicle which implies long calculation times. However, this process can be somehow alleviated by using simplified models. The main advantage of simplified models is that the computational cost is quite low compared with more complex models. The most common simplified approach to deal with the interaction of pressure waves and vehicle is the use of models based on the one-dimensional equations of gas dynamics.A method based on the propagation of the signature of the pressure waves inside the tunnel, which requires low CPU cost, is presented in [75]. This method of characteristics, 1D finite difference, is used as a numerical scheme, only predicts pressure changes, and is only valid for simple tunnels (without abrupt cross-section changes, air shafts, etc.), but allows investigation of trains crossings using different time delays in the entry time of the trains. In the zones where the flow is strongly three-dimensional (around tunnel portals, train ends) other methods have to be used.Concerning the interaction of the train head with the tunnel portal, an analytical method based on the aeroacoustic theory developed by Howe is applied to a short “acoustically compact” hood having a window on its side wall in [76]. The boundary element method is used to calculate the potential flow through a portal of arbitrary geometry. A parametric analysis and an experimental study are performed in order to show the great influence in its performance of the relative length of the train nose to the distance between the window and the hood entrance.The flow generated by a train in a tunnel with side branches is studied in [77]. The results confirm that the use of shafts provides a very effective method to reduce pressure fluctuations in tunnels. The shafts divide the waves into secondary waves which can be more effectively damped by viscous dissipation or acoustic diffusion.The pressure waves are reflected at tunnel discontinuities such as changes in cross-sectional area, airshafts, side branches, and also at the train’s head (or tail) when they hit its surface. These multiple reflections can create localized regions where constructive interferences produce higher pressure that can also affect the passengers and personnel in tunnel surroundings.


### 4.1. Mathematical Model for Prediction of Pressure Waves inside a Tunnel

The unsteady aerodynamic effects in railway tunnels, like pressure waves generated during the train passing through the tunnel, can be analyzed by Euler equations, which are particularly useful in preliminary design work, where information about pressure alone is desired with low CPU cost [78].The assumptions to obtain the one-dimensional wave equation are [79](1)
very small density fluctuations (linear acoustics),
(2)
inviscid fluid,
(3)
adiabatic process (this is quite reasonable except in the tunnel portal entrance/exit and near the walls),
(4)
fluid in local thermodynamic equilibrium (Stokes’ hypothesis) (this assumption may partially fail at high frequencies resulting in a dissipation related to the volume changes [80]),
(5)
homogeneous fluid,
(6)
stagnant and uniform fluid (quiescent fluid),
(7)
air as an ideal gas,
(8)
irrotational flow (this is quite reasonable except in the tunnel portal entrance/exit and in the train wake).With all the above considerations and in absence of sources, the continuity equation is(58)
∂
ρ
′
∂
t
+
ρ
0
∇
·
u
→
′
=
0
,where ρ
0 is the mean density in a stationary inviscid fluid and ρ
′, u
→
′ are the density and speed perturbations. Euler’s equation (momentum conservation law for a frictionless fluid in absence of external forces) is(59)
ρ
0
∂
u
→
′
∂
t
+
∇
p
′
=
0
,where p
′ is the pressure perturbation. Using the constitutive equation (p
′
=
c
0
2
ρ
′, as flow is isentropic) we obtain the wave equation for p
′
(60)
1
c
0
2
∂
2
p
′
∂
t
2
-
∇
2
p
′
=
0
.The general equation that approximates the time-harmonic acoustic waves propagation inside a hard-wall hollow cylindrical duct is(61)
1
r
∂
p
′
∂
r
+
∂
2
p
′
∂
r
2
+
1
r
2
∂
2
p
′
∂
φ
2
+
∂
2
p
′
∂
x
2
-
1
c
2
∂
2
p
′
∂
t
2
=
0
,where the instantaneous pressure inside the tube depends on the distance to the sound source x, the distance to the longitudinal axis of the duct r, the angle of the point position with the horizontal axis φ, the speed of sound in the medium c, and the time t. The solution of the homogeneous equation is (taking into account the boundary conditions (
∂
p
/
∂
r
)
r
=
a
=
0 and p
′
x
=
L
=
0)(62)
p
ω
′
(
x
,
r
,
φ
,
t
)
=
∑
m
=
-
∞
m
=
+
∞
∑
n
=
0
n
=
+
∞
c
m
n
J
m
q
m
n
r
r
i
e
i
k
m
n
x
+
m
φ
-
ω
t
,where ω is the angular frequency, k
m
n
=
±
(
ω
/
c
)
2
-
(
q
m
n
/
r
i
)
2 wave number, and r
i cylinder radius, J
m are Bessel functions of first kind, and q
m
n are the zeros of these Bessel functions J
m (the eigenvalues). The coefficients c
m
n are obtained from the boundary condition (p
′
(
0
,
r
,
φ
,
t
0
)
=
F
w
(
r
,
φ
)
e
-
i
ω
t), applying the orthogonality properties of Bessel functions:(63)
c
m
n
=
1
π
a
2
q
m
n
2
q
m
n
2
-
m
2
J
m
2
q
m
n
r
r
i
·
∫
0
2
π
∫
0
a
F
ω
r
,
φ
J
m
q
m
n
r
r
i
e
-
i
m
φ
r
d
r
d
φ
.


### 4.2. Condition for the Existence of a Plane Wave

The propagation mode(
m
,
n
) of an acoustic disturbance that propagates axially along the tube depends on the wave number. The first propagation mode of the disturbance is a plane wave (of wave number k
00
=
ω
/
c) in which there are not traversal fluctuations (
q
m
n
=
0
). The condition of propagation of an acoustic mode is that the wave number k
m
n must be real; otherwise the wave will decay exponentially, and this is known as an evanescent wave. The lowest zero of the Bessel function J
1 is q
11
=
1.8412; therefore the cut-off frequency is f
c
=
1.8412
c
/
(
2
π
a
). In a railway tunnel (a
~
5 m) for high-speed trains this frequency is close to 20 Hz.The wave generated at the tunnel entrance by the train can be considered as a compact source, which means that the tunnel cross-section characteristic length(
L
=
2
a
) of variation of the sources is small compared with the wave length λ, and the flow can locally be approximated as an incompressible potential flow. This can be done when the head of the wave is larger compared with the tunnel diameter. At the same time, the region of emission, in this case the portal, can also be considered compact. If a source region is compact compared to the sound wavelength it can be represented by the sum of basic sources (monopole, dipole, and quadrupole). For example, the train can be described as a compact source, as moving excitation (F
(
x
,
t
)
=
δ
(
x
-
U
tr
⁡
,
t
)), where δ is a Dirac function, and a Green function could be used to describe the pressure field.The pressure waves generated by a train entering inside a tunnel can be treated as plane waves after several diameters distance away from the tunnel portals (where the flow is strongly three-dimensional) [80].


### 4.3. Results: Propagation of the Signature

The train wave signature (TWS) method proposed in [75] has been used to predict the pressure wave during the train pass inside the tunnel. It is based on one-dimensional Green’s function and the results are validated with experimental measurements. The problem is similar to a piston moving in a finite pipe, but with both sides open because the blockage, ϕ
=
A
tr
⁡
/
A
tun (i.e., ratio between the cross-section area of the train, A
tr
⁡, and the cross-section area of the tunnel A
tun), is much smaller than unity and the waves that find the train experience a small reflection and practically all the energy is transmitted (pressure waves see the train as a small cross-section area reduction). The hypotheses assumed in this model are summarized as follows, in these groups, concerning the train, the tunnel, and the fluid.For the train, the hypotheses are(1)
constant speed,
(2)
slenderness enough for the flow not to be detached behind the nose,
(3)
leakage through the train walls modeled through a sealing constant (values vary from 0 (train window open) up to 20 seconds (in highly tight coaches of modern high-speed trains)).For the tunnel, the hypotheses are(1)
axisymmetric, circular cylinder duct,
(2)
hard walls (no transmission of sound through these walls),
(3)
constant perimeter and cross-sections area,
(4)
the altitude capable of varying with the distance to the entrance,
(5)
blockage being small(
ϕ
=
1
).For the fluid, the hypotheses are as follows.(1)
Pressure waves are assumed to be plane waves. We are considering a pressure pulse whose wavefronts (surfaces of constant phase) are parallel planes of constant peak-to-peak amplitude normal to the phase velocity vector.
(2)
Air behaves as an ideal gas, with isentropic, quiescent fluid (u
0
=
0). This hypothesis is modified in nonlinear propagation.The acoustic field induced by the near field of the train (TNS) is considered constant and does not generate any reflections with the wave that hits the train’s head or tail when passing through the tunnel, nor with cross-section area changes (or airshafts) along the aforementioned tunnel.The principle of the TWS method is to propagate the pressure profiles (TNS and TWS and its reflections at the tunnel portals). At each point of the tunnel the sum of the effects due to the passage of the pressure profiles is calculated to obtain the resulting pressure at this location. The amplitude of the reflections of the TWS at the tunnel portals is determined from experimental results.The one-dimensional d’Alembert’s solution for wave equation inx-axis direction is(64)
p
′
=
f
x
-
c
0
t
+
g
x
+
c
0
t
,
u
′
=
1
ρ
0
c
0
f
x
-
c
0
t
-
g
x
+
c
0
t
,where f and g are determined by boundary and initial conditions. They can be interpreted as longitudinal waves (fluid velocity is parallel to the direction of propagation) propagating in positive x direction f and g in opposite direction. The particle velocity induced by an acoustic plane wave in a stagnant medium is then u
′
=
p
′
/
ρ
0
c
0, where ρ
0
c
0 is the characteristic impedance of the fluid.The expressions used to describe the wave signatures have been extracted from [74, 81]. The expression used to estimate Δ
p
N is(65)
Δ
p
N
=
ρ
U
tr
⁡
2
1
-
M
2
ϕ
1
+
ϕ
,where U
tr
⁡ is the speed of the train, M is the Mach number based on the speed of the train, M
=
U
tr
⁡
/
c
0, and ϕ is the blockage. The total aerodynamic drag of a train in a tunnel is usually expressed [10] in terms of a nose loss coefficient, k
N, which implies a stagnation pressure loss as Δ
p
s
nose
=
1
/
2
k
N
ρ
U
nose
2 and a tail loss coefficient, k
T, which implies a stagnation pressure loss as Δ
p
s
tail
=
1
/
2
k
N
ρ
U
tail
2. Note that the velocities in the annulus between the train and the tunnel wall will not, in general, be the same at the nose and tail for a long train. An estimation of these coefficients when k
N
~
0.1 or less, for a streamlined nose, is given in [82].The expression used for the pressure gauge generated during the entry of the train body is [74](66)
Δ
p
f
r
=
γ
P
0
M
X
4
1
+
2
L
D
h
,
tun
c
f
tun
1
-
M
X
4
-
Δ
p
Nwith(67)
X
1
=
4
L
tr
⁡
D
h
,
tun
c
f
tun
1
M
-
1
-
1
-
ϕ
3
c
f
tr
⁡
ϕ
+
c
f
tun
-
ξ
n
1
1
-
ϕ
2
-
1
,
X
2
=
1
M
+
4
L
tr
⁡
D
h
,
tun
c
f
tun
1
M
-
1
-
1
1
-
ϕ
2
,
X
3
=
2
M
+
4
L
tr
⁡
D
h
,
tun
c
f
tun
1
M
-
1
-
1
1
-
ϕ
,
X
4
=
1
+
-
X
2
+
X
2
2
-
X
1
X
3
X
1and ξ
n
=
1
+
ξ
h
/
(
1
-
1
-
ϕ
2
). The coefficient ξ
n is a stagnation pressure loss coefficient at the train nose (common values between 0.03 and 0.05), and ξ
h is the stagnation pressure loss coefficient at the train in the tunnel. D
h
,
tun is the hydraulic inner diameter of the tunnel, P
tun the tunnel perimeter, and A
tun the area section of the tunnel and the friction coefficients of the train c
f
tr
⁡ (common values for high-speed train are between 0.003 and 0.005) and of the tunnel c
f
tun (values between 0.0050 and 0.0065).The expression used for the pressure drop due to the entrance of the train rear, according to [74], is(68)
Δ
p
T
=
-
Δ
p
N
1
-
ξ
n
1
-
1
-
ϕ
2
.Since the TWS is generated during the train entrance, the characteristic time is the same as that of the TNS. The Mach number based on the speed of the train indicates the ratio between the different traveling speed of the perturbation TWS and the train near-field TNS. The characteristic wavelength of the TWS,L
TWS, is therefore a combination of the train length, L
tr
⁡, the train-tunnel entrance, L
e, and the speed of sound:(69)
L
TWS
≃
L
tr
⁡
+
L
e
M
,where L
e is the length of the entrance. The entrance can modify the rise time of the pressure; this is a common method used to reduce the pressure gradient and with it the micropressure radiated at the tunnel portals.The attenuation is taken into account by using the relation [83](70)
Δ
p
=
Δ
p
0
e
-
α
a
c
0
t
,where the coefficient α
a is estimated experimentally. It is approximately 1.8 · 10−4 for a tunnel with ballasted track and around 3 · 10−5 for slab track tunnels.As the pressure at the tunnel exit must remain constant, the boundary condition is(71)
F
1
(
ξ
)
+
G
1
(
η
)
=
0
,where ξ
=
L
tun
-
c
0
t and η
=
L
tun
+
c
0
t. Thus the reflected wave is(72)
G
1
(
η
)
=
-
F
1
(
ξ
)
=
-
F
1
(
2
L
tun
-
η
)
=
-
f
2
L
tun
-
x
+
c
0
tand, taking into account the reflection coefficient at the exit R
exit from the experimental measurements, we can use the expression(73)
G
1
η
=
-
R
exit
·
f
2
L
tun
-
x
-
c
0
t
.For the reflection at the entrance (x
=
0, ξ
=
-
c
0
t, and η
=
-
ξ) we obtain(74)
F
2
ξ
=
-
G
1
(
η
)
=
-
G
1
(
-
ξ
)
=
R
entrance
·
f
2
L
tun
+
x
-
c
0
t
.A recurrent expression for the wave generated by the train entrance and the reflection at the exits would be considering exponential attenuation (to simplify the expressions, it is assumedR
exit
=
R
entrance
=
R, but, in the model, experimental results could be introduced by the geometric influence of the portal in the wave reflection):(75)
F
x
,
t
=
∑
s
g
f
k
L
L
tun
+
s
g
x
-
c
0
t
L
TWS
e
-
α
a
c
0
t
R
i
n
-
1
,where s
g
=
sign
(
-
1
)
i
n
+
1, i
n is the wave counter, and the wavenumber k
L
=
2
floor
(
i
n
/
2
) (as the function floor of Matlab calculates the nearest integer, the series of k
L is 0, 2, 2, 4, 4, 6, 6). For a wave generated at the exit, for example, the wave generated when the train leaves the tunnel, the general expression will be(76)
G
x
,
t
=
∑
-
s
g
1
g
k
L
1
L
tun
+
L
tun
c
0
/
U
tr
⁡
+
s
g
1
x
-
c
0
t
L
TWS
·
e
-
α
a
c
0
t
-
t
1
R
i
n
-
1
,where s
g
1
=
sign
(
-
1
)
i
n, the wavenumber k
L
1
=
2
floor
(
(
i
n
-
1
)
/
2
)
+
1 (the series of k
L
1 is then 1, 1, 3, 3, 5, 5, 7, 7), and t
1 is the instant in which the train reaches the exit (t
1
=
L
tun
/
U
tr
⁡). The TWS results from [75] with the TWS implementation with the selected expressions for Δ
p
N, Δ
p
f
r, and Δ
p
T are presented in Figure 9.Figure 9
Pressure evolution with time in the tunnel at 550 m from the tunnel entrance. Validation of the TWS implementation,M
=
0.24 and φ
=
0.135. Experimental results from [75].The pressure inside the train is calculated with the dynamic sealing constantτ
dyn (which represent the time for the train inside pressure p
i, to become a fraction of the pressure outside p
e):(77)
d
p
i
d
t
=
1
τ
dyn
p
e
-
p
i
.The pressure in the train exterior at 72 m from the nose is shown in Figure10. There are some situations in which the pressure gauge inside the train or in the tunnel reaches critical values or variations. Some of these cases are also mentioned in EN14067-5:2006:E standard [74].Figure 10
Pressure evolution with time in the train at 72 m from the nose. Validation of the TWS implementation,M
=
0.24 and φ
=
0.135. Experimental results from [75].The characteristic times of the situation of interest, denoted as points 1, 2, and 3 in Figure11, are(1)
t
c
t, when the end of the TWS reaches the train head (point 1)(78)
t
c
t
=
L
tr
⁡
U
tr
⁡
1
1
-
M
,
(2)
t
c
i, when the TWS reflected at the tunnel exit portal finds the train nose (point 2):(79)
t
c
i
=
2
L
tun
c
0
1
+
M
.Figure 11
Positionx of the pressure waves generated by a high-speed train inside a tunnel of length L
tun as a function of time t. C: compression; E: expansion wave.The TWS generated during the train exit (a compression wave) finds another compression wave (point 3). This is the situation in which the pressure inside the tunnel can reach its maximum; this situation is not as important for the train that just left the tunnel but for another train that can enter the tunnel just after the first one. A study of the time delays between the trains entering the tunnel must be done to design the minimum tunnel cross-section area. The critical situations depend mainly on the time delay between train entrances, train lengths, train speeds, tunnel lengths, reflection coefficients of the pressure waves at the portals, attenuation along the tunnel, existence of side branches, airshaft, and so forth.The implemented code can be used to quickly and approximately predict the induced particle velocity inside the tunnel and find the optimal locations for the damper devices based on aerodynamics instabilities, as a function of the train-tunnel operation parameters (train and tunnel lengths, train speed, sound speed, reflection coefficients at the portals, etc.).


## 4.1. Mathematical Model for Prediction of Pressure Waves inside a Tunnel

The unsteady aerodynamic effects in railway tunnels, like pressure waves generated during the train passing through the tunnel, can be analyzed by Euler equations, which are particularly useful in preliminary design work, where information about pressure alone is desired with low CPU cost [78].The assumptions to obtain the one-dimensional wave equation are [79](1)
very small density fluctuations (linear acoustics),
(2)
inviscid fluid,
(3)
adiabatic process (this is quite reasonable except in the tunnel portal entrance/exit and near the walls),
(4)
fluid in local thermodynamic equilibrium (Stokes’ hypothesis) (this assumption may partially fail at high frequencies resulting in a dissipation related to the volume changes [80]),
(5)
homogeneous fluid,
(6)
stagnant and uniform fluid (quiescent fluid),
(7)
air as an ideal gas,
(8)
irrotational flow (this is quite reasonable except in the tunnel portal entrance/exit and in the train wake).With all the above considerations and in absence of sources, the continuity equation is(58)
∂
ρ
′
∂
t
+
ρ
0
∇
·
u
→
′
=
0
,where ρ
0 is the mean density in a stationary inviscid fluid and ρ
′, u
→
′ are the density and speed perturbations. Euler’s equation (momentum conservation law for a frictionless fluid in absence of external forces) is(59)
ρ
0
∂
u
→
′
∂
t
+
∇
p
′
=
0
,where p
′ is the pressure perturbation. Using the constitutive equation (p
′
=
c
0
2
ρ
′, as flow is isentropic) we obtain the wave equation for p
′
(60)
1
c
0
2
∂
2
p
′
∂
t
2
-
∇
2
p
′
=
0
.The general equation that approximates the time-harmonic acoustic waves propagation inside a hard-wall hollow cylindrical duct is(61)
1
r
∂
p
′
∂
r
+
∂
2
p
′
∂
r
2
+
1
r
2
∂
2
p
′
∂
φ
2
+
∂
2
p
′
∂
x
2
-
1
c
2
∂
2
p
′
∂
t
2
=
0
,where the instantaneous pressure inside the tube depends on the distance to the sound source x, the distance to the longitudinal axis of the duct r, the angle of the point position with the horizontal axis φ, the speed of sound in the medium c, and the time t. The solution of the homogeneous equation is (taking into account the boundary conditions (
∂
p
/
∂
r
)
r
=
a
=
0 and p
′
x
=
L
=
0)(62)
p
ω
′
(
x
,
r
,
φ
,
t
)
=
∑
m
=
-
∞
m
=
+
∞
∑
n
=
0
n
=
+
∞
c
m
n
J
m
q
m
n
r
r
i
e
i
k
m
n
x
+
m
φ
-
ω
t
,where ω is the angular frequency, k
m
n
=
±
(
ω
/
c
)
2
-
(
q
m
n
/
r
i
)
2 wave number, and r
i cylinder radius, J
m are Bessel functions of first kind, and q
m
n are the zeros of these Bessel functions J
m (the eigenvalues). The coefficients c
m
n are obtained from the boundary condition (p
′
(
0
,
r
,
φ
,
t
0
)
=
F
w
(
r
,
φ
)
e
-
i
ω
t), applying the orthogonality properties of Bessel functions:(63)
c
m
n
=
1
π
a
2
q
m
n
2
q
m
n
2
-
m
2
J
m
2
q
m
n
r
r
i
·
∫
0
2
π
∫
0
a
F
ω
r
,
φ
J
m
q
m
n
r
r
i
e
-
i
m
φ
r
d
r
d
φ
.


## 4.2. Condition for the Existence of a Plane Wave

The propagation mode(
m
,
n
) of an acoustic disturbance that propagates axially along the tube depends on the wave number. The first propagation mode of the disturbance is a plane wave (of wave number k
00
=
ω
/
c) in which there are not traversal fluctuations (
q
m
n
=
0
). The condition of propagation of an acoustic mode is that the wave number k
m
n must be real; otherwise the wave will decay exponentially, and this is known as an evanescent wave. The lowest zero of the Bessel function J
1 is q
11
=
1.8412; therefore the cut-off frequency is f
c
=
1.8412
c
/
(
2
π
a
). In a railway tunnel (a
~
5 m) for high-speed trains this frequency is close to 20 Hz.The wave generated at the tunnel entrance by the train can be considered as a compact source, which means that the tunnel cross-section characteristic length(
L
=
2
a
) of variation of the sources is small compared with the wave length λ, and the flow can locally be approximated as an incompressible potential flow. This can be done when the head of the wave is larger compared with the tunnel diameter. At the same time, the region of emission, in this case the portal, can also be considered compact. If a source region is compact compared to the sound wavelength it can be represented by the sum of basic sources (monopole, dipole, and quadrupole). For example, the train can be described as a compact source, as moving excitation (F
(
x
,
t
)
=
δ
(
x
-
U
tr
⁡
,
t
)), where δ is a Dirac function, and a Green function could be used to describe the pressure field.The pressure waves generated by a train entering inside a tunnel can be treated as plane waves after several diameters distance away from the tunnel portals (where the flow is strongly three-dimensional) [80].


## 4.3. Results: Propagation of the Signature

The train wave signature (TWS) method proposed in [75] has been used to predict the pressure wave during the train pass inside the tunnel. It is based on one-dimensional Green’s function and the results are validated with experimental measurements. The problem is similar to a piston moving in a finite pipe, but with both sides open because the blockage, ϕ
=
A
tr
⁡
/
A
tun (i.e., ratio between the cross-section area of the train, A
tr
⁡, and the cross-section area of the tunnel A
tun), is much smaller than unity and the waves that find the train experience a small reflection and practically all the energy is transmitted (pressure waves see the train as a small cross-section area reduction). The hypotheses assumed in this model are summarized as follows, in these groups, concerning the train, the tunnel, and the fluid.For the train, the hypotheses are(1)
constant speed,
(2)
slenderness enough for the flow not to be detached behind the nose,
(3)
leakage through the train walls modeled through a sealing constant (values vary from 0 (train window open) up to 20 seconds (in highly tight coaches of modern high-speed trains)).For the tunnel, the hypotheses are(1)
axisymmetric, circular cylinder duct,
(2)
hard walls (no transmission of sound through these walls),
(3)
constant perimeter and cross-sections area,
(4)
the altitude capable of varying with the distance to the entrance,
(5)
blockage being small(
ϕ
=
1
).For the fluid, the hypotheses are as follows.(1)
Pressure waves are assumed to be plane waves. We are considering a pressure pulse whose wavefronts (surfaces of constant phase) are parallel planes of constant peak-to-peak amplitude normal to the phase velocity vector.
(2)
Air behaves as an ideal gas, with isentropic, quiescent fluid (u
0
=
0). This hypothesis is modified in nonlinear propagation.The acoustic field induced by the near field of the train (TNS) is considered constant and does not generate any reflections with the wave that hits the train’s head or tail when passing through the tunnel, nor with cross-section area changes (or airshafts) along the aforementioned tunnel.The principle of the TWS method is to propagate the pressure profiles (TNS and TWS and its reflections at the tunnel portals). At each point of the tunnel the sum of the effects due to the passage of the pressure profiles is calculated to obtain the resulting pressure at this location. The amplitude of the reflections of the TWS at the tunnel portals is determined from experimental results.The one-dimensional d’Alembert’s solution for wave equation inx-axis direction is(64)
p
′
=
f
x
-
c
0
t
+
g
x
+
c
0
t
,
u
′
=
1
ρ
0
c
0
f
x
-
c
0
t
-
g
x
+
c
0
t
,where f and g are determined by boundary and initial conditions. They can be interpreted as longitudinal waves (fluid velocity is parallel to the direction of propagation) propagating in positive x direction f and g in opposite direction. The particle velocity induced by an acoustic plane wave in a stagnant medium is then u
′
=
p
′
/
ρ
0
c
0, where ρ
0
c
0 is the characteristic impedance of the fluid.The expressions used to describe the wave signatures have been extracted from [74, 81]. The expression used to estimate Δ
p
N is(65)
Δ
p
N
=
ρ
U
tr
⁡
2
1
-
M
2
ϕ
1
+
ϕ
,where U
tr
⁡ is the speed of the train, M is the Mach number based on the speed of the train, M
=
U
tr
⁡
/
c
0, and ϕ is the blockage. The total aerodynamic drag of a train in a tunnel is usually expressed [10] in terms of a nose loss coefficient, k
N, which implies a stagnation pressure loss as Δ
p
s
nose
=
1
/
2
k
N
ρ
U
nose
2 and a tail loss coefficient, k
T, which implies a stagnation pressure loss as Δ
p
s
tail
=
1
/
2
k
N
ρ
U
tail
2. Note that the velocities in the annulus between the train and the tunnel wall will not, in general, be the same at the nose and tail for a long train. An estimation of these coefficients when k
N
~
0.1 or less, for a streamlined nose, is given in [82].The expression used for the pressure gauge generated during the entry of the train body is [74](66)
Δ
p
f
r
=
γ
P
0
M
X
4
1
+
2
L
D
h
,
tun
c
f
tun
1
-
M
X
4
-
Δ
p
Nwith(67)
X
1
=
4
L
tr
⁡
D
h
,
tun
c
f
tun
1
M
-
1
-
1
-
ϕ
3
c
f
tr
⁡
ϕ
+
c
f
tun
-
ξ
n
1
1
-
ϕ
2
-
1
,
X
2
=
1
M
+
4
L
tr
⁡
D
h
,
tun
c
f
tun
1
M
-
1
-
1
1
-
ϕ
2
,
X
3
=
2
M
+
4
L
tr
⁡
D
h
,
tun
c
f
tun
1
M
-
1
-
1
1
-
ϕ
,
X
4
=
1
+
-
X
2
+
X
2
2
-
X
1
X
3
X
1and ξ
n
=
1
+
ξ
h
/
(
1
-
1
-
ϕ
2
). The coefficient ξ
n is a stagnation pressure loss coefficient at the train nose (common values between 0.03 and 0.05), and ξ
h is the stagnation pressure loss coefficient at the train in the tunnel. D
h
,
tun is the hydraulic inner diameter of the tunnel, P
tun the tunnel perimeter, and A
tun the area section of the tunnel and the friction coefficients of the train c
f
tr
⁡ (common values for high-speed train are between 0.003 and 0.005) and of the tunnel c
f
tun (values between 0.0050 and 0.0065).The expression used for the pressure drop due to the entrance of the train rear, according to [74], is(68)
Δ
p
T
=
-
Δ
p
N
1
-
ξ
n
1
-
1
-
ϕ
2
.Since the TWS is generated during the train entrance, the characteristic time is the same as that of the TNS. The Mach number based on the speed of the train indicates the ratio between the different traveling speed of the perturbation TWS and the train near-field TNS. The characteristic wavelength of the TWS,L
TWS, is therefore a combination of the train length, L
tr
⁡, the train-tunnel entrance, L
e, and the speed of sound:(69)
L
TWS
≃
L
tr
⁡
+
L
e
M
,where L
e is the length of the entrance. The entrance can modify the rise time of the pressure; this is a common method used to reduce the pressure gradient and with it the micropressure radiated at the tunnel portals.The attenuation is taken into account by using the relation [83](70)
Δ
p
=
Δ
p
0
e
-
α
a
c
0
t
,where the coefficient α
a is estimated experimentally. It is approximately 1.8 · 10−4 for a tunnel with ballasted track and around 3 · 10−5 for slab track tunnels.As the pressure at the tunnel exit must remain constant, the boundary condition is(71)
F
1
(
ξ
)
+
G
1
(
η
)
=
0
,where ξ
=
L
tun
-
c
0
t and η
=
L
tun
+
c
0
t. Thus the reflected wave is(72)
G
1
(
η
)
=
-
F
1
(
ξ
)
=
-
F
1
(
2
L
tun
-
η
)
=
-
f
2
L
tun
-
x
+
c
0
tand, taking into account the reflection coefficient at the exit R
exit from the experimental measurements, we can use the expression(73)
G
1
η
=
-
R
exit
·
f
2
L
tun
-
x
-
c
0
t
.For the reflection at the entrance (x
=
0, ξ
=
-
c
0
t, and η
=
-
ξ) we obtain(74)
F
2
ξ
=
-
G
1
(
η
)
=
-
G
1
(
-
ξ
)
=
R
entrance
·
f
2
L
tun
+
x
-
c
0
t
.A recurrent expression for the wave generated by the train entrance and the reflection at the exits would be considering exponential attenuation (to simplify the expressions, it is assumedR
exit
=
R
entrance
=
R, but, in the model, experimental results could be introduced by the geometric influence of the portal in the wave reflection):(75)
F
x
,
t
=
∑
s
g
f
k
L
L
tun
+
s
g
x
-
c
0
t
L
TWS
e
-
α
a
c
0
t
R
i
n
-
1
,where s
g
=
sign
(
-
1
)
i
n
+
1, i
n is the wave counter, and the wavenumber k
L
=
2
floor
(
i
n
/
2
) (as the function floor of Matlab calculates the nearest integer, the series of k
L is 0, 2, 2, 4, 4, 6, 6). For a wave generated at the exit, for example, the wave generated when the train leaves the tunnel, the general expression will be(76)
G
x
,
t
=
∑
-
s
g
1
g
k
L
1
L
tun
+
L
tun
c
0
/
U
tr
⁡
+
s
g
1
x
-
c
0
t
L
TWS
·
e
-
α
a
c
0
t
-
t
1
R
i
n
-
1
,where s
g
1
=
sign
(
-
1
)
i
n, the wavenumber k
L
1
=
2
floor
(
(
i
n
-
1
)
/
2
)
+
1 (the series of k
L
1 is then 1, 1, 3, 3, 5, 5, 7, 7), and t
1 is the instant in which the train reaches the exit (t
1
=
L
tun
/
U
tr
⁡). The TWS results from [75] with the TWS implementation with the selected expressions for Δ
p
N, Δ
p
f
r, and Δ
p
T are presented in Figure 9.Figure 9
Pressure evolution with time in the tunnel at 550 m from the tunnel entrance. Validation of the TWS implementation,M
=
0.24 and φ
=
0.135. Experimental results from [75].The pressure inside the train is calculated with the dynamic sealing constantτ
dyn (which represent the time for the train inside pressure p
i, to become a fraction of the pressure outside p
e):(77)
d
p
i
d
t
=
1
τ
dyn
p
e
-
p
i
.The pressure in the train exterior at 72 m from the nose is shown in Figure10. There are some situations in which the pressure gauge inside the train or in the tunnel reaches critical values or variations. Some of these cases are also mentioned in EN14067-5:2006:E standard [74].Figure 10
Pressure evolution with time in the train at 72 m from the nose. Validation of the TWS implementation,M
=
0.24 and φ
=
0.135. Experimental results from [75].The characteristic times of the situation of interest, denoted as points 1, 2, and 3 in Figure11, are(1)
t
c
t, when the end of the TWS reaches the train head (point 1)(78)
t
c
t
=
L
tr
⁡
U
tr
⁡
1
1
-
M
,
(2)
t
c
i, when the TWS reflected at the tunnel exit portal finds the train nose (point 2):(79)
t
c
i
=
2
L
tun
c
0
1
+
M
.Figure 11
Positionx of the pressure waves generated by a high-speed train inside a tunnel of length L
tun as a function of time t. C: compression; E: expansion wave.The TWS generated during the train exit (a compression wave) finds another compression wave (point 3). This is the situation in which the pressure inside the tunnel can reach its maximum; this situation is not as important for the train that just left the tunnel but for another train that can enter the tunnel just after the first one. A study of the time delays between the trains entering the tunnel must be done to design the minimum tunnel cross-section area. The critical situations depend mainly on the time delay between train entrances, train lengths, train speeds, tunnel lengths, reflection coefficients of the pressure waves at the portals, attenuation along the tunnel, existence of side branches, airshaft, and so forth.The implemented code can be used to quickly and approximately predict the induced particle velocity inside the tunnel and find the optimal locations for the damper devices based on aerodynamics instabilities, as a function of the train-tunnel operation parameters (train and tunnel lengths, train speed, sound speed, reflection coefficients at the portals, etc.).


## 5. Conclusions

In the present paper analytical methods have been applied to study three different engineering problems. The high level of complexity in relation to the problems studied has not been a major obstacle in order to obtain valuable solutions. The methods developed in the present work preserve, to a certain point, the physics of the problems analyzed and can be modified in order to study other specific problems related to photovoltaic systems, cup anemometer performances, and the pressure distribution on high-speed train tunnels.



---

*Source: 897357-2015-04-16.xml*