F Tabkhi

L Pibouleau

C Azzaro-Pantel
email: catherine.azzaropantel@ensiacet.fr

S Domenech

Keywords: natural gas, network, pipeline, compressor, optimization, MINLP

Total Cost Minimization of a High-Pressure Natural Gas Network

This paper deals with a high-pressure gas pipeline optimization, where the problem is to find the design properties of the pipelines and necessary compressor stations to satisfy customer requirements, using available supply gas and storage capacities. The considered objective function is the total annualized cost, including the investment and operating costs. The binary variables used to represent the flow direction of pipelines lead to a mixed integer nonlinear programming problem, solved by using the standard branch and bound solver in GAMS. The optimization strategy provides the main design parameters of the pipelines (diameters, pressures, and flow rates) and the characteristics of compressor stations (location, suction pressure, pressure ratio, station throughput, fuel consumption, and station power consumption) to satisfy customer requirements.

Introduction

Natural gas ͑NG͒, viewed as a cleaner-burning alternative to coal and oil in terms of acidic and greenhouse gas pollution, is increasingly being used as an energy source and by most estimates, its global consumption will double by 2030 ͓1͔. The transport of large quantities of NG is carried out by pipeline network systems across long distances. For example, the European natural gas system, which is very well developed, consists, inter alia, of 1.4ϫ 10 6 km pipelines of which 145,000 km concern highpressure transmission pipelines. As the gas flows through the network, pressure ͑and energy͒ is lost due to both friction between the gas and the pipe inner wall, and heat transfer between the gas and its environment. Typically, natural gas compressor stations are located at regular intervals along the pipeline to boost the pressure lost through the friction of the natural gas moving through the steel pipe. They consume a part of the transported gas, thus resulting in an important fuel consumption cost on the one hand, and in a significant contribution to CO 2 emissions, on the other hand. Two main issues are generally highlighted when considering pipeline transmission networks, i.e., designing and operating a gas pipeline network. They generally involve technico-economic concerns, often based on capital cost minimization, throughput maximization, and fuel cost minimization.

Numerical simulations based on either steady state or transient models of the networks have been used to attempt to provide solutions to these problems. Following the development of the gas industry, gas pipeline networks have evolved over decades into very large and complex systems. A typical network today might consist of thousands of pipes, dozens of stations, and many other devices, such as valves and regulators. Inside each station, there can be several groups of compressor units of various vintages that were installed as the capacity of the system expanded. In that context, numerical simulation and optimization of gas pipeline can be of great help to design them, to predict their behavior, and to control their operation. Over the years, many researchers have attempted this issue.

The difficulties of such optimization problems arise from several aspects. First, compressor stations are very sophisticated entities themselves. They might consist of a few dozen compressor units with different configurations and characteristics. Each unit could be turned on or off, and its behavior is nonlinear. Second, the set of constraints that define feasible operating conditions in the compressors along with the constraints in the pipes constitutes a very complex system of nonlinear constraints. Third, the considered optimization problem often involves both continuous and integer variables and gives birth to a mixed integer nonlinear programming ͑MINLP͒ formulation. Moreover, the treatment of industrial-size problems may render the problem highly combinatorial.

The natural gas chain is generally constituted by various components as represented in Fig. 1. In this figure, the pressure regimes are just indicative and may differ from country to country.

The transport lying system between the natural gas deposits and the consumers is quite complex. After the gas has been extracted, so-called trunk lines are connected with pipeline compressor stations. The natural gas is then pumped into long distance pipelines called transmission lines and sent to the take-off stations for the consumers. From there, the gas is further transported to the control station of the regional distribution system. It then finally goes to industrial customers and households. A schematic view of a pipeline section is displayed in Fig. 2 with six compression stations, delivery, and supply points.

For transmission network problems considered in this study, compressor stations in a pipeline system can be subdivided into two classes: the originating stations, which are positioned at the inlet to the pipeline and are usually the most complex, and the booster stations, which are located along the pipeline to compensate for the pressure decrease due to friction and elevation losses. In principle, the longer the pipeline and the elevation of the terrain crossed, the more compressor horsepower is required to achieve the required delivery pressure at destination. However, under a fixed route and flow capacity, the number and size of booster stations can vary depending on circumstances and design. Although systems with fewer stations can be easier to operate, they have the disadvantage of introducing a need for high inlet pressures. Actual transmission systems represent a compromise between a very few powerful originating stations and a large number of small booster points. This paper deals with a high-pressure gas pipeline transmission network optimization, where the problem is to find the design properties of the pipelines and necessary compressor stations to satisfy customer requirements, using available supply gas and storage capacities. The considered objective function is the total annualized cost, including the investment and operating costs. The binary variables used to represent the flow direction of pipelines lead to a MINLP problem, solved by using the GAMS package. The approach is illustrated by an example taken from literature ͓2͔ inspired from the Belgian network. It comprises 20 nodes linked together with 20 arcs ͑see Fig. 3͒. The positions of nodes are selected such that they are consistent with their geographical location in the real network.

Previous Works

2.1 Transmission Pipeline Modeling. The qualitative presentation of the pipeline network and of its components highlights that engineering design and operation studies of pipeline network are involved with the optimization of the performance of the global system. The analysis of the dedicated literature shows that there is growing interest on the subject.

The optimization of the design of a pipeline to transmit fluids involves numerous variables, which include pipe diameter, pressure, temperature, line length, space between pumping or compressor stations, required inlet and delivery pressures, and delivery quantity. Each of these parameters influences the overall construction and operating cost in some degree. This is as true for the design of a system from a clean sheet of paper ͑grass roots͒ as it is for the development and upgrading of an existing system; the only real difference between these two examples is the extent to which some of the variables are already fixed.

Because of the number of variables involved, the task of establishing the optimum can be quite difficult, and in order to ensure a robust solution, many options may have to be investigated.

The two main approaches that are classically encountered in gas network representation are numerical simulation and optimization.

The main purpose of simulation is to determine the actual behavior of a gas network under given conditions. Simulation basically answers the following question: What happens if we run our grid with given control variables and known boundary flows? Typical questions like finding a control regime, which achieves several target values, usually require a series of simulation runs by expert users who are familiar with the network. Two disadvantages of numerical simulation are found. First, finding an adequate regime may even take a large number of runs and, second, it cannot ensure that the solution achieved is optimal.

This explains mainly why the searching process must be substituted with more sophisticated algorithms. Yet, optimization generally works with simplified models, but it yields optimum results where limits or certain target values will be achieved automatically if they are defined as optimization problem constraints. There has been a great deal of work on the optimization approach to gas pipe distribution networks.

Since 20 years, there has been an interest on the optimization of gas pipe distribution networks. A one-dimensional compressible fluid flow equation is proposed in Ref. ͓3͔. Lewandowski ͓4͔ implemented an object-oriented methodology for modeling a natural gas transmission network using a library of C ++ classes, and Osiadacz ͓5͔ presented a dynamic optimization of highpressure gas networks using hierarchical system theory. Surry et al. ͓6͔ formulated the optimization problem based on a multiobjective genetic algorithm. Mohitpour et al. ͓7͔ used a dynamic simulation approach for the design and optimization of pipeline ͓10͔ based their modeling approach on a hybrid network using being at the academic level, some ͑commercial or free͒ computational codes are available: the standard branch and bound ͑SBB͒, BARON, DICOPT+ +, and LOGMIP solvers within the GAMS modeling environment ͓40͔, and MINLP_BB ͓41͔ and aECP ͓42͔.

The second class, namely, metaheuristics or stochastic methods, is based on the evaluation of the objective function at different points of the search space. These points are chosen through the use of a set of heuristics, combined with generations of random numbers. Thus, metaheuristics cannot guarantee to obtain an optimum. They are divided into neighborhood techniques such as simulated annealing ͓43͔, tabu search ͓21͔, and evolutionary algorithms comprising genetic algorithms ͓44͔, evolutionary strategies ͓45͔, and evolutionary programming ͓46͔.

Since the number of constraints associated with the formulation of the problem related to the optimization of NG transmission networks may be important, the deterministic way is adopted in this study. To represent the deterministic class, solvers of the GAMS environment were chosen, since this optimization tool is widely used, and even stands as a reference for the solution of problems drawn from process engineering.

3 Model Formulation for Gas Pipeline Networks 3.1 Gas Pipeline Hydraulics. Due to operating problems, a gas transmission line is not normally designed to handle a twophase flow. Exceptions lie, for example, in oil/gas wells, gathering systems, and separation units. The formulation presented here is only valid for a single phase flow. The pressure drop in a gas pipeline, i.e., the essential parameter to determine the required compression power for the transmission, is derived from the differential momentum balance.

Writing momentum balance, the equation of movement is obtained as ͓47͔:

ץp ץx + f 2D ZRT pMA 2 m 2 Ϯ g pM ZRT sin ␣ + R A 2 M ץ ץx ͩ ZT p m 2 ͪ+ 1 A ץm ץt = 0

͑1͒

where A is the cross area of the pipe ͑m 2 ͒, R is the universal gas constant ͑8314 J/kmol K͒, T is temperature ͑K͒, m is the flow rate ͑kg/s͒, and ␣ is the angle between the pipe axis and the horizontal line. By using Eqs. ͑3͒-͑9͒ and integrating Eq. ͑1͒, Tabkhi ͓48͔ showed that between points 1 and 2 of a pipe of length L ͑m͒ and internal diameter D ͑m͒, the steady-state pressure drop can be expressed in the following form:

͑p 2 2 -p 1 2 ͒ - 32m 2 ZRT 2 D 4 M lnͩ p 2 p 1 ͪ+ 16f 2 D 5 ZRT M m 2 L = 0 ͑2͒
In this equation the various terms are defined as follows.

Gas density and pressure are linked by the following equation containing the compressibility factor Z:

= pM ZRT ͑3͒
where the average molecular mass of gas M is calculated using a simple mixing rule involving the mole fractions y i s and the molecular masses M i s:

M = ͚ M i y i ͑4͒
The compressibility factor Z is used to alter the ideal gas equation to account for the real gas behavior. Traditionally, the compressibility factor is calculated using an equation of state. Yet, for natural gas, it may be estimated from the empirical relationship proposed for simulation goals in literature ͓49͔. For example, this factor can be expressed as a function of the critical properties of the gas mixture, average pressure of the pipe segment, and the temperature that have been considered as constant: This literature analysis shows that there has been and continue to be a significant effort focused on the modeling of natural gas transmission networks. The objective of this work is to propose a framework able to embed formulations from design to operational purposes: this explains why only steady-state behavior of the gas flow is considered. The problem is to implement, for a given mathematical model of a pipeline network, a numerical method that meets the criteria of accuracy together with relatively small computation times.

Optimization Techniques.

A great variety of applications, drawn from a wide range of investigation areas, can be formulated as complex optimization problems. It covers, for instance, the famous traveling man problem intensively studied ͓18͔ as well as frequencies allocation for radio-mobile networks ͓19͔, process network optimization ͓20͔, physicochemical equilibrium calculations ͓21͔, or hydrology computing ͓22͔. In order to face these problems, a significant investigation effort has been carried out to develop efficient and robust optimization methods. At the beginning, this aim was concerned especially in the operational research and artificial intelligence areas. But the trend was subsequently followed by the process system engineering community, since this one provides a wide number of applications formulated as complex optimization problems. A typical reference is constituted by design problems related to heat or mass exchanger networks ͓23͔, supply chain design ͓24͔, multiproduct ͓25͔ or multipurpose ͓26͔, and batch plant design or retrofitting ͓27͔.

As a consequence, a great diversity of optimization methods was implemented to meet the industrial stakes and provide competitive results. But if they prove to be well-fitted to the particular case they purchase, these techniques' performance cannot be constant whatever the treated problem is. Actually, method efficiency for a particular example is hardly predictable, and the only certainty we have is expressed by the no free lunch theory ͓28͔: There is no method that outdoes all the other ones for any considered problem. This feature generates a common lack of explanation concerning the use of a method for the solution of a particular example, and usually, no relevant justification for its choice is given a priori. Among the diversity of optimization techniques, two important classes have to be distinguished: deterministic methods and stochastic ones. Complete reviews are proposed in literature ͓29-31͔. A thorough analysis of both classes was previously studied in Ref. ͓32͔ with the support of batch plant design problems.

The deterministic methods involve the verification of mathematical properties of the objective function and constraints, such as continuity or derivability. This working mode enables them to ensure to get an optimum, which is a great advantage. Among the deterministic class, the following ones stand out: the branch and bound methods ͓33-35͔, the generalized Bender decomposition ͓36͔ and the outer approximation algorithms ͓37͔, the extended cutting plane method ͓38͔, and disjunctive programming methods ͓39͔. Even though most of the above-mentioned methods keep

Z = 1 + ͩ0. 257 -0.533 T c T ͪ p ij p c ͑5͒
The pseudocritical temperature of natural gas T c and its pseudocritical pressure p c can be calculated using an adequate mixing rule starting from the critical properties of the natural gas components. The critical point of a material is the point where the distinction between the liquid and vapor phases disappears. In this work, average pseudocritical properties of the gas are determined from the given mole fractions of its components by Kay's rule, which is a simple linear mixing rule shown in:

T c = ͚ T ci y i ͑6͒ p c = ͚ p ci y i

͑7͒

Average pressure p ij can be calculated from two end pressures ͓49͔

p ij = 2 3 ͩpi+pj-p i p j p i + p j ͪ ͑8͒
Since the flow is considered fully developed here, which is the case concerned in gas pipelines, the friction factor f is estimated through the equation deduced by Prandtl-von Karman ͓50͔ in which the friction factor depends only on the relative roughness as follows:

f = ͩ-2 log /D 3.71 ͪ -2 ͑9͒

Fundamental Constraints

MAOP.

The internal pressure in a pipe causes the pipe wall to be stressed, and if allowed to reach the yield strength of the pipe material, it could cause permanent deformation of the pipe and ultimate failure. In addition to the internal pressure due to gas flowing through the pipe, the pipe might also be subjected to external pressure, which can result from the weight of the soil above the pipe in a buried pipeline and also by the probable loads transmitted from vehicular traffic. Therefore, the necessary minimum wall thickness will be dictated by the internal pressure in a gas pipeline. The pressure at all points of the pipeline should be less than the maximum allowable operating pressure ͑MAOP͒, which is a design parameter in the pipeline engineering ͑see Fig. 4͒. This upper limit is calculated using Eq. ͑11͒ ͓48͔ as follows:

p Ͻ MAOP ͑10͒ MAOP = SMYS 2t D -t f F f E f T ͑11͒
tor f F has been named the design factor. This factor is usually 0.72 for cross-country or offshore gas pipelines, but can be as low as 0.4, depending on class location and type of construction. The class locations, in turn, depend on the population density in the vicinity of the pipeline. The seam joint factor f E varies with the type of pipe material and joint type. Seam joint factors are between 1 and 0.6 for the most commonly used material types. The temperature factor f T is equal to 1 for gas temperatures below 120°C before it reaches 0.867 at 230°C.

Critical Velocities.

The increase in flow rate causes pressure drop increase. An important factor in the treatment of compressible fluid flow is the so-called critical flow. For a compressible flow, the increase in flow owing to the pressure drop increase is limited to the velocity of sound in the fluid, i.e., the critical velocity. Sonic or critical velocity is the maximum velocity which a compressible fluid can reach in a pipe. For trouble-free operation, velocities maintain under a half of sonic velocity. Sonic velocity in a gas c is calculated with a satisfactory approximation using Eq. ͑13͒. Here is the average isentropic exponent of the gas. C p is the heat capacity at constant pressure in J/͑kmol K͒.

v Ͻ c/2 ͑12͒ c = ͱ ZRT M ͑13͒ = ͚͑C pi y i ͒ M͚͑C pi y i ͒ -R ͑14͒
Increasing gas velocity in a pipeline can have a particular effect on the level of vibration and increase the noises too. Moreover, higher velocities in the course of a long period of time will cause the erosion of the inside surface of the tubes, elbows, and other joints. The upper limit of the velocity range should be such that erosion-corrosion cavitations or impingement attack will be minimal. The upper limit of the gas velocity for the design purposes is usually computed empirically with the following equation ͓51͔. In pipeline design domain, the erosional velocity v e falls always underneath the speed of sound in the gas:

v Ͻ v e ͑15͒ v e = 122ͱ
ZRT PM ͑16͒

Compressor Characteristics.

As shown in Fig. 5, a centrifugal gas compressor is characterized by means of its delivered flow rate and its pressure ratio, and the ratio between suction side pressure of the compressor and its discharge pressure. The compression process in a centrifugal compressor can be well formulated using isentropic process aiming for calculating horsepower The yield stress used in Eq. ͑11͒ is called the specified minimum yield strength ͑SMYS͒ of pipe material. SMYS is a mechanical property of the construction material of the gas pipeline. The fac-for a compressor station. The pressure ratio of a centrifugal compressor is usually linked with a specific term named "head" carried over from pump design nomenclature and expressed in meter even for compressors. The head developed by the compressor is defined as the amount of energy supplied to the gas per unit mass of gas.

The equation for power calculation can be expressed as follows:

P = m c h i i ͑17͒
where m c ͑kg/s͒ is the mass flow rate of compressed gas, h i ͑m͒ is the compressor isentropic head, and i is the compressor isentropic efficiency. This equation is obtained by considering compression adiabatic process that is a reasonable assumption because the heat transfer between gas and the outside is very low. For adiabatic compressor first the adiabatic efficiency is defined as:

i = P ideal P ͑18͒
As shown in the following equation, considering adiabatic compression, head is an index of the pressure ratio across the compressor. In this equation, p d is the discharge pressure of the compressor and p s is the suction pressure and is isentropic exponent and will be calculated using Eq. ͑14͒. The compressibility factor and the temperature here are considered at the suction side of the compressor ͓52͔

h i = Z s RT s M -1 ͫͩ p d p s ͪ ͑-1͒/ -1 ͬ ͑19͒
In this work, centrifugal compressors in the station are assumed to be driven by turbines whose supply energy is provided from a line of the gas derived from the pipeline passed through the station in order to be compressed as shown in Fig. 6. The flow rate of the consumed gas as fuel for the compression process in each compressor is obtained by dividing the required power for compression P by the mechanical efficiency m , driver efficiency d , and low heating value ͑LHV͒:

m f = 10 6 m c h i i m d LHV

͑20͒

LHV represents the quantity of energy released by mass unity of the gas during complete combustion. It is considered at 25°C and 1 bar in kJ/kg and is calculated from the mass lower heating values, LHV i of the molecules composing the gas:

LHV = ͚y i M i LHV i ͚y i M i ͑21͒ h i 2 = b 1 + b 2 Q s + b 3 ͩ Q s ͪ 2

͑22͒

Also, contours of constant isentropic efficiency could be fitted in the polynomial form of second degree shown in:

i = b 4 + b 5 Q s + b 6 ͩ Q s ͪ 2

͑23͒

The rotation speed of all compressors is comprised between the lower and upper boundaries as represented below. To prevent from surge phenomenon, by considering surge margin surge , the following constraint is introduced ͓54͔:

l Յ Յ u ͑24͒ surge Յ Q s -Q surge Q s ͑25͒
There is a surge flow rate Q surge corresponding to each compressor rotational speed ͓55͔:

Q surge = b 7ͩͩ Z s RT s Mp s 2 -1 h surge + ͩ Z s RT s p s M ͪ 2 ͪ /͑-1͒ -ͩ Z s RT s p s M ͪ 2 ͪ 1/2

͑26͒

In this equation, h surge is the surge head at a specified compressor speed and can be calculated using the following equation:

h surge 2 = b 1 + b 2 Q surge + b 3 ͩ Q surge ͪ 2

͑27͒

Considering a fixed value for the surge pseudo-efficiency, it will be introduced as a parameter during optimization procedure. Equation ͑27͒ represents a nonlinear correlation between surge flow rate and rotational speed of the compressor.

Incidence Matrices.

The different links between the elementary sections of a network can be defined using incidence matrices. So, all the relation between the variables of the system such as the material balances at steady state around the nodes of a pipeline network can be expressed under the very concise forms by using different types of incidence matrices such as the arc-node matrix ͓56͔. In the model, each pipe, each compressor, and each fuel stream are represented by an arc. Consider a network with n nodes, l pipe arcs, and m compressor arcs. Therefore, there will be m fuel streams since for each compressor unit there is a stream that carries fuel to it. Because in a compressor, compression process is carried out, a compressor unit can be named an active arc. In this way, a pipe segment, in where the pressure decreases, may be called a passive arc. Let us note that the fuel streams have been considered as inert arcs regarding pressure change through them. A flow direction is assigned preliminarily to each pipe that can be or not coincide with the real flow direction of the gas that runs through the arc.

Let A be a matrix with the dimension of n ϫ ͑l + m͒. Each of its elements, a ij is given by:

a ij = Ά 1 if arc j comes out from node i -1 if arc j goes into node i 0 otherwise

·

A is called the node-arc incidence matrix. Similarly, let B be another matrix with the dimension of l ϫ m whose element b ij is defined below and is named the pipe-compressor incidence matrix:

m f m f + m c m c
Centrifugal compressor Incorporated turbine C1

Fig. 6 Representation of the centrifugal compressor and its incorporated turbine

Applying standard polynomial curve-fit procedures for each compressor, the normalized head h i / 2 can thus be obtained under the form of the following equation ͓53͔:

b ij = Ά 1 if pipe i is connected to discharge node of compressor j -1 if pipe i is connected to suction node of compressor j 0 otherwise

·

The third defined matrix is the node-fuel incidence matrix, which describes the existing fuel stream derivations from a node, and it is called the compressor-fuel matrix. The dimension of this matrix is n ϫ m and its elements are defined as:

c ij = ͭ 1 if fuel stream i be derived from node j 0 otherwise

ͮ

This matrix is also defined as previously, indicating which fuel stream belongs to which compressor. These incidence matrices are used to write the material balances around each node i, the flow rate of the consumed gas as fuel for the compression process in each compressor and the equation of movement. For example, the material balance around node i is expressed as Eq. ͑28͒. In this equation S i represents the gas delivery or supply relative to this node. It is negative if the node is a delivery node and positive for a supply node where the gas is injected to the node. The network was already presented in Fig. 3. The properties of flowing gas through the network are given in Table 1. There are two supply points 1 and 8 to satisfy the delivery requirements. The daily gas inputs at supply nodes 1 and 8 should not exceed their lower and upper bounds as shown in Table 2. The pressure of the gas supplied from node 1 is not limited inferiorly.

Available gas flows from the nodes where gas can be injected from storage tanks listed in Table 3. Storage tanks are located in nodes 2, 5, 13, and 14. Each storage facility has a maximum daily delivery capacity to gas network reported in Table 3. Daily gas demand at the nine delivery nodes is given in Table 4. It is assumed that the mass flow rates of gas demands are fixed. Pressure boundaries imposed to deliver gases should be respected. There is no lower limit on the pressure related to interconnection nodes. The upper limits on these pressures are given in Table 5.

In addition, the pressure at all nodes including interconnection nodes is limited superiorly with the MAOP of the arcs, which link them together. In Table 6, upstream and downstream nodes for all arcs of the network are defined. For each arc, at the beginning of the optimization procedure, it is assumed that gas flows from the upstream node to the downstream node. Then, the direction of the flow may be changed by the optimizer.

MINLP Problem

General Formulation.

The goal is to search for the design properties of the pipelines and necessary compressor stations of the network presented in Fig. 3 to satisfy customer requirements using available supply gas and storage capacities. Yet, the topology of the transmission network is set, i.e., the links between two nodes and the distance of the corresponding arc. However, a degree of freedom is introduced in the procedure: No gas flow direction is assigned a priori. The pipe diameters, as well as the power of required compressor stations, are decided via a MINLP algorithm within GAMS modeling environment. In addition, locations of necessary stations are also determined by the optimizer.

Basic Design

Conditions. Design temperature is assumed equal to 281 K. All the other required properties of the components of natural gas are given in Table 1. For the used steel, SMYS, which is needed to determine pipeline MAOP, is supposed to be equal to 2000 bars that is within conventional typical values. Design factor, another parameter to calculate MAOP via Eq. ͑11͒ considered equal to 0.4, is a low enough and, consequently, safe value. Note that this factor depends on the class location of the pipeline, which depends also on the population density around it. Seam joint factor, which represents weld efficiency, is assumed to be equal to 1. The last factor present in Eq. ͑11͒, temperature derating factor, is considered equal to 1 since temperature is lower than 120°C. Absolute roughness of the interior surface of pipes is about 50 m for steel pipes.

In order to compute fuel consumption, the total efficiency of compressor stations is considered equal to 25% for all stations. It is assumed that a compressor station is at least 1 km far from the neighboring nodes.

Optimization Variables.

The continuous optimization variables are the pipeline diameters, the location of the compressor stations, as well as their horsepower, pressures, and gas flow rates at different nodes, pressure ratios, and gas consumption at compressor stations. In addition, a set of binary variables is used to define the flow directions in the pipelines. All of the continuous variables are positive except the flow rates of gas feeding or leaving the nodes of the network. This operating variable is negative for demand nodes, positive for supply or storage nodes, and zero for interconnection nodes.

There are basic design variables that are related to the equipment size including pipe diameters, the location of the compressor stations, as well as their horsepower. Note that for each pipeline, if a compressor station must be considered on a line, the diameter remains unchanged after the compressor station.

The other kind of continuous variables consists of operating variables such as pressures at different nodes, in addition to the gas flow rates feeding or leaving the nodes. None of the pressures subjected to nodes is fixed and all of them are bounded inferiorly and superiorly, so they must be decided during the optimization process.

Like as the mass flow rates through arcs, there are other operating variables such as pressure ratio at compressor stations, which are dependent variables that are computed when the previously nominated basic variables are fixed. Gas consumption at each station is another dependent variable. Because its amount is low enough compared with compressor throughput, the pipeline diameter will not be changed after a compressor station. Average gas velocity through a pipe, as well as its average compressibility factor, is also dependent on other variables such as the pressures at its end-point nodes.

Finally, a set of binary variables is defined to assign flow direction of pipelines. So a MINLP optimization problem will be encountered. If a binary variable would be equal to unity, flow direction will coincide with the one reported in Table 6, and if it would be equal to zero, gas flows in the direction opposite to that the one given in this table.

Objective Function.

The goal in this optimization problem is to minimize the total annual cost. The objective function, total annualized cost as given in Eq. ͑29͒ ͓57͔, includes the sum of the investment cost depreciated over 10 years and the operating annual cost for all arcs:

ATC = ͚ jarcs ͑ICP j + ICS j + OCS j ͒ ͑ 29͒
In this equation ATC is the annualized total cost expressed in euros/year. The investment cost consists of that of the pipelines named ICP and that of the required compressor stations named ICS. The investment cost corresponding to the pipelines depends obviously on their lengths and diameters as shown in Eq. ͑30͒ where C s , d, and L are yearly unit pipe capital cost, pipe diameter, and pipe length, respectively. C s is assumed 15,778 € / ͑km m year͒ as follows:

ICP j = ͑C s dL͒ j ͑30͒
The annualized investment cost of the compressor stations is divided into fixed and variable costs as shown in Eq. ͑31͒. The former term of the right-hand side is independent of station horsepower P, and takes into account installation costs, civil engineering, etc. The last term increases linearly WITH horsepower:

ICS j = ͑C f sgn͑P͒ + C b P͒ j ͑31͒
The terms C f and C b in Eq. ͑31͒ are fixed and variable unit capital costs, respectively, and are supposed equal to 7410 € / year and 70 € / kW in this order ͓57͔. If no compressor exists on an arc, its corresponding horsepower and so the term ICS j are equal to zero. The sign function is introduced in the formulation to take into account that for arc j, if the compressor station horsepower P is equal to zero, then ICS is also equal to zero. The only operating and maintenance cost that comes into account is due to the compressor stations. It is supposed that pipelines have no operating costs compared with compressor stations ͓57͔. The yearly operating cost of a compressor station is also a OCS j = ͑C o P͒ j ͑32͒ 4.2.5 Constraint Definition. Calculation of MAOP using Eq. ͑11͒ for each pipeline requires its wall thickness that is obtained from the following equation in meter. This equation is obtained using the scheduled dimensions provided by ASME B36.19M standard that concerns stainless steel pipes: t = 52 ϫ 10 -3 d + 989 ϫ 10 -7 ͑33͒

For each pipeline, using Eq. ͑2͒, continuity equation corresponding to the segment before compressor is written in the form of Eq. ͑34͒ by introducing the binary variable b that indicates the flow direction.

2 MD 5 16ZRTfL c ͑2b -1͒͑P i 2 -P 1 2 ͒ = ͓bm ˙2 + ͑1 -b͒͑m ˙-m ˙c͒ 2 ͔ ͫ 1 + 2D fL c lnͩ p i p 1 ͪͬ ͑34͒
The friction factor is calculated using Eq. ͑9͒. Similarly, continuity equation for the section of pipeline after a compressor station is written in the form of:

2 MD 5 16ZRTf͑L -L c ͒ ͑2b -1͒͑P 2 2 -P j 2 ͒ = ͓͑1 -b͒m ˙2 + b͑m ˙-m ˙c͒ 2 ͔ ͫ 1 + 2D f͑L -L c ͒ lnͩ p 2 p j ͪͬ ͑35͒
Using Eq. ͑28͒, the material balance around node i can be written as Eq. ͑36͒. In this equation S i represents the gas delivery or supply to this node. It is negative if the node is a delivery node, positive for a supply or storage node where the gas is injected to the gas network, and zero otherwise.

͚ karcs a i,j m j ͑2b j -1͒ -m cj ͑a i,j ͑2b j -1͒ -1͒/2 = S i ͑36͒

The existence of a compressor station on an arc is decided when pressure ratio is more than 1, its lower bound. Pressure ratio is obtained using following equation where p 1 and p 2 are compressor station end-point pressures. Note that if b = 1 then p 1 and p 2 are suction and discharge pressures, respectively, and if b = 0 vice versa:

Ratio = bͩ p 2 p 1 ͪ+͑1-b͒ͩ p 1 p 2 ͪ, 1Յ ratio Յ 2 ͑37͒ P = 0 ∨ , P -Յ P Յ P + , P -0 ͑38͒
Finally, the average gas velocity through pipelines is restricted by erosional velocity obtained following Eq. ͑16͒.

4.2.6 Optimization Procedure. Following the study of Ponsich ͓32͔, the code SBB of GAMS was chosen as a solution procedure. SBB is a GAMS solver based on a combination of the standard branch and bound method known from mixed integer linear programming ͑MILP͒ and some of the standard NLP solvers already supported by GAMS such as CONOPT and sparse nonlinear optimizer ͑SNOPT͒. SBB supports all the types of discrete variables including binary, integer, and some other specific ones ͓40͔.

The relaxed mixed integer nonlinear programming ͑RMINLP͒ model is initially solved using the starting point. A RMINLP model can contain both discrete variables and general nonlinear terms, but the discrete requirements are relaxed. In other words, the integer and binary variables can assume any values between their bounds; thus a RMINLP model will be just as a NLP problem. If all discrete variables in the RMINLP model are integer, SBB will return this solution as the optimal integer solution. Otherwise, the current solution is stored and the branch and bound procedure will start.

During the branch and bound process, the feasible region for the discrete variables is subdivided, and bounds on discrete variables are tightened to new integer values to cut off the current noninteger solutions. Each time a bound is tightened, a new tighter NLP submodel is solved starting from the optimal solution to the previous looser submodel. The objective function values from the NLP submodel are assumed to be lower bounds on the objective in the restricted feasible space ͑assuming minimization͒, even though the local optimum found by the NLP solver may not be a global optimum.

First Optimization

Case. In this first case, pipeline diameters are considered as continuous variables. Then the resulted diameters are rounded up toward commercial sizes presented in Table 7 ͑commercial pipe sizes are based on the standard ASME B36.19M stainless steel pipe͒, following logic relationship ͑39͒. Finally, the pressure variables, as well as the flow rates and all other variables, are computed again performing another optimization process in GAMS where all the diameters are fixed.

· ͑39͒

4.2.8 Second Optimization Case. Another way to proceed for considering only commercial sizes for pipe diameters is to code each diameter by an integer variable taking the values 1-19, each integer value being associated with a commercial size ͑for example, the value 1 represents the size 0.154 m, while the value 2 represents the size 0.178 m͒. Yet, in this case, the problem combinatorics explodes because the problem would involve, in addition to the 19 binary variables associated with the directions of As mentioned, the pressure ratio is comprised between a lower bound that is obviously equal to one and its upper bound considered equal to two. This value is a classical one adopted in the dedicated literature. The horsepower of a compressor station can be calculated starting from its pressure ratio and its throughput. Gas consumption rate in a station is obtained using Eq. ͑20͒. The gas consumption at each station depends linearly on its power consumption from Eqs. ͑17͒ and ͑20͒.

It has been imposed that if a compressor station must be considered on a line, it must work with a power greater than a lower value assumed here equal to P -=1000 kW. This constraint is taken into account by declaring a special variable type called semicontinuous in the GAMS input file. Semicontinuous variables are those whose values, if nonzero, must be above a given minimum level. This can be expressed algebraically as the logic relation ͓40͔: arcs, 19 additional integer variables, each one taking 19 values. The total number of possibilities is 2 19 ϫ 19 19 = 1.03 30 . So, in this case, the MINLP procedure either diverges or does not converge after several tens of CPU time.

4.2.9 Results. Because of the nonconvergence of case 2, only the results of case 1 are reported here. As already mentioned, after the post-treatment of the diameters, the pressure variables, as well as the flow rates and all the other variables, are computed again by solving another optimization problem in which the diameters are fixed, and the optimizer gives a total cost, discounted over 10 years, of 4.024ϫ 10 6 euros. The computed diameters are shown in Table 8. The optimum pressures subjected to the nodes are shown in Table 9. The average compressibility factor in each pipe can be obtained for the optimum solutions by applying the results presented in this table, using Eqs. ͑5͒ and ͑8͒.

The mass flow rate of delivery gas for the demand nodes is fixed as a constraint shown in Table 4 but the optimum flow rate of injected gas to the network at the supply nodes or the storages is calculated by the algorithm. They are shown in Table 10. The optimum mass flow rates across the different arcs are reported in Table 11. Using the information given in Tables 8,9, and 11, optimum average gas velocity and erosional velocity for each arc can be calculated simply, following the equation linking the mass flow rate, the density, and the velocity of gas ͑m = ͑ / 4͒D 2 v͒ and Eq. ͑16͒, respectively. For the pipelines on which a compressor station exists, the values presented in Table 11 express the gas flow rate before the station. So, there are compressor stations on the pipelines 10-11, 11-12, and 11-17. Their operating conditions are shown in Table 12.

The compressor stations with low-pressure ratio cannot be eliminated because their throughput is high enough to increase the station power consumption significantly.

The discharge pressure for each station can be easily obtained using Eq. ͑19͒. For each compressor station, the sum of the throughput and the fuel consumption is equal to its gas flow rate for each pipeline presented in Table 11. Obviously the mass flow rate of the gas through a pipeline after its compressor station is equal to the station throughput. The objective of this work was to develop a general methodology for gas transmission pipeline modeling based on an optimization-oriented framework, able to embed formulations from design to operational purposes. For this purpose, steady-state behavior of the gas is considered and assumed in the momentum and mass balances.

Although various optimization techniques can be used, the choice of a deterministic one has guided the solution strategy, since it is generally recognized that this kind of methods is particularly well-fitted to take into account the important number of constraints that are likely to be involved in the problem formulation. An adequate solver, namely, SBB, within GAMS environment was selected since this optimization toolbox is often considered as a standard in process systems engineering.

A strategy to solve the optimal design of gas transmission network of a given topology is proposed. The approach is illustrated by an example taken from literature, inspired from the Belgian network. It comprises 20 nodes linked together with 20 arcs. The considered objective function is the total annualized cost, including the investment and operating costs. The optimization strategy provides the main design parameters of the pipelines ͑diameters, pressures, and flow rates͒ and the characteristics of compressor stations ͑location, suction pressure, pressure ratio, station throughput, fuel consumption, and station power consumption͒ to satisfy customer requirements.

Fig. 2 Fig. 3

 23 Fig. 1 Schematic view of the different parts of a natural gas delivery system

  Fig. 4 Stress in pipeline subjected to internal pressure due to gas flow Qs (m 3 /s)

  compressors, types of compressor, and the cycling of each compressor. A reduction technique for natural gas transmission network optimization problems was implemented by Rios-Mercado et al. ͓12͔. Martinez-Romero et al. ͓13͔ used the software package "GAS NET." A MINLP model for the problem of minimizing the fuel consumption in a pipeline network was implemented by Cobos-Zaleta and Rios-Mercado ͓14͔. Mora and Ulieru ͓15͔ determined the pipeline operation configurations requiring the minimum amount of energy ͑e.g., fuel and power͒ needed to operate the equipment at compressor stations for given transportation requirements. Chauvelier-Alario et al. ͓16͔ developed CARPATHE,a simulation package GdF Suez, French Company for representing the behavior of multipressure networks and including functionalities for both network design and network operation. Optimization methods for planning reinforcement on gas transportation networks and for minimizing the investment cost of an existing gas transmission network were used byAndré et al. ͓17͔.

Table 2 Supply nodes' characteristics

 2

				Relative	Relative
		Minimum	Maximum	minimum	maximum
	Supply	value MMSCMD	value MMSCMD	pressure	pressure
	node	͑kg/s͒	͑kg/s͒	͑bar͒	͑bar͒
	1	08.870 ͑087.918͒ 11.594 ͑114.917͒	00.0	77.0
	8	20.344 ͑201.646͒ 22.012 ͑218.178͒	50.0	66.2
	Total	29.214 ͑289.564͒ 33.606 ͑333.095͒
	Note: Relative pressure means absolute pressure added to atmospheric pressure,
	MMSCMD is the abbreviation of million standard cubic meters per day, and standard
	condition is taken at 25°C and 101,325 Pa.

Table 3 Characteristics related to gas storage nodes

 3

			Relative	Relative
		Maximum	minimum	maximum
	Storage	outflow MMSCMD	pressure	pressure
	node	͑kg/s͒	͑bar͒	͑bar͒
	2	08.40 ͑083.259͒	0.0	77.0
	5	04.80 ͑047.577͒	0.0	77.0
	13	01.20 ͑011.894͒	0.0	66.2
	14	00.96 ͑009.515͒	0.0	66.2
	Total	15.36 ͑152.245͒

Table 4 Daily gas demands and their pressure restrictions at delivery nodes

 4

		Gas demand	Relative	Relative
	Delivery	MMSCMD	minimum pressure	maximum pressure
	node	͑kg/s͒	͑bar͒	͑bar͒
	3	03.918 ͑038.834͒	30.0	80.0
	6	04.034 ͑039.984͒	30.0	80.0
	7	05.256 ͑052.096͒	30.0	80.0
	10	06.365 ͑063.089͒	30.0	66.2
	12	02.120 ͑021.013͒	00.0	66.2
	15	06.848 ͑067.876͒	00.0	66.2
	16	15.616 ͑154.783͒	50.0	66.2
	19	00.222 ͑002.200͒	00.0	66.2
	20	01.919 ͑019.021͒	25.0	66.2
	Total	46.298 ͑458.896͒

Table 1 Properties of the flowing gas through the network

 1

	Gas component	CH 4	C 2 H 6	C 3 H 8
	Mole percent in gas	70	25	05
	Molecular mass ͑g/mol͒	16.04	30.07	44.10
	Lower heating value at 15°C and 1 bar ͑MJ/ m 3 ͒	37.706	66.067	93.936
	Critical pressure ͑bar͒	46.0	48.8	42.5
	Critical temperature ͑K͒	190.6	305.4	369.8
	Heat capacity at constant pressure ͑J/mol K͒	35.6635	52.848	74.916

Table 5 Upper limits on the pressures of interconnection nodes

 5 ͑32͒ where unit operating cost C o is assumed equal to 8.2 € / ͑kW year͒ ͓57͔ as follows:

		Relative		Relative
		maximum		maximum
	Interconnection	pressure	Interconnection	pressure
	node	͑bar͒	node	͑bar͒
	4	80.0	17	66.2
	9	66.2	18	63.0
	11	66.2

Table 7 Commercial pipe sizes

 7

	Pipeline internal diameters in meter according to standard schedule
			No. 40
	0.154	0.178	0.203	0.227	0.254
	0.279	0.305	0.336	0.387	0.438
	0.489	0.540	0.590	0.641	0.692
	0.743	0.794	0.844	0.895

Table 8 Optimum pipeline diameters

 8 : The wall thickness of the pipelines can be calculated following Eq. ͑33͒.

		Case 1, diameter
	Arc	͑m͒
	1-2	0.489
	2-3	0.743
	3-4	0.692
	5-6	0.305
	6-7	0.154
	7-4	0.305
	4-14	0.590
	8-9	0.844
	9-10	0.844
	10-11	0.692
	11-12	0.641
	12-13	0.590
	13-14	0.590
	14-15	0.844
	15-16	0.743
	11-17	0.305
	17-18	0.279
	18-19	0.336
	19-20	0.305

Note

Table 9 Pressure subjected to the nodes

 9

		Case 1, relative pressure
	Node	͑bar͒
	1	74.54
	2	72.56
	3	71.52
	4	67.11
	5	77.00
	6	31.33
	7	30.51
	8	65.60
	9	65.07
	10	62.91
	11	60.58
	12	64.52
	13	56.76
	14	55.64
	15	54.03
	16	50.00
	17	66.21
	18	55.61
	19	33.52
	20	30.92

Table 10 Mass flow rate of delivery gas or supply gas at the concerned nodes

 10

		Case 1, relative pressure
	Node	͑bar͒
	1	114.92
	2	83.26
	5	40.93
	8	201.65
	13	10.02
	14	9.52

Table 11 Gas mass flow rate through pipelines

 11

		Case 1, flow rate
	Arc		͑kg/s͒
	1-2		114.85
	2-3		198.04
	3-4		159.13
	5-6		40.86
	6-7		0.81
	7-4		51.36
	4-14		107.70
	8-9		201.58
	9-10		201.51
	10-11		138.36
	11-12		116.83
	12-13		95.70
	13-14		105.65
	14-15		222.80
	15-16		154.85
	11-17		21.51
	17-18		21.43
	18-19		21.36
	19-20		19.09
	Table 12 Optimum conditions related to necessary compres-
	sor stations
	Pipeline that requires a compressor station	10-11	11-12	11-17
	Suction pressure ͑bar͒	63.39	61.42	59.04
	Pressure ratio	1.022	1.147	1.142
	Station throughput ͑kg/s͒	138.34 116.73 21.49
	Fuel consumption ͑g/s͒	18.5	101.8	18.5
	Station power consumption ͑kW͒	1000	5520	1000
	Position from upstream node ͑km͒	3.56	1.00	9.50

Corresponding author. Contributed by the Petroleum Division of ASME for publication in the JOURNAL OF ENERGY RESOURCES TECHNOLOGY. Assoc. Editor: Andrew K. Wojtanowicz.


M Boulos

S Guillemet-Fritsch

Q Nguyen

Z Valdez-Nava

J Farenc

B Durand

Lanthanum Doped Barium Titanate Materials with Optimized Properties for High Capacity Materials Application

Keywords: Powders-chemical preparation, Ba 1 .,La,. TI 1 .x1 4 0 3, Sintering, Dielectric properties, Capacitors

A series of lanthanum doped barium titanate, of formula Ba 1 .,La,. TI 1 .x1 4 0 3 (0:,; x s 0.138), were prepared an d characterized. A soft chemistry technique, the coprecipitation method, was used to synthesize the powders. Fin;t, the oxalate precursors were obtained, starting from an aqueous solution containing TIO 3 , BaO 2 .2H 2 0 and LaO 3 .7H 2 0, coprecipitated with an alcoholic solution of oxalic acid. After calcination at 850 °C, the oxide powders of composition Ba 1 .,La,. TI 1 .x1 4 O 3 were obtained. The oxide powders crystallize in the cubic perovskite phase. SEM-FEG micrographs show that the powders are highly agglomerated. The size of the microcrystallites was estimated to 30 n m . Ceramics have been prepared by conventional sintering treatment (10 or 20 h at l250 °C or 1275 °C). The ceramics are well densified without any abnormal grain growth. The Curie temperature (T e l decreases linearly as the lanthanum content increases and, at the same time, an increase in the relative permittivity (E,) maximum at Tc occurs and reaches 10 000 in the most favourable case. This optimal composition is retained for possible applications.

. Introduction:

Barium titanate perovskite (BaTiO 3 ) is a ferroelectric and piezoelectric material that shows high dielectric perrnittivity and high resistivity. These characteristics allow such appli cations as multilayer capacitors (MLC) or positive tempera ture coefficient thermistors (PIC) ... [START_REF] Huybrechts | The positive temperature coefficient of resistivity in barium titanate[END_REF][START_REF] Choi | Characteristics of BaTiO 3 Particles prepared by spray coprecipitation method using titanium acylate based precursors[END_REF][START_REF] Hu | Homogeneous ( co )precipitation of inorganic salts for synthesis of monodispersed barium titanate particles[END_REF]. When pure, BaTiO 3 is highly resistive at room temperature, but its elec trical resistivity can be lowered by some dopants, like rare earth elements. Doping ofBaTiOrbased ceramics is of great importance in the fabrication of electric and electronic de vices. Owing to the intrinsic capability of the perovskite structure to host ions of different sizes, a large number of various dopants can be accommodated in the BaTiO 3 lattice. However, the effect of the specific impurity on the electrical properties depends on the substitution site (a trivalent ion be haves as an acceptor when substitution occurs at the Ti site or as a donor when it substitutes at the Ba site) as well as on the nature of the defect/electronic charge compensating de fect. For commercial applications of BaTiO 3 , it is necessary to shift the ferroelectric Curie temperature, Tc, towards room temperature and at the same time to broaden the maximum of permittivity, e' m a x , so that a higher permittivity can be obtained over a wider range of temperature. This effect is usually obtained by double-doping, in which for example some of the Ba is replaced by Ca and at the same time, some of the Ti is replaced by Zr. Such double-doped materials have been used for many years as capacitor materials.

Aliovalent, A-site dopants such as La 3+ are known to in duce n-type semiconductivity in BaTiO 3 when added at rel atively low concentration [ 4-7). The origin of the electrical properties of the so-called donor-doped BaTiO 3 continues to be a matter of controversy. Jonker et al. [4] have established that with increasing donor content, i.e. partially replacing Ba by a trivalent rare-earth ion, the room temperature resistivity goes through a deep minimum.

West et al. [8-1 O

] recently reported that La-doped BaTiO 3 retains its insulating character when processed so as to avoid oxygen Joss and that Tc and E' max values depend upon La content. They showed that pure BaTiO 3 cerarnics are insu lating for ~ 0 3 at. % substitution of Ba by La in ceramics fired at temperatures higher than 1350 °C. In this material, the resistivity passes through a minimum value of~ 10 Q.cm and then increases again to insulating levels with further sub stitution of La. This behaviour is further complicated by loss of oxygen in ceramics heated above l 300 °C and subsequent partial reoxidation on cooling in air [4,[START_REF] Saburi | Properties of semi-conductive barium titanate[END_REF]. Thus, the electrical microstructure of ceramics strongly depends on the sintering and post sintering conditions. For instance, it is possible to prepare materials presenting semi-conducting grain cores with insulating grain boundaries [6].

Most of the La-doped BaTiO 3 sarnples reported in the lit erature have been obtained by solid/solid reactions. Only Costa and Mantas [ 11] used the oxalate method to produce the sole composition B"<:)_ 997 Lao_ 00 3 Ti 1 , 0045 0 3 . Their pow ders had small grain size and were disaggregated. However, the sintering process was difficult to optimize in this material due to the presence of large pores uniformly distributed in the matrix and the resulting grain size was not uniform. They showed also that the transition temperature decreases from 120 °C to l 15 °C. Recently, Beltran et al. [12] synthesized pure and La doped-BaTiO 3 by the sol-gel route with x = 0, 0.003, 0.03 and 0.10.

The purpose of this paper is to present new results con ceming the elaboration of La-doped BaTiO 3 powders using a very simple way, the coprecipitation method i.e. an alter native route to the solid/solid one allowing a better control and a sharp distribution of the grains size. The physicochem ical and the electrical properties of the sintered ceramics are then described. The properties of the materials presented here are at least equal or higher to the ones previously re ported.

Experimental

Sample preparation

The La-doped BaTiO 3 powder was synthesized by the oxalate method followed by a calcination treatment. Sarnples of composition Ba 1 .xLaxTi J .x/ 4 O 3 were prepared with dif ferent x values (0, 0.026, 0.036, 0.048, 0.58, 0.090, and 0.138). The starting materials were BaC1 2 .2H 2 O (Prolabo), TI Cl 3 (Prolabo,d = 1.20, % min= 15) and LaCl 3 .7H 2 O (Pro labo, 99%). They were dissolved in water in various propor tions and the coprecipitation was performed by addition of a solution of oxalic acid dissolved in ethanol. TICI 3 was cho sen because it was easier to handle than TICJ 4 . The full ox idation of Ti 3+ in 1i 4+ in the solution was ensured by air bubbling during the reaction which was evidenced by the Sketch of the co precipitarion method color variation from brown to yellowish. The volume of wa ter was taken much lower than the one of ethanol, in order to decrease the dielectric permittivity of the precipitation me dia. Since the nucleation of the particles is favoured com pared to their growth, smaller particles, homogeneous in size and composition, are likely to be obtained [START_REF] Villette | Sur l'.élaporation par chimie douce de ferrites de cuivre A forts champs coercitifs[END_REF]. The solution was aged for a couple of hours, and the obtained precipitate was centrifuged. The precursors were then pyrolized in air at 850 °C for 4 h to obtain the oxides. The sketch of the re action is given in figure 1. These powders were then mixed with an organic binder, then uniaxially pressed into disk (di arneter: 6mm and thickness: 2 mm) under 250 MPa pressure. The sintering temperature and the dwell limes have been op timized in order to get a high densification value and at the sarne time to avoid the formation of large grains. Hence green cerarnics were sintered in static air at l 250 °C for 10 and 20 h or at 1275 °C for 10 h.

Characterization

The chemical composition of the oxide powder was de termined accurately by induced coupled plasma spectrosco PY (ICP AES Thermo-Optec ARL 3580). Electron microscopy JEOL JSM 6700F SEM-FEG and scanning elec tron transmission JEOL 2010 were used to observe the oxide powder morphology. The specific surface area was deter mined thanks to a Micrometrics accusorb 2100 E.

The structure was determined by X-ray diffraction analy sis: a SEIFERT XRD-3003-TT diffractometer using the Cu Ku radiation (À.= 0.15418 nm) allowed to collect the data.

For electrical measurements, silver electrodes were screenprinted on the opposite surfaces of the pellets and heat treated at 700 °C for 30 min. The complex effective permit tivity € * was measured thanks to a Hewlett Packard HP 4284A precision LCR meter for frequency ranging between I 00 Hz to l MHz and for temperature between -150 to l 40 °C. The voltage magnitude was I V. The chemical composition, the structure and the specific surface area of the different powders prepared by pyrolysis of oxalate precursors at 850 °C for 4 h are reported in Table 1. The X-ray diffraction patterns of "undoped" Ba1iO 3 and La-doped-Ba1iO 3 are given in Fig. 2 (a andb). If undoped BaTiO 3 prepared by coprecipitation crystallises in the te-tragonal form, the cubic phase is observed when the powders are doped with La 3+ (x è!: 0. 026). The lattice parameters of the different powders were cal culated using the Rietveld method (Table 2). Changes in the lattice parameters at room temperature as a function of x are shown in Fig. 3 . For0 s x s 0. 026, a decrease in tetragonality is observed when increasing x up to 0. 026. The structure of B3cJ.97Lllo.026TJ o.99 3 O 3 is cubic with a= 4. 0037 (4)À. still increasing x (x > 0. 026) the cubic symmetry is retained, and the lattice parameter decreases down to a = 3 . 994 l A for ~ x = 0. 090. For x è!: 0. 090, increasing La content has little ef fect on the unit cell which remains cubic and, for x = 0. 13 8, a = 3 . 9940 A. The same tendency, even if they didn 't follow the same route, has been obtained by West et al. [START_REF] Morrison | Electrical and structural characteristics of lanthanumdoped barium titanate ceramics[END_REF], and can be explained by a cation size effect. The smaller La 3+ cation stabilizes the cubic form as predicted by Gold schmidt's tolerance factor for undistorted and distorted per ovskites [ 15 ]. The reasons why the unit cell parameter, a (and hence unit cell volume) should appear to reach a lower limit of a -3 .995 À with increasing x are still unclear. The li mit may signify the highest possible packing density for a (Ba,La)TiO 3 cubic close packed lattice [ 16 J. Measurements of the specific surface area show that the particle size de creases when increasing the amount of lanthanum (Table 1). It is well known that the addition of donor dopants above the critical concentration (0 2-0 5 al % ) which corresponds 10 the semi conductor 10 insula1ion transition, greatly inhibits grain growth [START_REF] Chan | Compensating defects in highly donor-doped BaTiOJ °[END_REF](18)[START_REF] Xue | The infl uence of ionic radii on the incorporation of trivalent dopants into BaTiO 3[END_REF][START_REF] Belous | Microstructure and semiconducting properties of Barium titanate containing heterovalent substituents on titanium site[END_REF][START_REF] Bussem | Effects of grain growth on the distri bution of BaTiO 3[END_REF][START_REF] Kurata | Semiconducting-lnsulating Transition for highly-Donor Doped Barium titanate ceramics[END_REF]. The variation of the specifi c surface area as a function of x is shown in Fig. 4. A quasi-linear evolution is observed up to x = 0.090. SEM-FEG observations of calcined powders (Fig. 5) con finn that the La doped-BaTiO 3 powders have smaller particle size (< 100 nm) than the undoped ones. For undoped barium titanate (Fig. Sa), the particle size is -250 nm. When doped, the powders become highly agglomerated and partial sinter ing has occurred. TEM observations (Fig. 6) on doped BaTiO 3 powders (x = 0.138) showed the agglomerated state of the powders.

Ceramics

The sintering treatment has been performed in air at 1250 °C for lOh and 20h or at 1275 °C for 10 h. The densities of the sintered ceramics are shown in Table 3. In ail cases, the oeramic's densification is over 89%. An extension of the dwell time of sintering from lO to 20 h increases the density whereas increasing the sintering temperature from 1250 to 1275 °C seems to lead to the opposite effect.

The structure of the ceramics is similar to the one of the oxides of same composition. Tetragonal perovskite was an alysed by X-ray diffraction for undoped BaTiO 3 ceramics and cubic perovskite was identified for La-doped-BaTiO3.

The average grain size determined from SEM images is 1 µm for pure BaTiO 3 ceramics. Fig. 7 shows the microstructure of various ceramics con taining different arnounts of La 3+ and sintered at l 250 °C for 20 hours. The compositions comprised between 0 sx sû.058 of a -3.995 A with increasing x are still unclear. The limit may signify the highest possible packing density for a (Ba,La)TIO 3 cubic close packed lattice [START_REF] Buscaglia | Influence of foreign ions on the crystal structure of BaTiO 3[END_REF] . Measurements of the specific surface area show that the particle size de creases when increasing the amount of lanthanum (Table I ). It is well known that the addition of donor dopants above the critical concentration (0.2-0.5 at %) which corresponds to the semi conductor to insulation transition, greatly inhibits grain growth [START_REF] Chan | Compensating defects in highly donor-doped BaTiOJ °[END_REF][18][START_REF] Xue | The infl uence of ionic radii on the incorporation of trivalent dopants into BaTiO 3[END_REF][START_REF] Belous | Microstructure and semiconducting properties of Barium titanate containing heterovalent substituents on titanium site[END_REF][START_REF] Bussem | Effects of grain growth on the distri bution of BaTiO 3[END_REF][START_REF] Kurata | Semiconducting-lnsulating Transition for highly-Donor Doped Barium titanate ceramics[END_REF] . The variation of the specific surface area as a function of x is shown in Fig. 4. A quasi-linear evolution is observed up to x = 0 .090. SEM-FEG observations of calcined powders (Fig. 5) con finn that the La doped-BaTIO 3 powders have smaller particle size (< 100 run) than the undoped ones. For undoped barium titanate (Fig. Sa), the particle size is -250 nm. When doped, the powders become highly agglomerated and partial sinter ing has occurred. TEM observations (Fig. 6) on doped BaTiO 3 powders ( x = 0.138) showed the agglomerated state of the powders.

Ceramics

The sintering treatment has been perfonned in air at 1250 °C for !Oh and 20h or at 1275 °C for 10 h. The densities of the sintered ceramics are shown in Table 3. In ail cases, the ceramic's densifi cation is over 89%. An extension of the dwell time of sintering from 10 to 20 h increases the density whereas increasing the sintering temperature from 1250 to 1275 °C seems to lead to the opposite effect.

The structure of the ceramics is similar to the one of the oxides of same composition. Tetragonal perovskite was an alysed by X-ray diffraction for undoped BaliO3 ceramics and cubic perovskite was identifi ed for La-doped-BaliO3.

The average grain size detennined from SEM images is 1 µm for pure BaliO3 ceramics. Fig. 7 shows the microstructure of various ceramics con taining different amounts of La 3+ and sintered at l 250 °C for 20 hours. The compositions comprised between 0 :sx :SÜ.058

Electrical properties

A summary of the changes in the relative permittivity, meas ured at l kHz, venms the temperature (-l20 °C :s; T :s; 160 °C), is given in Fig. 8 for La--Ooped-Ba1iû:3 ceramics sintered at l250 °C/20h. The permittivity behaviour of La-doped composi tions, results in an initial increase of the values of E r , followed by a broadening and a decrease of Er, It is important to note that the values of the reaI part of the permittivity are at least one third of the values obtained following the solid/solid route. On the other hand, the los.ses remain lower than l % whatever the samples (this information being rarely given in the bibliography). Such results may appear surprising since it is often con sidered that decreasing the grain size (at least down to 0.8 µm) leads to an increase of the permittivity. As mentioned before, the grains sizes obtained in the ceramics under study are low er than the ones obtained following other routes. Therefore, we may consider that the main driving phenomena able to explain this behaviour would probably not be the same.

The undoped barium titanate ceramics present a classic behaviour with a Curie temperature, Tc -l20 °C, and a shal low maximum at a temperature T oit ~ 25 °C associated with the orthorhombicltetragonal (oit) transition. When doping with La, the Er -temperature profiles for compositions 0 s x s 0.058 retained the essentiaJ signature ofundoped BaTiO 3 , but were displaced to lower temperatures with increasing x. Er, reached a maximum of 10000 at I kHz for x = 0.058. Moreover, for this value, Er associated with the oit transition was barely discernible as a shoulder on the low temperature side of the curve. For x > 0.058, Er, decreased and became increasingly broad with increasing x. The T oit temperature aJso decreased with increasing x, but in a Jess dramatic rate than Tc. For x:.:: 0.09, it was not possible to detect the Curie temperature anymore.

Variation of Tc and T011 as afunction ofxfor doped bariunt titanate

Changes in Tc and To it as a function of x (Fig. 9) indicate a rapid decrease in both phase transition when x increases. Tc decreased approximately linearly at a rate of ~ 22 °C per at. % La 3+ substitution which is in good agreement to the 21 °C per at.% obtained in the case of Ce 3+ [23,24]. The La doping aJso has the greatest reported effect on T o1t , even if it decreased Jess rapidly.

In summary, La-doping "pinches" Tc and T oit together in a similar manner as Zr doping [26]. The main difference is that La-doping decreases both Tc and Toi t • whereas Zr doping decreases Tc and T 01t until they coaJesce. In the present case the phase transition temperatures coaJesce for a dopant con centration of~ 9 at. % La since for this vaJue, the phase tran sition temperatures are no more detectable. The temperature at which Tc and T 011 coaJesce is --130 °C.

The impact of the dwell time has been studied. Fig. lO shows the variation of Er for ceramics (x = 0.058) sintered at l 250 °C for lO and 20 hours. The permittivity vaJues present a slight increase which may be explained by an in crease of the density when increasing the dwell time of sin tering (see Table 3). In order to identify the most Iikelihood phenomena able to explain the observed behaviour, impedance spectroscopy is used. This method is aJI the more interesting since it allows giving a simple and accurate representation of the material behaviour. Hence, the measured complex permittivity (Er): may be transformed in complex impedance Z* using: Z*(w)= 1/ jOJCoE*(w) = Z'+ j Z" Where w is the angular frequency, Co the geometrical ca pacitance (=foA/t, with A the electrode surface and t the ma terial thickness). It is sometimes interesting to plot Z" vs Z', Cole-Cole plot, or Z" ( t) T= c te (or Z" (T) f= c te ), "De bye like" plots.

In this last case the behaviour may be depicted by: 'Z:' = (V (Cw (!+ (tgô) 2 )

Where tgô represents the dielectric losses defined as tgÔ= i::"IE'= 1/RCw in the parallel RC mode!. It is interesting to note that when Wmax,:: = 2nfmax ,:: = WmaxR.C = 1, the value The results for the samples with x= 0.058 are presented in Fig. 11 and 12. For a measuring temperature of 250 °C, two semicircles are observed in the Cole -Cole plot. The impedance of the equivalent electrical circuit may therefore be written as:

R1 Rz Z ( w) -l + jw-c 1 + 1 + jw. 2
The high frequency semicircle is ascribed to the material grain bulk whereas the low frequency semicircle corresponds to the response of the grain boundary. In Cole-Cole plots, it is also interesting to calculate the two angles 81 and 82 by which semicircles are depressed below the real axis. These angles are related to the width of the relaxation time distri bution.

The calculated values of the O angles are very small: 81 =0.034 rad and 82 = 0.0495 rad respectively for the bulk and for the grain boundary. Narrow distribution function of time relaxations are therefore observed confirrning the homoge neity of the ceramics structure as already observed by SEM. Since, the mean relaxation time Tm is equal to l .0467x 10-5 s for the bulk and of 3.8007 X w-3 s for the grain boundary it is easy from the maximum value of Z" in the Z"(T) plot (Fig. 14) to deterrnine the value of R. Knowing 'tm and R, the electrical characteristics of each elementary cell may then be deterrnined. They are equal to:

-271.17 kQ (for R 1 ) and 38.6 pF (for Cj)(bulk) -3.13 MQ (for R 2 ) et 1.23 nF (for C 2 ) (grain boundary) These values are different to the ones reported up to now. Hence the value� of R and C deterrnined for the bulk may explain the large value obtained for Er.

The same calculations may be undertaken for the samples sintered at the same temperature but for a longer dwell time. The results reported in Table 4 prove that increasing the dwell time has a little, if not any, effect on the different char acteristics issued from impedance spectroscopy analysis.

Oxygen Joss has been avoided by heating in air at lower temperature eg l 200 °C where these losses are negligible. j

Conclusion

The oxalate method used in this work is a very simple way to prepare La doped BaTIO 3 materials. The obtained ox ide powders had particle size lower then 100 nm and were agglomerated. The effect of La-doping on the pennittivity temperature profile showed that the replacement of Ba 2+ by smaller La 3+ cation on the A-site and the presence of B-site vacancies results in a rapid decrease in T c and T o i t • Doped samples showed maximum permittivity values higher than 10000 for a value of x=0.058. These values are mainly re lated to the quality and homogeneity of the ceramics. La cat ions seem a good choice of dopants for improving barium titanate ceramics properties for specific applications.

Figure 1 :

 1 Figure 1:

Figure 2

 2 Figure 2:

Figure 5 :

 5 Figure 5:

  Figure 5:

Figure 6 :

 6 Figure 6: TEM micrographs of doped barywn. titanate powders of composition x = 0.138

Figure 7 :

 7 Figure 7: SEM micrographs of doped bariWll ritanate Ba 1 . .,La,1i 1 _,14O 3 ceramics sintered al 1250°C for 20 hours x = 0.048 (a), x = 0.058 (b), x = 0.09 (c) and X = 0.J 38.

  Figure 8: Relative pennittivity as function of tempera/ure of doped bariWll titanate ceramics Ba 1 . .,la,1i 1 .v 4 O 3 with various compositions, 0 s x s 0.138

  Figure 9:

150

  Figure JO:

Z

  Figure Il: Cole-Cole plot measured at different temperatures for a doped barium ti tanate Ba 1 _,La_ ,Ti 1 _,1 4 O 3 ceramics x=0.058 Z"(ohm) 1.0x10 7 ....... �.....-""'..-� ....... ....--.-� ...... ���..--......,....,..-.-,

Table / :

 / Structure. chemica/ composition and specific su,face area of differenl doped barium titanate powders Ba 1 .,La,1i 1 .tJ 4 O 3 .

						Specific
	Rcfcrence	X	Structure	Chcmical composition •	surface arca
						(m 2 /g)
		DopO	0.00	tctragonal		B a T i0 3	6.00
		Dop0 3	0.026	cubic	Bao.974Lao. 026 Tfo.99303	7.55
		Do p04	0.036	Cubic	Bao.964Lao. o36 T'0.99103	8.16
		Dop05 Dop06 Dop09	O.û48 0.058 0.090	Cubic Cubic Cubic	Bao.9s1Lao04 s T10.9s703 Bao.942Lao. oss Tio.98s03 Bao.90sLao. 090 Tio.91103	9.10 9.50 11.0
		Dopl4	0.138	Cubic	Bao.862L"o.138 T'0.96503	11.6
		4.024
	� 4.020 I!! 4.016
		4.012
	4.008 3 4.004 a.
	�	4.000 3.996
		3.992
		0.00	0.03	0.06	0.09	0.12
					X
					Figure 3:

*Formula ca/culated after normali:ing Ba +, La = 1.

0.15

Variation in lattice parameters at room tempera/ure as fonction of x of doped barium titanate powders Ba 1 . .Lax1i 1 _, 14 03.

Table 2 :

 2 Structure and lattice parameters of different /anthanum doped barium titanate powders Ba 1 jA_ ,1i 1 .v�0 3 .

					Lanice
	Structure	a(À)	c(À)	c/a	volume
					(ÀJ)
	Tctr•gonal	4.0006	4.0215	1.005	64.364
	Cubic	4.0037	4.0037	1	64.137
	Cubic	4.0028	4.0028	1	64.133
	Cubic	4.0020	4.0020	1	64.095
	Cubic	4 0009	4.0009	1	64.045
	Cubic	3.9941	3.9941	1	63.717
	Cubic	3.9940	3.9940	1	63.712
	0

.00 0.04 0.08 X 0.12 0.16 Figure

4: Variation ofswface specific area ar room temperature as function of x ( L15 .., = 3%)

Table 2 :

 2 Stnu:ture and latrice parameters of diff erent lanthanum doped barùun titannte powders Ba1.J.,11,Ti 1 . .,NO1.

					Latticc
	Structure	a(ÂJ	C(Â)	c/a	volume
					(À3)
	Tctragonal	4.0006	4.0215	1.005	64.364
	Cubic	4.0037	4.0037	1	64.137
	Cubic	4.0028	4.0028	1	64.133
	Cubic	4.0020	4.0020	1	64.095
	Cubic	4.0009	4.0009	1	64.045
	Cubic	3.9941	3.9941	1	63.717
	Cubic	3.9940	3.9940	1	63.712
	0

.00 0.04 0.08 X 0.12 0.16

Fil/ure 4: Variation of s11rface specific area at room temperature asfunction ofx (Af_. = 3%)

Table 3 :

 3 Density of doped barium titanate Ba 1 _,La,1i 1 .,14 0 3 ceramics, as a function sintering conditions.

	X	Sintering T Sintcring dwell ('C) timc(h)	Density (g/Cc)	Densification (%)
	0.00	1250	10	5.33 :t 0.()1	89
		1250	20	5.64 :t 0.02	94
	0.026	1250	20	5.51 ± 0.07	91
	0.036	1250	20	5.75 :t 0.07	96
		1250	IO	5.52 :t O.ü7	92
	O.�	1250	20	5.58 ± O.IO	93
		1275	10	5.56 ± 0.10 :.	93
		1250	10	5.57 ± 0.02	93
	0.058	1250	20	5.78 ± 0.07	96
		1275	10	5.35 ± 0.06	90
		1250	IO	5.32±0.11	90
	0.090	1250	20	5.70 ± 0.04	95
		1275	10	5.62 ± 0.12	94
		1250	10	5.63 ± 0.08	94
	0.138	1250	20	5.81 ± 0.08	97
		1275	IO	5.52 ± 0.01	92

Table 4 :

 4 Summary of the main electrical characteristics of the material Ba 1 _.J-a,Ti 1 ., 14 O 3 for two sintering dwell times and x=0JJ58

		X= 0.058	X= 0.058
		Ts= l 250'C, t= !Oh	Ts= l250 °C, t= 20h
	8 1 [rad]	0.034	0.030
	82 [rad]	0.0495	0.050
	Mean relaxation	1.0467 xto• 5	1.416 x10• 5
	timeTmJ [s]
	Mean relaxation	3.8007 x10• 3	4.0807 xto• 3
	timeTm2 [s]
	Resistance R 1	0.271	0.441
	[MQJ
	Resistance R2	3.13	3.47
	[MQJ
	Capacitance C 1	38.6	32.13
	[pF]
	Capacitance C2	1230	1176
	[pFJ
	Equivalent circuit
		-èJc:J-1\	R,	1\

Acknowledgrnents

The authors thank Ch. Calmet for realizing the FEG-SEM pictures, A. Barnabé for lattice parameters measurements and JJ. Demai for TEM pictures.