Mathias Berthod

Commitment and efficiency-inducing tax and subsidy scheme in the development of a clean technology

Introduction

This paper analyses the optimal environmental policy design in situations where the regulator -hereafter a government -can not strongly commit to announcements about future tax and subsidy levels. The motivation is that long-run perspectives of environmental policies often face short run concerns. One consistent illustration in Europe is related to the french government which cancelled the carbon tax increase for year 2019 following recent demonstrations of the yellow vest movement. Other examples include the case of the australian government abolishing the carbon tax in 2014, or the spanish government abruptly cancelling the renewable energy subsidies in 2012. I specifically consider environmental policies which aim at supporting the transition from the use of dirty technologies to clean technologies by subsidizing innovation. The interplay between innovation and environmental policies has been extensively addressed. 1 However, a large share of the literature abstracts from the issue of commitment. In most papers, the analysis consists in comparing the optimal policy and a business-as-usual scenario (see [START_REF] Bosetti | Optimal energy investment and R&D strategies to stabilize atmospheric greenhouse gas concentrations[END_REF]Edenhoffer et al., 2006 ;[START_REF] Popp | ENTICE-BR: The effects of backstop technology R&D on climate policy models[END_REF]). Yet several authors point out that the government lack of commitment may lead to inefficient environmental innovation (see [START_REF] Wirl | Taxes versus permits as incentive for the intertemporal supply of a clean technology by a monopoly[END_REF][START_REF] Montero | A note on environmental policy and innovation when governments can not commit[END_REF].

The question thus arises: if a government can not strongly commit to announcements about future tax and subsidy levels, is there an efficient policy design? And, if so, how does it differ from the case of strong commitment?

To tackle this issue, I analyse a model where two perfectly substitutable technologies are available to produce a final good. The "dirty" technology, used by a competitive industry, is associated with a flow of pollution. The "clean" technology, owned by a monopoly, is nonpolluting but also initially more expensive than the dirty technology. The monopoly can invest to improve its clean technology. The improvement is modeled as the accumulation of a stock of knowledge, which reduces the cost associated with the clean technology as in [START_REF] Tsur | Optimal transition to backstop substitutes for nonrenewable resources[END_REF]. The government can tax the pollution flow and subsidy monopoly investments. In the timeline of the game setting, the government is the first-mover (leader), it sets up its policy before the monopoly, which then chooses its production and investment levels as the secondmover (follower). This paper then characterizes the efficiency-inducing tax and subsidy schemes for two different commitment levels adopted by the government.

The analysis is divided into three steps. First, I describe the social optimum as a benchmark.

Second, I characterize the Open-Loop Stackelberg Equilibrium (OLSE) of the game. This models a strong commitment level: the government announces and fully commits to the future path of the tax and subsidy rates. Third, I characterize the Global Stackelberg Equilibrium (GSE) of the game [START_REF] Basar | Dynamic noncooperative game theory[END_REF]. This models a weaker level of commitment than the OLSE: the government rather commits to a policy rule that characterizes the levels of tax and subsidy rates as a function of the stock of knowledge.

I first show that only a tax is necessary to induce efficiency under strong commitment (OLSE).

Under weak commitment (GSE), I find that the policy rule that induces efficiency also features a subsidy to investment. The subsidy is used to deal with the anticipation by the monopoly that its choices of investment in knowledge will affect the future tax rates. The improvement of the clean technology being modeled via the stock of knowledge, a natural extension of the basic framework is to analyse what happens when two firms can contribute to this stock under technological spillovers. In this case, I show that a tax and a subsidy are necessary to induce efficiency under strong commitment (OLSE). Yet, compared to the case of the monopoly, a subsidy is needed to correct the free-rider problem that emerges from the existence of spillovers.

Under weak commitment (GSE), the firms take into account that any new unit of knowledge decreases the future tax rate but also increases the future subsidy rate. This new effect implies a more complex form of the subsidy rate. A surprising result is that, in some situations, the subsidy rate may be negative in the short run, which implies that it may be optimal to initially tax investments in knowledge.

Regarding the related literature, the present framework is close to [START_REF] Wirl | Taxes versus permits as incentive for the intertemporal supply of a clean technology by a monopoly[END_REF], who studies the effect of government lack of commitment on the intertemporal monopolistic supply of a clean technology. In this paper, the inability of the government to commit leads to a suboptimally low improvement in the clean technology and, thus, to an inefficient allocation of investments. The result is equivalent for emission permits and emission taxes. There are two main differences with regard to my contributions. First of all, Wirl (2013) compares the OLSE to the Markovian Nash equilibrium of the game. As such, he compares what I refer to as a case of strong commitment level to a case where the government can not commit at all. By considering a GSE, I consider an intermediate case where the government commitment level is weak. I can characterize a policy rule that induces efficiency in the absence of a strong commitment, which is the first main contribution. Second, the present framework also differs from [START_REF] Wirl | Taxes versus permits as incentive for the intertemporal supply of a clean technology by a monopoly[END_REF] because he considers production capacity rather than knowledge. The stock variable relates to a physical capital in [START_REF] Wirl | Taxes versus permits as incentive for the intertemporal supply of a clean technology by a monopoly[END_REF]. In this paper, it relates to a stock of knowledge that lowers the cost related to the clean technology. The immediate consequence is that, by contrast to Wirl (2013), the stock is not associated with well-defined property rights in the present model, and it can spill out between firms. The second main contribution of this paper is thus to analyse the effect of a lack of commitment for policy design under technological spillovers and development of an alternative clean technology. Also related to the present framework, [START_REF] Montero | A note on environmental policy and innovation when governments can not commit[END_REF] tackles the issue of environmental policy when the government can not commit. In a two-period model, one potential innovator invests, in the first period, to develop a clean technology that can be used in the second period.

The analysis consists in comparing situations when the government does update (no commitment) or does not update (commitment) its policy in the second period. It shows that the performance of taxes and permits differs in such a setting. The present framework allows to further analyze the innovator's intertemporal investment dynamics in the development of the clean technology.

In Montero (2011), innovator's decision only takes place in the first period.

Ulph and Ulph (2013) focus on a specific source of time-inconsistency of environmental policies, namely the fact that government are in power for a limited period of time and that next governments may give more or less weight to environmental damage costs. In this context, they find that a government may need to use additional R&D subsidies to induce immediate investments in the development of the clean technology. In the present article, I consider a different time-inconsistency, which emerges as a government may be tempted to update its policy and deviate from the originally decided policy path depending on the reaction of the firms. In this different context, I find, as [START_REF] Ulph | Optimal climate change policies when governments cannot commit[END_REF], that a R&D subsidy is needed to make up for the government lack of commitment. However, this is more complex in the presence of multiple firms owning clean technologies under spillover effects. I show that, in some situations, the government may even need to tax firms investments in knowledge in the short run to slow down the pace of development of the clean technology.

In the remainder of this article, I introduce the basic framework and characterize the social optimum in section 2. Then, I characterize the OLSE and the GSE of the Stackelberg dynamic game in section 3. I extend the model to the case of a duopoly under technological spillovers in section 4. Last, I provide some concluding remarks in section 5.

The model and the social optimum

Basic framework

First, in what follows, I always use linear quadratic specifications. This is standard in the literature on differential games [START_REF] Dockner | Differential games in economics and management science[END_REF][START_REF] Engwerda | LQ dynamic optimization and differential games[END_REF]. As is known, this structure yields tractable analytical solutions. Second, I am being parsimonious in the use of parameters to ease the reading. Each parameter is introduced hereafter only if it does matter in the analysis.

Last, the model is suitable to analyse any type of transition from the use of polluting techniques -or technologies -to non-polluting ones for the production of a good or a service. However, throughout the article, I focus on the energy sector to illustrate my framework. 2 For instance, the reader may think about the transition from technologies such as thermal power plants -which 2 Among others, an alternative application of the model could be, for instance, the agricultural sector. Nuclear technologies, which are supported and promoted by the Joint FAO/IAEA Division of Nuclear in Food and Agriculture (see [START_REF] Hendrichs | Improving the cost-effectiveness, trade and safety of biological control for agricultural insect pests using nuclear techniques[END_REF][START_REF] Mabit | Promoting the use of isotopic techniques to combat soil erosion: An overview of the key role played by the SWMCN Subprogramme of the Joint FAO/IAEA Division over the last 20 years[END_REF], are examples of alternative and less polluting technologies.

involve the burning of fossils fuels such as coal or gas -to technologies such as hydro, wind or solar power in the production of electricity.

Consumers enjoy a gross surplus from consuming a homogeneous good q at a price p:3

U (q) = q - q 2 2 (1)
This expression implies a linear inverse demand with a choke price normalized to one:

p = 1 -q (2)
Output q can be produced either using a clean, x, or a dirty technology, y. These two sources of supply are perfect substitutes, so that q = x + y. This modeling of the demand is in line with [START_REF] Wirl | Taxes versus permits as incentive for the intertemporal supply of a clean technology by a monopoly[END_REF]. One may think of thermal power plants using coal as the dirty technology y and of solar energy as the clean technology x, both of them providing electricity q.

Dirty technology y is used on a competitive industry: the marginal cost of production is constant and normalized to zero. Its use induces a flow of pollution, which results in an external cost D(y):

D(y) = γ 2 y 2 , γ > 0 (3)
This constitutes a negative flow externality. For instance, the burning of coal required to operate a thermal power plant results in the emission of sulfur dioxide, a toxic gas which reacts with water molecules in the atmosphere, and causes acid rains.

Clean technology x is owned by a monopoly. This modeling of two different market structures for the technologies -competition for the dirty technology and monopoly for the clean one -is in line with [START_REF] Montero | A note on environmental policy and innovation when governments can not commit[END_REF] and [START_REF] Wirl | Taxes versus permits as incentive for the intertemporal supply of a clean technology by a monopoly[END_REF]. The rationale is that the development of a clean technology proceeds in two steps. In the first step, many firms compete to find a new technology.

In the second step, the owner of the new technology is protected by a patent. In what follows, the analysis is restricted to the second step, and the owner is a monopoly. The cost of using the clean technology depends on a stock of knowledge K:

C x (x, K) = 1 2 x + c x0 -c x1 K x, c x0 > 0, c x1 > 0, c x0 -c x1 K > 0 (4)
This expression implies that the marginal cost associated with the clean technology equals

K = I, K(0) = 0 (5)
This approach, referred to as R&D-induced technological change [START_REF] Popp | Energy, the environment, and technological change[END_REF], is one of the main approaches used to endogenize technological change. An alternative is to consider Learning-by-doing (LBD), where the cumulative experience of using a technology induces cost reductions. My choice to focus on R&D-induced technological change is supported by the fact that, according to empirical studies, R&D activities seem to contribute more than LBD to technology improvements. Using a panel data for wind power installations in four western European countries, Söderholm and Sundqvist (2007) find, a LBD learning rate of around 5%, which is to be compared to 15% for R&D. 4 The absence of depreciation and the normalization of the initial stock of knowledge to zero are assumptions made for simplicity, and do not change the results.

For any investment in the stock of knowledge, the monopoly faces the following adjustment cost:

C I (I, K) = 1 2 I + c I K I (6)
First, this expression is quadratic in the level of investment I, which is a standard assumption in the theory of investment of the firm (see for instance [START_REF] Gould | Adjustment of Costs of in the the Theory Investment Firm[END_REF]). Second, it depends on the stock of knowledge K. This models the fact that improvements of a new technology become increasingly more difficult as the level of R&D advances increases. Indeed, as the stock of accumulated knowledge increases, the marginal adjustment cost of investment rises, which means that it is increasingly costly to increase the stock of knowledge K by a given amount I. This specification is in consistent with the "fishing out effect" introduced by Jones (1995) in growth theory: the easiest ideas are always discovered first. Therefore, as technology becomes more and more complex, it takes more time and more effort to make further improvements. 5Parameter c I in (6) can be interpreted as a measure of the magnitude of the fishing out effect.

To summarize, my framework applies to a situation where a clean technology already exists and is owned by a firm, but it still needs R&D-induced technological change to be improved.

The improvement decreases the cost, which makes the clean technology less expensive relative to the dirty technology. This implies a smooth transition from the dirty technology to the clean one.

Again, this applies well to the energy sector since governments use various environmentally related taxes and subsidy R&D activities in favor of renewable energies. OECD's policy instruments for the environment database (PINE) covers more than 1400 different environmentally related taxes (OECD, 2017). IEA's annual Energy Technology RD&D Budget Database estimates a total budget dedicated to RD&D activities in favor of renewable energy sources to 2,972 million dollars for IEA members (IEA, 2019).

Timeline of the game, equilibrium concepts and players' objectives

To tackle the issue of government commitment level, I consider a game with hierarchical play between a government and the monopoly owning the clean technology. The government (first-mover) can announce its policy before the monopoly (second-mover) makes its decisions. The government policy consists possibly in a tax τ on the dirty technology level and a subsidy σ to investment in the stock of knowledge. I compare two commitment levels from the government, each of them corresponding to a specific equilibrium concept.

(K) = τ 0 + τ 1 K, σ(K) = σ 0 + σ 1 K), K ∈ [0, ∞)}.
This equilibrium concept is in general considered as more realistic, in the sense that the government is not asked to follow a specific path but it rather only observes the evolution of the stock of knowledge and updates the tax rate and the subsidy rate. It models a weak commitment level from the government.

In line with [START_REF] Wirl | Taxes versus permits as incentive for the intertemporal supply of a clean technology by a monopoly[END_REF], the tax τ induces the monopoly to set the following price:

p = τ (7)
Perfect substitution between the clean technology and the dirty technology implies that the monopoly must charge the same price than the competitive industry. Since the cost related to the dirty technology is normalized to zero, the competitive industry profit is py -τ y which yields the expression of the price in [START_REF] Edenhofer | Mitigation strategies and costs of climate protection: The effects of ETC in the hybrid model MIND[END_REF].

At each period, the monopoly instantaneous profits equal the revenue -with a price equal to the tax -plus the subsidy to investment, net of the production cost and the investment adjustment cost:

π = τ x - 1 2 x + c x0 -c x1 K x - 1 2 I + c I K I + σI (8)
Taking into account this announcement, the monopoly chooses production level x and investment level I. The monopoly chooses clean technology level x and investment I to maximize the net present value of its instantaneous profits (8) subject to the accumulation of the stock of knowledge (5):

max x,I ∞ 0 πe -rt dt, subject to K = I (9)
The benevolent government accounts for the profit of the monopoly (π), the net surplus of the consumers (U (q) -pq), the tax revenue (τ y) and the external cost of pollution (D(y)). Using [START_REF] Edenhofer | Mitigation strategies and costs of climate protection: The effects of ETC in the hybrid model MIND[END_REF], I can write that

U (q) = (1 -τ ) -(1-τ ) 2 2 and D(y) = γ 2 y 2 = γ 2 (1 -x -τ ) 2 .
Therefore, the social welfare equals, at each time t:

S = (1 -τ ) - (1 -τ ) 2 2 - 1 2 x + c x0 -c x1 K x - 1 2 I + c I K I - γ 2 (1 -x -τ ) 2 (10)
The government chooses tax rate τ and subsidy rate σ to maximize the aggregate welfare, i.e. the net present value of the stream of social welfare [START_REF] Grimaud | Climate change mitigation options and directed technical change: A decentralized equilibrium analysis[END_REF] subject to the accumulation of the stock of knowledge (5):

7 max τ,σ ∞ 0 Se -rt dt, subject to K = I (11)

The social optimum

Before analyzing the game between the government and the monopoly, I characterize the social optimum as a benchmark. To determine this allocation, I consider a benevolent social planner who chooses tax rate τ , clean technology level x and investment I to maximize the aggregate welfare subject to the evolution of the stock of knowledge:

max τ,x,I ∞ 0 Se -rt dt, subject to K = I (12)
The current-value Hamiltonian corresponding to [START_REF]Energy Technology RD&D Budgets: Overview[END_REF] is:

H = (1 -τ ) - (1 -τ ) 2 2 - 1 2 x + c x0 -c x1 K x - 1 2 I + c I K I - γ 2 (1 -x -τ ) 2 + λI (13)
where λ denotes the shadow value of the stock of knowledge. For an interior solution, necessary conditions are:

-τ + γ(1 -x -τ ) = 0 (14) -x -c x0 + c x1 K + γ(1 -x -τ ) = 0 (15)
I + c I K = λ (16)
The adjoint equation is:

λ = rλ -c x1 x + c I I (17)
Last, the transversality condition is lim t→+∞ e -rt λ(t)K(t) = 0.

Conditions ( 14) and ( 15) yield:

τ (K) = γ(1 + c x0 ) 1 + 2γ - γc x1 1 + 2γ K (18) x (K) = γ -(1 + γ)c x0 1 + 2γ + (1 + γ)c x1 1 + 2γ K (19)
On the other hand, by combining ( 16) and ( 17), I get the following equation:

İ = rI + rc I K -c x1 x (20)
Then, substituting [START_REF] Popp | ENTICE-BR: The effects of backstop technology R&D on climate policy models[END_REF] into [START_REF] Popp | Energy, the environment, and technological change[END_REF] and rearranging the terms, I get:

İ = rI + rc I - (1 + γ)c 2 x1 1 + 2γ K -c x1 γ -(1 + γ)c x0 1 + 2γ , K = I (21)
Under conditions γ - 21) has a positive saddle-point stable steady state. 8 The solution is:

(1 + γ)c x0 > 0 and rc I -γ 1+2γ c 2 x1 c x0 > 0, equation (
I (K) = θ (K -K ∞ ) (22)
where

K ∞ = c x1 (γ -(1 + γ)c x0 ) rc I (1 + 2γ) -(1 + γ)c 2 x1 (23)
denotes the steady state of the stock of knowledge and:

θ = r 2 - r 2 2 + rc I - 1 + γ 1 + 2γ c 2 x1 ( 24
)
relates to the speed of convergence to the steady state.

The characterization of the social optimum goes as follows. In order to internalize the external cost of pollution, the benevolent social planner must use the tax rate [START_REF]Policy INstruments for the Environment database[END_REF]. This tax yields a positive price for the output, which allows the development of a clean technology as given by [START_REF] Popp | ENTICE-BR: The effects of backstop technology R&D on climate policy models[END_REF]. Thus, the social planner invests according to the feedback rule [START_REF] Soderholm | Empirical challenges in the use of learning curves for assessing the economic prospects of renewable energy technologies[END_REF] to lower the marginal cost related to the clean technology. The development is implemented over time due to the adjustment costs of investment. The fishing out effect on R&D implies that investments are less and less profitable through time, so the stock of knowledge reaches a maximum level given by [START_REF] Tsur | Optimal transition to backstop substitutes for nonrenewable resources[END_REF]. Following the progressive development of the clean technology, the tax rate must decrease while the stock of knowledge accumulates. 8 The economic interpretation of these assumptions is the following. Condition γ -(1 + γ)cx0 > 0 means that the net benefit of the first unit of knowledge is positive (recall that K(0) = 0). Without this first assumption, the R&D is never profitable, and there is no way to reduce the cost differential between the two technologies without lowering the welfare. Condition rcI -γ 1+2γ c 2 x1 c x0 > 0 ensures two properties. First, it means that, at each date, the net benefit of any infinitesimal quantity of knowledge is decreasing in K. It is consistent with the assumption of a fishing out effect: the stock of knowledge becomes less and less profitable as it accumulates. Second, it ensures that cx0 -cx1K ∞ > 0, which means that the cost can not become negative because of the accumulation of knowledge.

3 Analysis in a decentralized economy

Strong commitment: Open-loop Stackelberg equilibrium

In this section, I characterize the Open-Loop Stackelberg Equilibrium of the game, which corresponds to the case of strong commitment. In this equilibrium concept, the government announces a policy at t = 0, i.e. {(τ (t), σ(t)), t ∈ [0, ∞)}, and commits to not to deviate from it. Given such a policy, the current-value Hamiltonian associated with the monopoly's maximization problem [START_REF] Gould | Adjustment of Costs of in the the Theory Investment Firm[END_REF] is:

H = τ (t)x - 1 2 + c x0 -c x1 K x - 1 2 I + c I K I + σ(t)I + λI (25)
where λ denotes the shadow value of the stock of knowledge. Necessary conditions are:

τ (t) -x -c x0 + c x1 K = 0 (26)
I + c I K -σ(t) = λ (27)
The adjoint equation is:

λ = rλ -c x1 x + c I I (28)
Last, the transversality condition is lim t→+∞ e -rt λ(t)K(t) = 0.

The goal is to find a policy {(τ (t), σ(t)), t ∈ [0, ∞)} such that equations ( 26) to (28) induce the efficient path characterized in the previous section. I notice that if the government sets the tax rate at each date as follows:

τ (t) = τ (K(t)) = γ(1 + c x0 ) 1 + 2γ - γc x1 1 + 2γ K(t) (29)
Then, the clean technology level is:

x(t) = γ -(1 + γ)c x0 1 + 2γ + (1 + γ)c x1 1 + 2γ K(t) (30)
Expression [START_REF] Popp | ENTICE-BR: The effects of backstop technology R&D on climate policy models[END_REF] corresponds to socially optimal clean technology level. In other words, x(t) =

x (K(t)).

On the other hand, by combining ( 27) and (28), I get:

İ = rI + rc I K -c x1 x + σ(t) -rσ(t) (31)
Therefore, it is straightforward that if σ(t) = 0, by substituting (30) into (31), I get the same equations than in [START_REF] Popp | Environmental policy and innovation: A decade of research[END_REF] whose resolution yields a stock of knowledge following the efficient path.

I summarize this result in the following Lemma:

Lemma 1. Under a strong commitment level, the government can induce the efficient path in the decentralized economy by using the following tax and subsidy rates:

τ (t) = τ (K(t)), σ (t) = 0 (32)
Wirl (2013) finds a similar outcome (see Proposition 2, p. 7), keeping in mind that he only considers production capacities.

Weak commitment: Global Stackelberg equilibrium

Now, let us characterize the Global Stackelberg Equilibrium of the game, which corresponds to the case of weak commitment. In this equilibrium concept, the government announces a policy rule at t = 0. Such a policy is defined as a tax rate and a subsidy rate that both depend linearly on the stock of knowledge K, i.e. {(τ

(K) = τ 0 + τ 1 K, σ(K) = σ 0 + σ 1 K), K ∈ [0, ∞)}.
Given such a policy, the current-value Hamitonian associated with monopoly's maximization problem [START_REF] Gould | Adjustment of Costs of in the the Theory Investment Firm[END_REF] is:

H = τ (K)x - 1 2 + c x0 -c x1 K x - 1 2 I + c I K I + σ(K)I + λI (33)
where λ denotes the shadow value of the stock of knowledge. Necessary conditions are:

τ (K) -x -c x0 + c x1 K = 0 (34)
I + c I K -σ(K) = λ (35)
The adjoint equation is:

λ = rλ -c x1 x + c I I -τ (K)x -σ (K)I (36)
Last, the transversality condition is lim t→+∞ e -rt λ(t)K(t) = 0.

The goal is to find a policy rule {(τ

(K) = τ 0 + τ 1 K, σ(K) = σ 0 + σ 1 K), K ∈ [0, ∞)} such
that equations (34) to (36) characterize the social optimum. Again, I notice that if:

τ (K) = τ (K) = γ(1 + c x0 ) 1 + 2γ - γc x1 1 + 2γ K (37)
Then, equation (34) yields:

x(K) = γ -(1 + γ)c x0 1 + 2γ + (1 + γ)c x1 1 + 2γ K (38)
This corresponds to the expression [START_REF] Popp | ENTICE-BR: The effects of backstop technology R&D on climate policy models[END_REF] of the clean technology level under the social optimum. In other words, x(K) = x (K). On the other hand, by combining ( 35) and (36), I get:

İ = rI + rc I K -c x1 x -rσ(K) -τ (K)x (39)
I notice that if rσ(K) = -τ (K)x, equation (39) simplifies to:

İ = rI + rc I K -c x1 x (40)
It is straightforward that substituting (38) into (40) yields the same equations ( 21) than under the social optimum, which ensures that the stock of knowledge follows the efficient path.

I summarize this result in the following Proposition:

Proposition 1. Under a weak commitment level, the government can induce the efficient path by setting both the following tax rule

τ (K) = τ (K) = γ(1 + c x0 ) 1 + 2γ - γc x1 1 + 2γ K (41)
and the following subsidy rule:

σ (K) = - τ (K)x (K) r = γc x1 (γ -(1 + γ)c x0 ) r(1 + 2γ) 2 + γc 2 x1 (1 + γ) r(1 + 2γ) 2 K (42)
This result states that if the government can not commit to announcements about future tax rates (strong commitment level) but can rather commit to announcements about a policy rule (weak commitment level), then it can induce the efficient path.

By comparison with Lemma 1, Proposition 1 highlights that, under weak commitment, a subsidy to investment in knowledge is required in addition to the tax. The reason is that the government lack of commitment implies that the firm takes into account that any new unit of knowledge reduces its revenue by decreasing the tax rate, as shown by the term -τ (K)x in equation (39). Therefore, a subsidy is needed to make up for this future marginal loss in revenue. [START_REF] Wirl | Taxes versus permits as incentive for the intertemporal supply of a clean technology by a monopoly[END_REF] considers the lowest possible level of commitment and finds the government can not induce the efficient path. By contrast, it is possible to induce the efficient path in this setting.

Returning to the case of the energy sector, the general message of Proposition 1 is as follows.

Without strong commitment, the government may commit to follow a policy rule and, in order to make up for its lack of commitment, it should subsidy R&D activities in favor of the development of non-polluting sources of energies.

The case of a duopoly

An issue emerging in a game with a stock of knowledge is the issue of spillover effects. In this section, I extend the previous model by considering two identical firms, i = 1, 2, which own two clean technologies. Each firm i owns its private stock of knowledge K i and I denote by I i firm i's individual investment level. I assume that there are spillover effects in the evolution of both stocks. That is, the stocks of knowledge evolve as:

Ki = I i + βI -i , 0 ≤ β ≤ 1, i = 1, 2 (43)
This is a standard approach to knowledge spillovers in a dynamic context, see e.g [START_REF] Cellini | Dynamic R&D with spillovers: Competition vs cooperation[END_REF]. Coefficient β denotes the technological spillover that each firm derives from the R&D activity of the opponent. This applies to a situation where there exist different clean technologies, each of them protected by patents, and when the improvement of each technology benefits to the others. For instance, there exist several types of solar photovoltaic technologies (silicon solar cells, thin film solar cells...).

The clean technology level of firm i, denoted by x i , is associated with a cost depending on the stock of knowledge as follows:

C i x (x i , K i ) = 1 2 x 2 i + (c x0 -c x1 K i )x i , c x0 > 0, c x1 > 0, c x0 -c x1 K i > 0, i = 1, 2 (44)
By investing in the stock of knowledge, firms face adjustment costs of investment:

C i I (I i , K i ) = 1 2 I i + c I K i I i , c I > 0, i = 1, 2 (45)
The instantaneous profit of firm i is equal to the revenue plus the subsidy to investment, net of the production cost and the investment adjustment cost:

π i = px i - 1 2 x i + c x0 -c x1 K i x i - 1 2 I i + c I K i I i + σI i , i = 1, 2 (46)
where p = τ . Thus, its objective is:

max x i ,I i ∞ 0 π i e -rt dt, subject to K1 = I 1 + βI 2 , K2 = I 2 + βI 1 , i = 1, 2 (47)
The benevolent government accounts for the firms' profits, the consumers' surplus, the tax revenue and the external cost of pollution:

S d = (1-τ )- (1 -τ ) 2 2 - i=1,2 1 2 x i + c x0 -c x1 K i x i - i=1,2 1 2 I i + c I K i I i - γ 2 (1-x 1 -x 2 -τ ) 2 (48)
The government chooses tax rate τ and subsidy rate σ to maximize the aggregate welfare:

max τ,σ ∞ 0 S d e -rt dt, subject to K1 = I 1 + βI 2 , K2 = I 2 + βI 1 (49)
As in the case of a monopoly, I consider a game with hierarchical play between the government (first-mover) and the duopoly (second-mover). In the remainder of this section, I characterize the social optimum, the OLSE and the GSE.

The social optimum

A benevolent social planner chooses tax rate τ , clean technology levels (x 1 , x 2 ), investment levels (I 1 , I 2 ) in order to maximize the aggregate welfare subject to the evolution of the stocks of knowledge:

max τ,x 1 ,x 2 ,I 1 ,I 2 ∞ 0 S d e -rt dt, subject to K1 = I 1 + βI 2 , K2 = I 2 + βI 1 (50)
The current-value Hamiltonian corresponding to (50) is:

H = (1 -τ ) - (1 -τ ) 2 2 - i=1,2 1 2 x i + c x0 -c x1 K i x i - i=1,2j 1 2 I i + c I K i I i -γ 2 (1 -x 1 -x 2 -τ ) 2 + λ 1 (I 1 + βI 2 ) + λ 2 (I 2 + βI 1 ) (51)
where λ 1 and λ 2 respectively denote the shadow values of the stocks of knowledge K 1 and K 2 . For an interior solution, necessary conditions are:

-τ + γ(1 -x 1 -x 2 -τ ) = 0 (52)
-x i -c x0 + c x1 K i + γ(1 -x 1 -x 2 -τ ) = 0, i = 1, 2 (53)
I i + c I K i = λ i + βλ -i , i = 1, 2 (54)
Adjoint equations are:

λi = rλ i -c x1 x i + c I I i , i = 1, 2 (55)
Last, the transversality conditions are lim

t→+∞ e -rt λ i (t)K i (t) = 0 (i = 1, 2).
Equations ( 52) and (53) yield:

τ (K 1 , K 2 ) = γ(1 + 2c x0 -c x1 (K 1 + K 2 )) 1 + 3γ (56)
x i (K 1 , K 2 ) = γ -(1 + γ)c x0 + c x1 ((1 + 2γ)K i -γK -i ) 1 + 3γ , i = 1, 2 (57)
On the other hand, by combining (54) to (55), I get:

İi = rI i + rc I K i -c x1 x i -βc x1 x -i , i = 1, 2 (58)
Since firms are identical, I only consider symmetric equilibria where all firms behave identically. It implies that K 1 = K 2 = K/2, I 1 = I 2 = I/2 and x 1 = x 2 = x/2 where K, I and x denote the total stock and control variables. The expression of the tax rate (56) and the total clean technology level (57) become:

τ (K) = γ(1 + 2c x0 ) 1 + 3γ - γc x1 1 + 3γ K (59) x (K) = 2γ -2(1 + γ)c x0 1 + 3γ + (1 + γ)c x1 1 + 3γ K (60)
Futhermore, equation (58) can be written as:

İ = rI + rc I K -c x1 (1 + β)x (61)
Finally, substituting (60) into (61), I get:

İ = rI + rc I -(1 + β) (1 + γ)c 2 x1 1 + 3γ K -c x1 (1 + β) 2γ -2(1 + γ)c x0 1 + 3γ (62)
I (K) = 1 1 + β θ (K -K ∞ ) (63)
where

K ∞ = 2(1 + β)c x1 (γ -(1 + γ)c x0 ) rc I (1 + 3γ) -(1 + β)(1 + γ)c 2 x1 (64)
denotes the steady state stock of knowledge and:

θ = r 2 - r 2 2 + (1 + β) rc I - (1 + β)(1 + γ) 1 + 3γ c 2 x1 ( 65
)
relates to the speed of convergence to the steady state.

Strong commitment: Open-loop Stackelberg equilibrium

In this section, I characterize the Open-Loop Stackelberg Equilibrium of the game, which corresponds to a strong commitment level from the government. The government announces a policy at t = 0, i.e. {(τ (t), σ(t)), t ∈ [0, ∞)}. Given such a policy, and given the strategy of the other firm (q -i , I -i ) the current-value Hamiltonian associated with the firm maximization problem (47) is:

H i = τ (t)x i - 1 2 x i + c x0 -c x1 K i x i - 1 2 I i + c I K i I i +σ(t)I i +λ i (I i +βI -i )+µ i (I -i +βI i ), i = 1, 2 (66)
where λ i denotes the shadow value of the stock of knowledge K i and µ i the shadow value of the stock of knowledge K -i . The necessary conditions are:

τ (t) -x i -c x0 + c x1 K i = 0, i = 1, 2 (67)
I i + c I K i -σ(t) = λ i + βµ i , i = 1, 2 (68)
The adjoint equations are:

λi = rλ i -c x1 x i + c I I i , i = 1, 2 (69)
μi = rµ i , i = 1, 2 (70)
Last, the transversality conditions are lim t→+∞ e -rt λ i K i (t) = 0 and lim t→+∞ e -rt µ i K -i (t) = 0.

The goal is to find a policy {(τ (t), σ), t ∈ [0, ∞)} such that equations (67) to (70) correspond to the social optimum. First of all, condition (70) has no use in this case since the transversality

H i = τ (K 1 , K 2 )x i -1 2 x i + c x0 -c x1 K i x i -1 2 I i + c I K i I i + σ(K 1 , K 2 )I i +λ i (I i + βI -i ) + µ i (I -i + βI i ), i = 1, 2 (77)
where λ i denotes the shadow value of firm i's stock of knowledge and µ i the shadow value of firm -i's stock of knowledge. The necessary conditions are:

τ (K 1 , K 2 ) -x i -c x0 + c x1 K i = 0, i = 1, 2 (78)
I i + c I K i -σ(K 1 , K 2 ) = λ i + βµ i , i = 1, 2 (79)
The adjoint equations are:

λi = rλ i -c x1 x i + c I I i -τ i 1 x i -σ i 1 I i , i = 1, 2 (80)
μi = rµ i -τ -i 1 x i -σ -i 1 I i , i = 1, 2 (81)
The transversality conditions are lim t→+∞ e -rt λ i (t)K i (t) = 0 and lim t→+∞ e -rt µ i (t)K -i (t) = 0.

The goal is to find a policy rule {(τ

(K 1 , K 2 ) = τ 0 + τ 1 1 K 1 + τ 2 1 K 2 , σ(K 1 , K 2 ) = σ 0 + σ 1 1 K 1 + σ 2 1 K 2 ), (K 1 , K 2 ∈ [0, ∞))
} such that equations (34) to (36) correspond to the social optimum. Again, I notice that if:

τ (K 1 , K 2 ) = τ (K 1 , K 2 ) = γ(1 + 2c x0 -c x1 (K 1 + K 2 )) 1 + 3γ (82)
Then the clean technology levels are:

x i (K i , K -i ) = γ -(1 + γ)c x0 + c x1 ((1 + 2γ)K i -γK -i ) 1 + 3γ , i = 1, 2 (83)
This levels correspond to the expression (57) of the clean technology levels in the social optimum:

x i (K i , K -i ) = x i (K i , K -i ).
On the other hand, by combining (79), ( 80) and (81), I get:

İi = rI i + rc I K i -(c x1 + τ i 1 + βτ -i 1 )x i + (σ -i + βσ i -βc I )I -i -rσ(K 1 , K 2 ), i = 1, 2 (84)
The firms being identical, equation (84) becomes:

İ = rI + rc I K -(c x1 + (1 + β)τ 1 )x + (-βc I + (1 + β)σ 1 )I -2r(σ 0 + σ 1 K), K = (1 + β)I (85)
where

τ 1 = τ i 1 = τ -i 1 and σ 1 = σ i 1 = σ -i 1 .
I notice that if subsidy rate σ(K) = σ 0 + σ 1 K satisfies:

(1 + β)σ 1 I(K) -2r(σ 0 + σ 1 K) = β (c I I(K) -c x1 x(K)) + (1 + β)τ 1 x(K) (86)
Then, equation (85) simplifies to the same equation under the social optimum (61), which solution yields a stock of knowledge following the efficient path. Finally, I obtain the following Proposition:

Proposition 2. Under a weak commitment level, the government can induce the efficient path by setting both the following tax rule

τ (K) = τ (K) = γ(1 + 2c x0 -c x1 K) 1 + 3γ (87)
and the following subsidy rule:

       σ (K) = σ (K) -(1+β) 2 τ 1 r x (K) + φ(K) with φ(K) = -θ /2 r-θ /2 β 2(1+β) c I + θ 1+β + (1+β) 2 τ 1 r (1+γ)c x1 1+3γ (K -K ∞ ) (88)
The proof of Proposition 2 is provided in Appendix B.

When the government commits to a policy rule, the subsidy rate (88) is made of three components. The first component, i.e. σ (K), corrects the free-rider problem as under strong commitment, see (76). The second component, i.e. -(1+β) 2 τ 1 r x (K), corrects the fact that firms take into account that any new unit of knowledge reduces their future revenues because of the evolution of the tax rate as in the case of a monopoly, see (42).

The third effect, related to φ(K), is novel. It is used to correct an indirect effect on the subsidy. As for the tax rate, it results from the fact that both firms anticipate that any new unit of knowledge affects the future subsidy rate. This induces the firms to invest faster in the short run. The component φ(K) of the subsidy rate corrects this indirect effect in the short run. One may check that this component equals zero in the long run, i.e. φ(K ∞ ) = 0.

A consequence of this novel effect is that, in some situations, one may have σ (K) < σ (K) in the short run. Surprisingly, it can even be the case that σ (K) be initially negative to slow down the inefficiently high pace of investments. This is what I summarize in the following corollary:

Corollary 1. Under a weak commitment level, the government may initially tax investments in knowledge made by the firms in order to slow down an otherwise too fast accumulation of knowledge.

Figure 1 illustrates this last result. In the case depicted here, the pace of knowledge investment is too high. Then, the government subsidies less under weak commitment (GSE) than under strong commitment (OLSE) in the short run. Finally, in the very short run, it actually taxes investment in the GSE case.

Concluding remarks

I show that if a government can not commit to a policy in the long-run, one can find a policy rule that induces a monopoly to efficiently accumulate the stock of knowledge required for the development of a new clean technology. This policy is made of a tax and subsidy mix. The tax is used to internalize the external cost due to the pollution related to the dirty technology. The subsidy corrects the distorsion that results from the firm anticipation of how its decisions affect the magnitude of the tax rate.

Then, extending the model to the case of a duopoly with knowledge spillovers, I characterize the policy rule that induces efficiency: the subsidy rate may be initially lower under weak commitment than under strong commitment. The reason is that each firm, taking into account that the other firm investment affects the subsidy rate, may be induced to invest at a too high pace initially. A surprising result is that, under weak commitment, the subsidy rate may even be negative to disinduce firms to invest too much in the short run. A general take away of this paper is that abstracting from the issue of the government commitment level yields suboptimal policies. In particular, this may lead one to either underestimate or overestimate subsidy levels allocated to R&D activities.

Of course, this contribution calls for further work. In particular, I do not analyse the case of a closed-loop information structure in the duopoly setting. Such an analysis is very likely to rely on numerical simulations only, due to the presence of two stock variables. I also do not consider the case of a stock of pollution: again, such an extension would require extensive simulation work due to the presence of two stock variables. Such extensions are left for future research work.

A Proof of Lemma 2

I have shown that σ(t) must be solution to (75) on the efficient path. Then, substituting (63) and ( 60) into (75), I get:

σ(t) -rσ(t) = β 2 c I θ 1 + β (K(t) -K ∞ ) -c x1 2γ -(1 + γ)(2c x0 -c x1 K(t)) 1 + 3γ (89)
Now, along the efficient path the stock of knowledge evolves as K(t) = K ∞ (1 -e θ t ) so that equation (90) can be written as:

σ(t) -rσ(t) = β 2 -c I θ 1 + β K ∞ e θ t -c x1 2γ -(1 + γ)(2c x0 -c x1 K ∞ (1 -e θ t )) 1 + 3γ (90)
Rearranging the terms:

σ(t) -rσ(t) = β 2 -c I θ 1 + β + (1 + γ)c 2 x1 1 + 3γ K ∞ e θ t - βc x1 2 2γ -(1 + γ)(2c x0 -c x1 K ∞ ) 1 + 3γ (91)
Standard methods lead to the following solution:

σ(t) = β 2 -c I θ 1 + β + (1 + γ)c 2 x1 1 + 3γ K ∞ θ -r e θ t + βc x1 2r 2γ -(1 + γ)(2c x0 -c x1 K ∞ ) 1 + 3γ (92)
The expression of σ(t) given in (92) can be simplified. Indeed, using (65), I note that:

(1 + γ)c 2 x1 1 + 3γ = θ (r -θ ) (1 + β) 2 + rc I 1 + β (93)
Moreover, using (64), I notice that:

2c x1 (γ -(1 + γ)c x0 ) 1 + 3γ = K ∞ 1 + β rc I - (1 + β)(1 + γ)c 2 x1 1 + 3γ = -K ∞ θ (r -θ ) (1 + β) 2 (94)
Finally, substituting (93) and ( 94) into (92), I get:

σ(t) = - βK ∞ 2(1 + β) c I + θ 1 + β e θ t + βK ∞ c I 2(1 + β) (95)
Or, in feedback form:

σ(t) = - β 2(1 + β) θ 1 + β K ∞ + c I + θ 1 + β K(t) (96)

B Proof of Proposition 2

Plugging the expressions of I (K) and x (K) given in ( 63) and ( 60) into (86), I get a polynomial expression of degree one in K. Then, I collect terms involving K and terms not involving K and get the following system of equations:

       (θ /2 -r) σ 1 = β 2 c I θ 1+β - (
+ (1+γ)c x1 1+3γ K ∞ (99)
Rearranging the terms:

       σ 1 = β 2(1+β) c I + θ 1+β -(1+β)

CEE-M Working Papers

x

  + c x0 -c x1 K. The term x accounts for diseconomies of scale in the production process. The supply of a large amount of electricity requires the operation of a large number of plants, including the less efficient ones, which increases the cost. Consistently with Tsur and Zemel (2003), the term c x0 -c x1 K models the progressive improvement of the clean technology. The stock of knowledge increases with the investment I undertaken by the monopoly:

2 x1c

 2 Under conditions γ -(1+γ)c x0 > 0 and rc I -(1+β) γ 1+3γ c x0 > 0, equation (62) has a positive saddle-point stable steady state. The solution is:

Figure 1 :

 1 Figure 1: Comparison of the subsidy rates under strong commitment level (OLSE) and under weak commitment level (GSE) in the case of a duopoly with c x0 = 0.1, c x1 = 0.05, γ = 0.2, c I = 1, r = 0.05 and β = 0.8.

-θ K ∞ σ 1 2 --rσ 1 K 2 ( 1 r-θ / 2 (1+β)τ 1 2 (

 212122 rσ 0 = -∞ -rσ 0 = -β ) c I + θ 1+β -1+γ)c x1 1+3γ -rσ 1 K ∞ -rσ 0 = -β 2(1+β) rc I K ∞ + (1+β)τ 1 2 2γ-2(1+γ)c x0 1+3γ

σ 0 =

 0 -σ 1 K ∞ + β 2(1+β) c I K ∞ -the terms, I obtain expression (88).

  The first equilibrium concept is the Open-Loop Stackelberg Equilibrium(Dockner et al., 2000, p. 113). It refers to a situation where the government announces at t = 0 a time-path strategy for the tax and the subsidy, i.e. {(τ (t), σ(t)), t ∈ [0, ∞)}. It requires that the government can commit to the future path of tax and subsidy rates. It models a strong commitment level from the government.The second equilibrium concept is the Global Stackelberg equilibrium(Basar and Olsder, 1999, p. 373).6 It refers to a situation where the government announces a policy rule. Specifically, it is standard to consider a (stationary) Linear Markov Rule[START_REF] Dockner | Differential games in economics and management science[END_REF][START_REF] Benchekroun | efficiency-inducing taxation for polluting oligopolists[END_REF][START_REF] Karp | On efficiency-inducing taxation for a non-renewable resource monopolist[END_REF], i.e. {(τ

To keep notations as simple as possible, note that I generally omit the time variable "t".

The learning rate denotes the pourcentage change in cost reductions resulting from a doubling of either the stock of knowledge (for R&D) or the cumulative capacities (for LBD).

Thomas Edisson invented the light bulb on his own. Today, improving a technology as solar panel requires research teams.

See alsoDockner et al. (2000, p. 142) where this equilibrium concept is refered to as a "nondegenrate Markovian Stackelberg equilibrium".

At this stage, the reader may object that the subsidy rate σ is absent from[START_REF] Grimaud | Climate change mitigation options and directed technical change: A decentralized equilibrium analysis[END_REF]. However, keep in mind that my analysis will further focus on a game with hierarchical play. The government being the first-mover, it takes into account the response of the monopoly to its policy. In such a context, subsidy rate σ does step in the government objective function.

This work was publicly funded through ANR (the French National Research Agency) under the "Investissements d'avenir" programme with the reference ANR-16-IDEX-0006, and by the ANR project GREEN-Econ (Grant ANR-16-CE03-0005).

condition implies that µ i = 0 for all t. Furthermore, I notice that if:

Then the clean technology levels are:

This corresponds to the expression (57) of the clean technology levels under the social optimum:

. Now, by combining (68) and (69), I get:

The firms being identical, equation (73) becomes:

I notice that if σ(t) is solution to:

Then, equation (74) simplifies to the same equation than under the social optimum (61), which solution yields a stock of knowledge following the efficient path. Finally, I obtain the following Lemma:

Lemma 2. Under a strong commitment level, the government can induce the efficient path in the decentralized economy by using both a tax rate and a subsidy rate satisfying:

The proof of Lemma 2 is provided in Appendix A. By comparison with the case of a monopoly, the government policy includes a subsidy even under a strong commitment level. This latter is needed to internalize the free-rider problem that emerges from the existence of spillovers. One may check that in the absence of spillover effects, i.e. if β = 0, no subsidy is needed.

Weak commitment: Global Stackelberg Equilibrium

Now, let us characterize the Global Stackelberg Equilibrium of the game, which corresponds to the case of weak commitment. The government announces a policy rule at t = 0: such a policy is defined as a tax rate and a subsidy rate which both depend linearly on the two stocks of knowledge

Given such a policy, and given the strategy of the other firm (x -i , I -i ) the current-value Hamiltonian associated with the firm maximization problem ( 9) is:


Bernard Makhraz

Hubert Piquet •

Xavier Roboam

Jerome Mavier

Matrix approach based on Quadripole for quality analysis in Aircraft Electrical Power Distribution System

Keywords: Matrix approach, More-Electric Aircraft, Power Distribution, Quality, HVDC Standards

HVDC bus based electrical network architectures have shown to be relevant for More-Electric Aircraft. Today, adequate standards of power quality and stability must be set and even optimized: subsequent design constraints have to allow the system integrator to guarantee the safe operation of the network, without excessive impact on the design of the equipment provided by the suppliers, especially regarding weight penalty. In this context, this paper proposes a specific methodology to analyze the power quality constraints in an HVDC network. A matrix approach, based on Quadripole analysis of the power distribution is considered. It aims at achieving a direct sensitivity analysis of any quantity (regarding its harmonic content) of any electrical network with respect to selected inputs and based on the network topology and parameters. One major advantage of this methodology is that the modeling approach is analytically developed, based on symbolic calculation. All the devices (sources, cables, filters, loads, etc) of the system are represented as four-terminal (quadripole) devices, connected together according to the network topology. The calculations are performed at very high speed in the frequency domain. An application of this approach on a relevant aircraft electrical power distribution is shown to highlight its benefit on the path of optimizing new HVDC quality standards.

Introduction

The More-Electric Aircraft embeds more electric power systems aiming at minimizing both sizes and weights of components (especially in the filter device) by fulfilling quality and stability standards [START_REF] Liu | Power Quality and Stability Issues in More-electric Aircraft Electrical Power Systems[END_REF] [START_REF] Girinon | Analytical input filter design in DC distributed power systems approach taking stability and quality criteria into account[END_REF], as well as energy management optimization by using the same power inverter to supply different loads during different flight phases [START_REF] Giraud | Optimal design of the Integrated Modular Power Electronics Cabinet Aerospace Science and Technology[END_REF] [START_REF] Etayo | Methods and tools for the optimization of modular electrical power distribution cabinets in aeronautical applications[END_REF]. New distribution architecture, for instance using HVDC buses to save weight [START_REF] Roboam | More electricity in the air: Towards optimized electrical networks embedded in "more electrical aircraft[END_REF] -though local HVDC bus bars are used in the current state of art [START_REF] Sinnett | Saving Fuel and Enhancing Operational Efficiencies[END_REF]. Using a +/-270V DC network is a promising way to save weight as well [START_REF] Brombach | Optimizing the Weight of an Aircraft Power Supply System through a +/-270 VDC Main Voltage[END_REF]. These implications entail new challenges on the quality and the stability of the DC network [START_REF] Girinon | Analytical modeling of the input admittance of an electric drive for stability analysis purposes[END_REF] [START_REF] Charrada | Equipments characterization methods for stability analysis of DC networks IEEE International Symposium on Industrial Electronics[END_REF] due to the use of non-linear loads, such as power converters, [START_REF]IEEE Recommended Practice and Requirements for Harmonic Control in Electric Power Systems[END_REF] [START_REF] Blooming | Application of IEEE STD 519-1992 Harmonic Limits[END_REF]. Interoperability between the sources and the loads falls under the responsibility of the system integrator through the HVDC standard, which defines the power quality requirements that need to be satisfied under different operating conditions. Today, standards for new distribution networks (HVDC) have to be studied and even optimized in order to avoid weight surplus.

In this context, this paper propose an analysis methodology based on a quadripole (matrix) approach. One major advantage of this method is that the modeling approach is analytically developed, based on symbolic calculation which allows to directly analyze network couplings in the frequency domain, leading to fast and direct analysis.

In this paper, section 1 introduces the four-terminal model used for all the building blocks of the HVDC network and B. MAKHRAZ • J. MAVIER Airbus Operation SAS Toulouse, France E-mail : bernard.makhraz@airbus.com, jerome.mavier@airbus.com H. PIQUET • X. ROBOAM • B. MAKHRAZ Laboratoire Plasma et Conversion d'Energie (LAPLACE) Toulouse, France E-mail : piquet@laplace.univ-tlse.fr, roboam@laplace.univ-tlse.fr its acausal formulation, with generalization at system level. The topology of the network is taken into account in section 2, using matrix formulation. Calculations presented in section 3 allow drawing the correlation and sensitivity analysis between any pair of devices ("disruptive" and "victim" devices), to evaluate the impact at the victim level of the environment pollution, especially fulfillment of standard requirements in the design of the disruptive element. An illustrative example at system level is presented in section 4. Finally, an approach integrating standard's requirements currently under study is proposed in section 5.

Quadripole based model definition 1.Four-terminal model definition

The modelling approach for quality analysis is based on representing of each component of the electrical distribution network with a four-terminal circuit (also sometimes called "two-port network"), which is represented by a "quadripole" matrix. This matrix, involving two linear equations, links currents and voltages on both sides of the concerned element as shown in figure 1. As the matrix representing each quadripole can have more than one format (based on the different input/output combinations of currents and voltages), the first matrix introduced in the equation ( 1) is the "Transfer Matrix" which connects the voltage and the current on the right side with those on the left side.

v o i o ] = [ T 11 T 12 T 21 T 22 ] * [ v i -i i ] (1)
One should notice that the terms of the quadripole matrix might take into account nonlinear phenomena in the frequency domain.

Quadripole acausal Matrix (QM)

Thereafter, the linear system based on the transfer matrix is transformed into an acausal set of equations showing both variables and a zero vectors, as shown in equation (2 & 3).

[ 1 0 -T 11 T 12 0 1 -T 21 T 22 ] * [ v o i o v i i i ] = [ 0 0 ] (2)
[QM] * [

v o i o v i i i ] = [ 0 0 ] (3)
For instance, equation (4) represents the Quadripole acausal Matrix (QM) of a LC filter configuration with the addition of a shunt Rd-Cd as a damping network. This configuration is represented in figure 2.

𝑅 𝑑 + 1 𝐶 𝑑 𝑝 ) + 1] * [ v o i o v i i i ] = [ 0 0 ] (4)
1.3. System Quadripole Matrix (SQM) With the Quadripole acausal Matrix of each component, the "System Quadripole Matrix" is formed by gathering all of these quadripole matrices into one global matrix representing all the building blocks of the distribution network under study. In the System Quadripole Matrix, only the block diagonal terms are different from zero. For instance, the System Quadripole Matrix of an electrical distribution network that includes 'n' elements looks like the matrix shown in equation ( 5) which gathers all voltage current variables of the system.

[ QM 1 0 ⋯ 0 0 QM 2 0 ⋯ ⋯ 0 0 0 0 0 0 QM n ] * [ v o 1 i o 1 v i 1 i i 1 v o 2 i o 2 v i 2 i i 2 ⋮ v o n i o n v i n i i n ] = [ 0 ⋮ 0 ] (5)
Note that the System Quadripole Matrix only integrates quadripole elements of the network without considering the topology which results from their connections; this point is considered in the following section.

Structural Relationships Matrix (SRM)

This second step takes into account the electrical distribution network architecture. Thereafter, structural relations are the equations that identify the connections between the different quadripoles. These relations are established using both Kirchhoff's current (KCL) and voltage (KVL) laws. Therefore, the Structural Relations Matrix is mapped for the same variable's vector used to define the System Quadripole Matrix.

The number of Structural Relationships

The number of structural relations depends on two factors: -The first one is 'p', the number of nodes in the network; -The second one is the number of quadripoles connected at each node. For instance, for a node connecting 'k' quadripoles, the number of structural relations is one for the currents (KCL) and 'k-1' for the voltages (KVLs). Therefore, the total structural relation matrix is obtained by summing all structural relations for each node of the system, as shown in equation [START_REF] Sinnett | Saving Fuel and Enhancing Operational Efficiencies[END_REF].

𝑇𝑜𝑡𝑎𝑙 𝑠𝑡𝑟𝑢𝑐𝑡𝑢𝑟𝑎𝑙 𝑅𝑒𝑙𝑎𝑡𝑖𝑜𝑛𝑠 = ∑ 𝑘 𝑛𝑜𝑑𝑒 𝑖 𝑝 𝑖=0 (6)
An example of Electric Power Distribution, as shown in figure 3, is introduced in this section as a simple application. Two nodes are identified (#1 and #2) in a three levels electric power distribution. The first node connects 4 quadripoles; the second one connects 3 quadripoles. In this case, the number of structural relationships is 4 for the first node and 3 for the second node, resulting in a total of 7 structural relationships, as illustrated by the table 1.

Node #1 Node #2 v 0 -v 1 ′ = 0 v 3 -v 4 ′ = 0 v 0 -v 2 ′ = 0 v 3 -v 5 ′ = 0 v 0 -v 3 ′ = 0 i 3 + i 4 ′ + i 5 ′ = 0 i 0 + i 1 ′ + i 2 ′ + i 3 ′ = 0

Electric Power System Matrix (EPSM)

In order to get a matrix describing the electrical power system under analysis, both the System Quadripole Matrix and the Structural Relations Matrix, which share the same voltage-current vector, are vertically merged in equation [START_REF] Brombach | Optimizing the Weight of an Aircraft Power Supply System through a +/-270 VDC Main Voltage[END_REF].

[EPSM] * [ v o 1 i o 1 v i 1 i i 1 v o 2 i o 2 v i 2 i i 2 ⋮ v o n i o n v i n i i n ] = [ SRM SQM ] * [ v o 1 i o 1 v i 1 i i 1 v o 2 i o 2 v i 2 i i 2 ⋮ v o n i o n v i n i i n ] = [ 0 ⋮ 0 ] ( 7
)
To go further, the aim is now to identify the "causality for the analysis", by establishing the input/output matrix related to the input variables chosen for analysis with all other variables, which can be calculated as the system's outputs.

The rank of the Electrical Power System Matrix determines the number of output variables, which are arbitrary selected. Accordingly, the input variables are the remaining system's variables. Thereafter, the matrix connecting the inputs with the output variables is identified based on the merged Electrical Power System Matrix (equation ( 7)). Thus, the causal input-output matrix, equation [START_REF] Girinon | Analytical modeling of the input admittance of an electric drive for stability analysis purposes[END_REF], is introduced as the "Transfer Function Matrix" relating output vector with the desired input vector.

Supposing that 'r' is the rank of the Electrical Power System Matrix, in a system containing 'm' variables (Var), which can be either voltage or current, the "Transfer Function Matrix" is derived in equation ( 8).

(

) 8
For instance, the Electrical Power System Matrix's rank of the 24 variables example introduced in figure 3 is 19. Therefore, 19 variables are to be considered as output and 5 variables as input.

Note that all calculations developed in the sections 2 and 3 can be derived analytically by means of symbolic calculation solvers; for that purpose, we have selected the Matlab's Symbolic Toolbox. In our case, the Transfer Function Matrix, in equation ( 8), is the output of the 'solve'

[ Var out 1 ⋮ Var out i ⋮ Var out r ] = [TFM] * [ Var in 1 ⋮ Var in m-r ]
function (equations and systems solver from the Symbolic Toolbox) that solves a system of equations: it is applied for solving the Electrical Power System Matrix (equation ( 7)) considering specified output and input vectors.

Application

The application shown in this part is a simplified HVDC electric power distribution for More-Electric Aircraft integrating +/-270 VDC voltage bus with only two loads connected in a star configuration, as shown in figure 4. In this application, each load connected to the bus through an input filter associates an inverter fed Permanent Magnet Synchronous Motor (3-Phase PMSM-1kW).

This topology is simplified in order to clarify the relevance of our modelling approach based on quadripole and frequency domain simulation and its ability to facilitate the power quality analysis. More typically, the issue is to clarifiy how this quadripole based approach allows directly analyzing the influence of power quality standard requirements on the quality of waveforms (harmonic content, rms values, resonances, etc) at the inputs/outputs of sources and loads and consequently on the weights of the connected electrical devices (especially the filters). This application case is built of three quadripoles. The source quadripole includes one capacitor (C 0 = 2mF) connected to two inductors (L 0 = 2.8uH) representing the harness. Equation ( 9) represents its quadripole matrix (QM0).

[ 1 0 -(C 0 L 0 p 2 + 1) L 0 p + L 0 p(C 0 L 0 p 2 + 1) 0 1 -C 0 p C 0 L 0 p 2 + 1 ] * [ v′ 0 i′ 0 v 0 i 0 ] = [ 0 0 ] (9)
On the other hand, both load quadripoles represent two different damped LC input filters which the structure is on Fig. 2. and which values after a weight optimization are:

(L 1 = 190uH, C 1 = 2.3uF, Cd 1 = 0.57uF, Rd 1 = 40Ω, L 2 = 306uH, C 2 = 1.6uF, Cd 2 = 0.7uF, Rd 2 = 42Ω).
QM1 and QM2 are the same matrix as in equation ( 4).

In the process developed in previous sections, the next step is to set up the System Quadripole Matrix, using the three Quadripole Matrices (QM0, QM1 and QM2) as shown in equation ( 10).

(

) 10
The network of our case study is then transformed into a quadripole schematic, as shown in figure 5. For instance, choosing the source voltage and both load currents [v 0 , i 1 , i 2 ] as the analysis input vector implies that the remaining nine variables are considered as outputs, as in the equation (15) which sets the causality of the analysis.

[

QM 2 0 0 0 QM 1 0 0 0 QM 0 ] * [ v 2 i 2 v ′ 2 i ′ 2 v 1 i 1 v ′ 1 i ′ 1 v ′ 0 i ′ 0 v 0 i 0 ] = [ 0 ⋮ 0 ] v 0 ′ -v 1 ′ = 0 v 0 ′ -v′ 2 = 0 i 0 ′ + i 1 ′ + i 2 ′ = 0 (14)
(15)

Among the capabilities of the analysis method one can emphasize: first, the analytical (symbolic) derivation of the transfer functions (Mi,j) in the frequency domain which clarifies the analysis of main couplings (typically resonant effects, etc.) between network elements. Secondly, based on the harmonic content (magnitude and phases) of the input data (here v0,i1,i2), the harmonic content of selected outputs here v2) can be numerically estimated to quantify the coupling analysis as displayed in figure 6. The electric network of figure 4 was also implemented using PLECS ® , a time domain simulator, in order to compare and validate the results. Each load is constituted of a power drive system is here represented by a current source which harmonic content represents the harmonic pollution that would be drawn by a 1kW permanent magnet synchronous machine fed by a voltage source inverter.

The simulation validation through the frequency analysis of the load voltage v2 is represented in figure 6, comparing both PLECS (red spectrum) and Quadripole model (in blue) for a frequency range of [0 ; 100kHz]: the results are completely overlaid, with a negligible difference lower than 0.4% for the main harmonics.

quality analysis based on power quality standards

Civil and military aircraft manufacturers define standards, such as the MIL-STD-704 for military aircraft, in order to maintain the power quality of AC and DC electric power network and to ensure interoperability at system level. These standards set the interface between sources and loads to ensure compatibility between the aircraft electric system, external power, and airborne utilization equipment. Therefore, it give designers the specifications that they need to respect in order to ensure that an airborne system will be compatible with the power supply provided by the aircraft manufacturer.

Standards define the characteristics for different operational conditions such as power up, steady state and reconfiguration in normal, abnormal and emergency condition. However, no appropriate international standards are currently set to evaluate the power quality of the HVDC network, especially for recent concepts integrating +/-270 VDC voltage distribution in aircraft [START_REF] Moir | Aircraft Systems: Mechanical, Electrical, and Avionics Subsystems Integration[END_REF].

In this part, the aeronautical standard's requirements are "directly" integrated as the inputs of our analysis approach: the issue is to assess, thanks to the quadripole based analysis method, the influence of standards on the outputs selected to be analyzed. For this purpose, we first consider that the power source and the load 1 fulfill the quality standards: "it is like if these devices would pollute at the limit level allowed by the quality standards". Thereafter, [v 0 ′ , i 1 ′ ] are then imposed as the input vectors by considering that the output voltage (v 0 ′ ) exactly involves the frequency spectrum of fig 7 while the input current of the 1 st load (i 1 ′ ) has the envelope of the current requirements of Fig. 8. On the other hand, the actual emitted current pollution by the load 2 (which is under analysis) is considered as a third input.

The current i 1 ′ is replaced by the maximum current spectrum authorized by the power quality requirement as shown in figure 8. This current envelope represents the maximum pollution that this load "may emit". In the same way, the voltage v 0 ′ is considered as an input with the maximum voltage distortion spectrum; a typical envelope is represented in figure 7. Thus, the idea is to analyze the constraint on the device N°2, especially its voltage v 2 (shown in fig. 8) when its network environment (source and load 1) fulfills the power quality requirements.

[ By this way, the approach allows "directly studying" the worst-case scenario with the impact of any quality requirements on the electrical device under analysis: in figure 9, the output voltage of load 2 is analyzed regarding the requirements previously considered and the current pollution of the load 1. This information is relevant for the designer of the power drive system of the analyzed load.

0 0 0 0 0 0 -1 0 1 0 0 0 0 0 -1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0 1 0 0 ] 𝑡 * [ v 2 i 2 v ′ 2 i ′ 2 v 1 i 1 v ′ 1 i ′ 1 v ′ 0 i ′ 0 v 0 i 0 ] = [ 0 ⋮ 0 ] [ v 2 v 2 ′ i 2 ′ v 1 v 1 ′ i 1 ′ i 0 v 0 ′ i 0 ′ ] = [ M 1,1 M 1,2 M 1,3 ⋮ ⋱ ⋮ M 9,1 ⋯ M 9,1 ] * [ v 0 i 1 i 2 ]

Conclusion

The main assets of our strategy are the capability to perform symbolic and analytical calculations in the frequency domain based on selected inputs directly showing the main couplings on outputs: the description of the network topology as well as the first steps of the calculations are achieved using symbolic equations which are solved by specific symbolic calculation programs (e.g. Matlab and Maple). According to various scenarios, the mathematical relationships correlating input -output quantities can be analytically built out. In a second step, the actual implementation of each quadripole is taken into account, using the numerical values of the electrical components. Using this approach to merge electrical distribution network with aeronautical standards may simplify the power quality analysis and may be ultimately used to reconsider and even to optimize power quality requirements, since integrating frequency spectrum into simulation software can be a tedious work to do. Using the proposed approach will simplify -eventually -the compliance to the requirements analysis.

Fig. 1 -

 1 Fig. 1 -Schematic representation of a generic electric Quadripole.

[

Fig. 2 -

 2 Fig. 2 -LC filter with Rd-Cd damping branch.

Fig. 3 -

 3 Fig. 3 -Schematic representation of a simplified electric power distribution with three distribution levels.

Fig. 4 -

 4 Fig. 4 -A simplified HVDC electric power distribution for More-Electric Aircraft integrating a +/-270 VDC voltage bus with 2 loads.

Fig. 5 -

 5 Fig. 5 -Schematic representation of a simplified HVDC electric power distribution with one source and two loads.For this n=1 node architecture associating k=3 quadripoles, the Structural Relations Matrix reflects Kirchhoff's laws for currents (1 equation) and voltages (3-1=2 equations): (11) (12) (13)The equations above are gathered in one matrix (14). The aggregation of both System Quadripole Matrix with the Structural Relations Matrix gives the Electrical Power System Matrix with 12 relationships. The rank of this matrix being 9, 12-9=3 variables have to be selected as the analysis inputs and nine as outputs.

Fig. 6 -

 6 Fig. 6 -Validation (by simulation) of the load voltage v 2 : blue harmonics represent the frequency domain calculation; these latter are superposed with harmonics (in red) issued from the PLECS® time simulator. Zooms are proposed at 20kHz, 40kHz and 60kHz

Fig. 7 -Fig. 8 -

 78 Fig. 7 -Frequency spectrum characteristics for load Voltage specification

Fig. 9 -

 9 Fig. 9 -Voltage 𝑣 2 for a frequency range of [0 ; 150kHz] when the source and the load 1 fulfill the quality requirements.

Table 1 -

 1

list of the structural relationships for nodes #1 and #2.


Nellyana Gonzalo Flores

Eric Goncalvès

Regiane Fortes

Julien Rolland

Claude Rebattet

Head Drop of a Spatial Turbopump Inducer

A computational fluid dynamics model for cavitation simulation was investigated and compared with experimental results in the case of a three-blade industrial inducer. The model is based on a homogeneous approach of the multiphase flow coupled with a barotropic state law for the cool water vapor/liquid mixture. The numerical results showed a good prediction of the head drop for three flow rates. The hydrodynamic mechanism of the head drop was investigated through a global and local study of the flow fields. The evolution of power, efficiency, and the blade loading during the head drop were analyzed and correlated with the visualizations of the vapor/liquid structures. The local flow analysis was made mainly by studying the relative helicity and the axial velocity fields. A first analysis of numerical results showed the high influence of the cavitation on the backflow structure.

Introduction

Cavitation occurs frequently in the axial inducer stage of rocket engine turbopumps. It is initiated by a pressure decrease due to a high-speed local liquid velocity and can lead to a fatal failure in pump performance. In order to improve the design method and to evaluate the performance and application limits of inducers working under cavitating conditions, experimental and numerical works have been carried out by some research teams, as for example Refs. ͓1-8͔. In complement to experimental observations, numerical approaches enable flow local analyses and the prediction of global performances.

Concerning three-dimensional ͑3D͒ numerical studies, cavitating flows in turbomachinery are generally modeled by a homogeneous fluid assumption through the one-fluid Reynolds-averaged Navier-Stokes equations ͑RANS͒. Different methods have been proposed to model the mixture and mass transfer between the liquid and vapor. We can cite the interface tracking method proposed by Hirschi et al. ͓9͔, the use of a state law to close the system ͓3,4,10,11͔, the introduction of an additional equation including source terms for vaporization, and the condensation processes applied by Medvitz et al. ͓12͔,Ait Bouziad et al. ͓5,6͔,Mejri et al. ͓2,13,14͔,or Athavale et al. ͓15͔.

Although the numerical modeling of such a phenomenon has received a great deal of attention, it is still a very difficult and challenging task to predict 3D complex cavitating flows with an acceptable accuracy. Physical and numerical calculations should be improved and validated. Methodologies enabling 3D flow analyses and design criteria should be proposed.

In this context, the aim of the present study is to endeavor to improve previous works ͓4,10,11,16͔ by ͑a͒ carrying out complementary observations and obtaining new experimental data in the case of a three-blade inducer, ͑b͒ applying and testing implemented physical and numerical models in the case of this inducer geometry running under different steady cavitating flow conditions, and ͑c͒ proposing an appropriate 3D analysis methodology.

The numerical simulation of cavitating flows in turbomachinery has been developed at the turbomachinery and cavitation team of Laboratoire des Ecoulements Géophysiques et Industriels ͑LEGI, Grenoble, France͒ in collaboration with the French space agency, Centre National d'Etudes Spatiales ͑CNES͒, and the rocket engine division of Snecma ͑Vernon, France͒, the civil and military aircraft manufacturer, a SAFRAN Group company. Studies are based on a homogeneous approach coupled with a barotropic state law to model the cavitation phenomenon ͓17͔. The physical model has been previously integrated by the LEGI laboratory in the computational fluid dynamics ͑CFD͒ code FINETM/TURBO ͓11,18,19͔, developed and commercialized by Numeca International ͑Bruxelles, Belgium͒. Experimental works have been conducted at the Centre d'Essais de Machines Hydrauliques ͑CREMHyG͒ laboratory ͑Grenoble, France͒. This paper presents obtained results concerning the head drop, blade load, power, and hydraulic efficiency evolution in the case of a three-blade inducer. Local analyses of the cavitation influence on secondary flows and backflow are also developed.

Experimental Data

The present study refers to a three-blade industrial inducer, CREM1, illustrated in Fig. 1. It was designed with the goal to suppress the rotating cavitation instabilities ͓16͔. It is worth noting that, for this geometry, the nominal flow rate is much smaller than the zero-incidence angle flow rate.

Water tests were performed at CREMHyG. The TM2 bench is devoted to the characterization of high-speed inducers. Timeaveraged pressure measurements were performed at the shroud with two pressure taps, whose position is indicated in Fig. 2, connected to Druck 910-type pressure piezoresistive transducers with a relative accuracy of 0.2%. The volumetric flow rate is measured by an automatic oil tool ͑AOT͒-type turbine flowmeter, with a relative accuracy of 0.3%. The rotational speed is measured with an HBM T32-type tachymeter with a relative accuracy of 0.1%. From these measurements, the pressure coefficient and the cavitation parameter can be evaluated with an accuracy smaller than 0.4%.

The inducer is placed in a transparent liner housing to enable the observations of the cavitating flow. Photographs and video movies are taken with stroboscopic light synchronized with the rotating shaft. These are used to observe the topology of vapor structures. Details of procedure and results are presented in Ref.

͓16͔.

A large range of flow rates was investigated around the nominal point of operation in noncavitating and cavitating conditions. Noncavitating experiments consist of setting a high static pressure level, to avoid any apparition of vapor and then to vary the flow rate. In cavitating conditions, the flow rate is kept constant, and the static pressure is decreased slowly to enhance vapor development in the inducer and to reach the performance breakdown.

In Fig. 3, breakdown curves are presented for three flow rates. Parameter ⌿ is evaluated from static pressures P inlet and P outlet measured at the shroud, as indicated in Fig. 2. ء is calculated as a function of the total pressure at the inlet level. The breakdown curves appear very gradual for this inducer geometry. Head drops equivalent to ϳ10% of the noncavitating performance correspond to ء ϳ 0.13 ͑for 1.25Qn͒, and to ء ϳ 0.11 ͑for 1.35Qn and 1.4Qn͒. In the final drop zone, corresponding to about 20% of the head drop and ء ϳ 0.07, each point is stabilized with Ϯ1% of the variation in the rate flow coefficient.

According to experimental observations ͓16͔, the cavitating zone is mainly located in a peripheral region and is fed by backflow circulation. From the spectral analysis of inlet pressures, and of axial and radial forces on the shaft, axial instability is observed during breakdown. No rotating cavitation signature appears on spectral cascades of the CREM1 inducer.

Flow Model and Numerical Tool

3.1 Governing Equations and the Barotropic Model. The applied two-phase flow model is based on a mixture assumption ͑one-fluid model͒. The compressible RANS system, written in a rotating frame, associated with the k-turbulent model ͓20͔, has been solved. The barotropic model proposed by Delannoy and Kueny ͓17͔ is applied to model the cavitation phenomenon. In the model, 1 and v are, respectively, the density of the pure liquid ͑corresponding to ␣ =0͒ and pure vapor ͑␣ =1͒. Both densities, corresponding to pure phases, are taken as constant. In the cavitation zone, the density of the vapor/liquid mixture is given as a function of the void ratio by the relation: This law is mainly controlled by its maximum slope, which is related to the minimum speed of sound c min in the mixture. In all computations presented in this paper, c min Ϸ 3m / s ͑Fig. 4͒. The influence of the law state parameter c min on the behavior of the inducers has been evaluated in previous works ͓11,20͔.

= ␣ v + ͑1-␣͒ l ͑1͒
The liquid density is set to 1000 kg/ m 3 and the vapor density is fixed to its physical value of 0.02 kg/ m 3 . Because of convergence difficulties in the case of the smallest cavitation numbers, some computations have been performed with a maximum value of v = 15 kg/ m 3 .

Numerical Methods.

Steady numerical simulations were carried out using the 3D CFD code FINE TM /TURBO ͓19͔, developed by Numeca International in collaboration with the LEGI laboratory. The discretization is based on a cell-centered finite-volume approach for multidomain structured meshes. The convective numerical fluxes are computed with a second order central scheme stabilized by the artificial dissipation proposed by Jameson et al.

͓21͔. An explicit four-stage Runge-Kutta time stepping procedure is used to advance the solution to the steady state. For incompressible areas, a low-speed preconditioner is introduced ͓18,23͔, enabling the reduction of the eigenvalue stiffness and the improvement of convergence.

The Yang-Shih k-turbulent model ͓22͔ with extended wall functions is used. The turbulence transport equations are integrated with a first order upwind scheme. The physical model and the numerical code are detailed in Refs. ͓11,23͔.

In the proposed cavitation model, thermal effects are neglected in the vaporization and condensation phenomena ͓24͔. Complementary works are in progress to improve the cavitation model in order to take into account thermodynamic effects ͓25,26͔.

Grids and Boundary Conditions.

Only one blade-toblade channel is considered and periodicity conditions are applied to its frontiers in azimuthal directions to simulate the presence of the contiguous blades. The mesh is composed of four 3D blocks ͑Fig. 5͒: the upstream part of the bulb, the blade-to-blade channel, and two blocks, respectively, upstream from the blade leading edge and downstream from its trailing edge. The whole mesh contains around 700,000 nodes for the first considered grid V1 ͑Fig. 1͑b͒͒. The tip leakage is not taken into account in the present study.

The y+ values from the center of the adjacent cells to walls vary from 0.8 to 100, for all computations, and the largest y+ values are mainly placed in the leading edge near the shroud of the inducer.

The calculations are performed with water flow conditions. The mass flow is imposed at the inlet, and the static pressure is imposed at the outlet. The hub and the blades are rotating, the shroud is immobile. Wall functions are imposed along solid boundaries.

In order to develop the phenomenon of cavitation in the numerical simulations and to obtain different values of ء during head breakdown, the imposed vapor pressure is increased progressively. The calculated inlet pressure is used to evaluate the numerical cavitation parameter ء . Parameter is calculated by the difference between the outlet and inlet static pressures, as indicated in Fig. 5.

In order to evaluate the grid influence, two additional grids ͑V2 and V3͒ have been tested. The grid specifications are shown in Table 1. Steady numerical simulations have been performed in cavitating conditions at 1.4Qn flow rate with the same conditions applied in the reference grid ͑V1͒. Numerical head drop charts are presented in Fig. 6. It is worth noting that each calculated steady state solution was carefully converged. The relative differences between noncavitating pump heads cal-culated with the three meshes are smaller than 5%. Concerning ء equivalent to a 10% head drop, the maximum relative error between calculated values is about 30%, which corresponds to a ⌬t ء ϳ 0.02. Numerical results have been analyzed in order to improve knowledge about head drop mechanisms, local flow characteristics, and secondary flows in the inducer. The development of cavitation structures ͑Sec. 4.4͒ calculated with these grids are quite similar, and the analyses of secondary flows, presented in Sec. 4.5, appear almost independent to the applied mesh.

The aim was also to propose a methodology of 3D calculations and analyses of complex cavitating flows. In relation to prior works presented in Ref. ͓11͔ concerning the inducer geometry, meaningful numerical work has been performed by Pouffary et al. ͓4͔,Numeca International,and Rolland ͓26͔ to improve the preconditioner and the stability of the numerical code for calculations of cavitating flows. It is worth noting that these numerical modifications led to major improvements in the accuracy and numerical stability in comparison to our previous computations ͓10͔.

Numerical Results and Discussion

Steady numerical simulations of the inducer have been performed with grid V1 in cavitating conditions at a given rotational speed and three flow rates: 1.4Qn, 1.35Qn, and 1.25Qn. First, we present a quantitative analysis by comparing the experimental and numerical results of head drop charts.

The power and efficiency evolutions obtained by numerical simulations are illustrated for the three flow rates. Some qualitative results, concerning the blade load and visualizations of vapor/ liquid structures are also presented for several ء values. Finally, a study of secondary flows during the ء decrease is proposed. For 1.4Qn and the flow rate, the numerical prediction at the 10% head drop were found for values of ء close to 0.1. For 1.35Qn and 1.25Qn flow rates, it is close to ϳ0.11, these estimations correspond well to experimental observations. The most important differences between the numerical and experimental values of / 0 ͑around ϳ4%͒ are observed in the range of 0.2 Ͻ ء Ͻ 0.4.

We underline that for lower flow rates, the angle of attack of the flow at the blade leading edge is higher than for 1.4Qn; the cavity Globally, the simulated head drops are gradual, and numerical results show a good correlation with the experimental ones for the three studied cases. The breakdown is clearly associated with a decrease in the torque and the hydraulic efficiency, as illustrated below.

Power and Efficiency Evolution.

As illustrated in Fig. 8, the computed power and efficiency of the inducer are calculated at considered flow rates. The power is done from the calculated numerical torque and the fixed rotational speed. The torque is computed by the integration of the pressure and the stress tensor along the blade surface. The Power ء represents the ratio of calculated power to reference power corresponding to the noncavitating condition for 1.4Qn. The efficiency is defined as the ratio of the provided power to the available power.

The maximal values of Power ء are observed for ء around 0.38, higher than the value in the noncavitating condition. This may be attributed to the change in the flow incidence angle due to the appearance of cavitation attached to the blade that induced a better blade load, mainly near the hub ͑see Sec. 4.3͒. For smaller ء values, the computed power decreases, this is attributed to the development of cavitation structures and a lower blade load.

The efficiency decreases with the cavitation parameter reduction. The efficiency variation is about 3% for 1.4Qn, 4% for 1.35Qn, and more than 5% for 1.25Qn in the range of analyzed ء . We conclude that the head drop is associated with both the decrease in torque and the hydraulic efficiency.

Static Pressure Distributions Around a Blade Section.

Figure 9 shows the different sections considered for the study of the static pressure distributions around the blade: midspan plane, two planes placed between the midspan and shroud, and one plane placed near the hub.

The blade load at midspan and three additional planes ͑h 1 , h 2 , and h 3 ͒ at 1.4Qn flow rate are drawn in Fig. 10 during the ء decrease. For these studied planes and analyzed ء values, it can be observed that the blades are loaded mainly near the leading edge. In the cavitating regime, the development of the sheet cavity from the leading edge of the suction side modifies progressively the blade load and, consequently, the inducer head decreases gradually, as illustrated in Fig. 7͑a͒.

The blade load is slightly modified for ء = 0.38. For the small- est ء value, the figures show clearly a short sheet cavity attached from the leading edge of the suction side of the blades. For the considered conditions, the appearance of the cavitation does not influence the pressure side load. We do not observe the strong influence of the cavitation on the trailing edge pressure distribution because the cavity attached is not developed toward the interior of the blade-to-blade channels for the considered ء values.

The pressure values for plane h 2 ͑near the shroud͒ at the pressure side are higher than the ones evaluated for the other planes. Moreover, the blade load seems to be more important at the shroud.

To improve the analyses during the ء decrease, the evolution of the blade load at additional plane h 3 near the hub ͑Fig. 10͑d͒͒ has been studied. For this plane, a better blade load for ء = 0.38 was observed, if compared with the noncavitating condition; the appearance of the cavitation influences the pressure side load and The attached cavitation sheet slightly increases from the hub to the shroud and merges with the cavitating backflow structures. Head drop occurs progressively with the increase in vapor structures near the shroud. It is worth noting that at ء ϳ 0.1, the vapor structures do not reach the throat and do not enter into the bladeto-blade channel.

The cavitation behavior obtained by numerical simulations presents a very good qualitative agreement with experimental visualizations ͑Figs. 12 and 13͒. We can observe that a consequence of the CREM1 design is the development of the ringed vapor structure at the shroud, which limits the development of the cavitation sheet on the blades.

Calculations performed with grids V2 and V3 lead to similar vapor/liquid structures.

Local Flow

Analyses. An analysis of the flow within the inducer has been performed in three different sections in the blade-to-blade channel ͑Fig. 14͒. The local analysis of the influence of cavitation on the flow structure in the different sections has been made using the convention in Fig. 15.

Helicity.

In the present paper, the analysis of the secondary flows is first based on the dimensionless relative helicity.

Secondary flows can be characterized by the dimensionless relative helicity, which is defined as The normalized relative helicity indicates the value of the cosine of the angle between the velocity vector and the swirl vector. Consequently, the swirling centers will be associated with values +1 and -1, the sign determining the direction of the rotation of the swirl. This quantity is drawn in Fig. 16 for the three defined sections during the breakdown, for the 1.4Qn flow rate. The black color represents anticlockwise vorticity and yellow to white indicate clockwise vorticity.

In Sec. 1, the secondary flow in the rotating frame develops mainly along the blade suction side ͑white color͒, especially in the noncavitating condition. With the cavitation coefficient decrease, the clockwise vorticity on the suction side tends to be less important. A clockwise vorticity structure of low intensity is located on the pressure side. This behavior is due to the development of the boundary layer, which is less important on the pressure side than on the suction side because of the leading edge proximity on the pressure side. It is also observed that the flow is disturbed in the zone near the hub; the flow displays a high clockwise vorticity.

It is important to recall that tip leakage was not considered in the numerical simulations ͑which explains the absence of zones with a strong vorticity near the shroud͒.

In the noncavitating regime in Sec. 2 ͑Fig. 16͑b͒͒, we observe a clockwise vorticity ͑white color͒ on the suction side. For the cavitating conditions the secondary flows are similar to the structures observed in Sec. 1.

Finally, in Sec. 3 ͑Fig. 16͑c͒͒, in the noncavitating regime, the clockwise vorticity is located on the suction side and some regions of the anticlockwise vorticity appear close to the shroud on the pressure side. The vorticity effects are less important near the hub for this section due to the meridian curvature of this inducer.

In the analyzed sections in cavitating conditions, we observed some vorticity structures ͑white color͒ located near the shroud, which seem to be associated with the ringed vapor structure at the shroud ͑Figs. 11 and 12͒. Boundary layers are developed in the zones close to solid surfaces of the blades in the pressure side and suction side for the three different sections.

Figure 17 shows the distribution of the radial velocity in the three sections. For noncavitating conditions we observe a strong centrifugation of the boundary layers, mainly at the suction side ͑radial velocity is around 12 m/s͒, which is associated to the secondary structures observed from helicity analyses.

It is worth noting that comparisons between results obtained from grids V1, V2, and V3 indicate slight modifications on calculated secondary flows, as illustrated in Fig. 18.

Axial

Velocity. In order to analyze the flow rate distribution in the inducer, results obtained for the 1.4Qn flow rate concerning the axial velocity in the three sections previously defined ͑Fig. 14͒ are illustrated in Fig. 19. Results correspond to different ء values. In Sec. 1 ͑Fig. 19͑a͒͒, a backflow zone near the shroud, characterized by negative axial velocities can mainly be observed. This zone is clearly observed in the meridian view by the analysis of the streamlines ͑Fig. 20͒. An anticlockwise vorticity ͑Fig. 16͑a͒͒ can also be observed in this zone ͑Sec. 1, noncavitating condition͒.

The figures corresponding to Sec. 2 ͑Fig. 19͑b͒͒ show that the zones of the negative axial velocity are mainly near the hub. For the smallest values of the cavitation parameter, the appearance of zones of low velocities close to the shroud can be observed. For Sec. 3 ͑Fig. 19͑c͒͒, the negative axial velocity area is placed in the pressure side near the hub and it seems to be related to the high meridian curvature of the hub.

In a general overview, it can be noticed that the flow rate is not uniformly distributed in the channel and that axial velocities are more important at the blade suction side. By analyzing results presented in Fig. 20, we observe that the backflow zone near the shroud tends to decrease during breakdown; the appearance of cavitation structures and the reduction of pressure gradients between the suction and pressure sides lead to a decrease in the thickness of the backflow area ͑in the meridian plane͒. The width ͑Fig. 19͑a͒͒ and the length ͑Fig. 20͒ of this zone are also modified by cavitation structures.

Conclusion

In the present work, 3D steady RANS computations were performed and compared with experiments concerning water cavitating flows in a three-blade inducer at different flow rates. In complement to prior works ͓4,10,11͔, the main targets of the present study were to pursue model validation and to improve analyses of complex 3D cavitating flows in rocket engine turbopump inducers.

The qualitative development of the cavitation in the inducer was simulated and the extension of vapor near the shroud was favorably compared with experimental results. Experimental and numerical performance charts concerning three different flow rates were also compared. Thanks to recent numerical improvements, predictions of the inducer breakdown were adequate and much better than the ones presented in previous works ͓11͔.

Analyses of the numerical results provided interesting information about the power, blade load, and hydraulic efficiency evolution during head breakdown. Secondary flows were analyzed from dimensionless relative helicity maps corresponding to three different sections in the blade-to-blade channel. Moreover, local analyses of axial velocity fields enabled the determination of backflow zones and of the flow rate distribution in the inducer channels. For the considered cavitation number range, the study pointed out a relevant influence of the cavitation structures on the flow fields around the leading edge, near the shroud zone.

The methodology of analysis proposed in this paper showed that numerical simulations can provide useful information for the design of turbomachinery, and more particularly for inducers. Based on the experimental and numerical results obtained, the inducer geometry was improved. A new inducer was designed and is being tested for future works concerning mainly instable cavitation analyses and prediction.

Complementary works are also in progress, mainly

• to improve the cavitation model to take into account thermodynamic effects, mainly in the case of cryogenic propellants ͓25,26͔ and • to implement and to apply a multigrid strategy to accelerate the convergence of cavitating flow calculations in inducer geometries. This work is very important for carrying out 3D unsteady calculations and for analyzing cavitation instabilities ͓8,20,27͔.

It is worth underlining that further local and global experimental studies are also required in order to validate and calibrate the physical models applied in an improved way.

Nomenclature

Fig. 1 Fig. 2

 12 Fig. 1 Details of the inducer: "a… photograph of the inducer CREM1 and "b… a view of the inducer meshing

Fig. 5

 5 Fig. 5 Inducer views: "a… Inducer meridian view. Representation of mesh blocks and boundary conditions. "b… Representation of inducer blade-to-blade view mesh blocks.

4. 1

 1 Head Drop. Numerical head drop charts are compared with experimental ones for three flow rates in Fig. 7. The absolute comparisons between the experimental and numerical parameters are not possible because outlet pressure values are not taken at the same points ͑see Figs. 2 and 5͒. The head drop / 0 ͑ ء ͒ is shown in relation to the maximal values ͑ 0 ͒ obtained for each case ͑experimental and numerical͒ and for the different flow rates.

Fig. 7

 7 Fig. 7 Experimental and numerical head drop charts. "a… 1.4Qn, "b… 1.35Qn, and "c… 1.25Qn.

Fig. 9

 9 Fig. 9 Planes h 3 , h 2 , and h 1 and the midspan represented in the meridian view considered for the blade load analysis

Fig. 11

 11 Fig. 11 Isoline of density " = 950 kg/ m 3 … during the ء de- crease at 1.4Qn; 1 is no cavitation, 2 is ء = 0.31, 3 is ء = 0.14, and 4 is ء = 0.11

Fig. 16

 16 Fig. 16 Secondary flows: dimensionless relative helicity representation at three locations in the blade-to-blade channel: "a… Sec. 1, "b… Sec. 2, and "c… Sec. 3. Blue represents the maximum anticlockwise vorticity zone and red the maximum clockwise vorticity zone. The results are shown for the noncavitating condition "1… and typical cavitating regimes: ء = 0.31, ء = 0.14, and ء = 0.11 "2, 3, and 4….

c

  min ϭ minimum speed of sound in the vapor/liquid mixture ͑m/s͒ L ref ϭ geometry reference length ͑m͒ P ϭ local static pressure ͑Pa͒ P tot ϭ P + ͑V 2 / 2͒ local total pressure ͑Pa͒ P vap ϭ vapor pressure ͑Pa͒ power ϭ T ͑W͒ power ء ϭ power/ power noncavitating dimensionless power Q ϭ flow rate ͑m 3 / s͒ Qn ϭ nominal flow rate ͑m 3 / s͒ T ϭ torque ͑Nm͒ V ϭ flow velocity ͑m/s͒ V ref ϭ reference velocity: tangential tip velocity ͑m/s͒ W ϭ relative flow velocity ͑m/s͒ ⌬P tot ϭ P tot downstream -P tot upstream total pressure variation in inducer ͑Pa͒ ⌿ ϭ ͑P static outlet -P static inlet ͒ / • V ref 2 shroud elevation coefficient ⌿ 0 ϭ shroud elevation coefficient for the noncavitating condition ϭ fluid density ͑kg/ m 3 ͒ l , v ϭ liquid/vapor density ͑kg/ m 3 ͒ a ϭ ͑ -l ͒ / ͑ v -l ͒ void ratio ϭ⌬ P tot Q / T hydraulic efficiency ϭ P tot inlet -P vap / • V ref 2 cavitation parameter ء ϭ ͑ / const͒ dimensionless cavitation parameter w ϭ angular velocity ͑rad/s͒ y + ϭ dimensionless distance to the wall

Fig. 17

 17 Fig. 17 Secondary flows: radial velocity evolution in the noncavitating regime in Secs. 1-3 for 1.4Qn

Fig. 19

 19 Fig. 19 Secondary flows: axial velocity distribution in the blade-to-blade channel: "a… Sec. 1, "b… Sec. 2, and "c… Sec. 3.The results are presented for the noncavitating regime "1… and cavitating regimes ء = 0.31, ء = 0.14, and ء = 0.11"2, 3, and 4….

Fig. 20

 20 Fig. 20 Streamlines in the meridian plane: "a… noncavitating regime, "b… cavitating regime corresponding to ء = 0.31, "c… cavitating regime corresponding to ء = 0.14, and "d… cavitating regime corresponding to ء = 0.11

Table 1 Tested meshes

 1

		Reference grid V1	V2	V3
	Totals numbers of grid point	697564	838672	962844
	y + max	97.6	87.2	78.3
	y + min	0.81	0.65	0.61

Acknowledgment

The authors wish to express their gratitude to French space agency CNES, to Snecma for their continuous support, and to Numeca International for their cooperation in the development of the numerical code. This research was also supported by a doctoral grant from the Los Andes University, Venezuela.