Didier THEILLIOL Philippe Weber
email: philippe.weber@cran.uhp-nancy.fr

Christophe Aubrun
email: christophe.aubrun@cran.uhp-nancy.fr

Component Reliability in Fault Diagnosis Decision-Making based on Dynamic Bayesian Networks

Keywords: Model-based fault diagnosis, Bayesian Networks, Reliability, Markov chains, Decision-making

The decision-making in fault diagnosis methods generally relies on the analysis of fault signature vectors. This paper presents a new approach of decision-making in the case of the signature vectors for various identical or similar faults. The main contribution of the paper consists in the fusion between the reliability and the evaluation of the residuals in order to increase the fault isolation efficiency. The decision-making, formalised as a bayesian network, is established with a priori knowledge on fault signatures, false alarm and missing detection probability, on line component state estimation computed by a bayesian fusion of the component reliability and measurements. The effectiveness and performances of the method are illustrated on a heating water process corrupted by various faults.

INTRODUCTION

Modern control systems are becoming more and more complex. Consequently, there is a growing demand for fault diagnosis to increase the reliability of such systems. In this context, fault diagnosis domain has gained increasing consideration. A fault is considered as malfunction in the actual system that tends to degrade the overall system performances. In this paper, our attention is focused on model-based fault diagnosis which makes use of the analytical relationships between measured variables. A short historical survey on fault diagnosis can be found in [1] and various approaches have been reported in [START_REF] Patton | Issues of fault diagnosis for dynamic systems[END_REF]. The fault diagnosis procedure consists in three stages:

• Residuals generation: it consists in associating a model-observation pair to evaluate difference with respect to the normal operating conditions;

• Residuals evaluation: the residuals are compared to some predefined threshold according to a test and, at this stage, symptoms are produced;

• Decision making: the role of the decision making is to decide according to the symptoms, which elements are faulty, that is isolation. This requires the design of residuals that are close to zero in the fault-free situations while clearly deviating from zero in the presence of faults. These residuals possess the ability to discriminate between all possible modes of faults, which explain the use of the term decision-making. Classically, decision-making is realized according to an elementary logic. Nevertheless, when multiple faults or false alarms occur, the faults cannot be isolated [START_REF] Korbick | Fault Diagnosis, Models, Artificial Intelligence, Applications[END_REF]. Some specific mathematics algorithms can improve the efficiency of the decision-making, for instance: [START_REF] Chen | Robust model-based fault diagnosis for dynamic systems[END_REF], have proposed methods which increase the robustness of residual by decoupling the effects of faults from each other and from the modelling errors and uncertainties.

• [4],
• [START_REF] Sauter | Fault diagnosis in systems using fuzzy logic[END_REF] have developed a method based on adaptive threshold approach to reduce the sensitivity of the residuals evaluation against false alarms. However, in any case, the binary data produced by residuals evaluation are poor in information. Moreover, the degree of isolability based on Hamming distance as suggested [START_REF] Cassar | Approche structurelle de la conception de systèmes de surveillance pour des procédés industriels complexes[END_REF] is very low. Consequently, some other knowledge related to the residuals should be considered for isolation. [START_REF] Theilliol | A hierarchical fault diagnosis method using a decision support system applied to a chemical plant[END_REF] and [START_REF] Evsukoff | Recurrent neuro-fuzzy system for fault detection and isolation in nuclear reactors[END_REF] have combined qualitative and quantitative knowledge to improve the fault diagnosis efficiency. In the spirit of [10], fault isolation performance may be increased by integrating extra information in the diagnosis process. Thus, reliability, classically computed by means of stochastic process model as Markov Chains (MC), defines the a priori behavior of the probability distribution over the nominal states and faulty states of the system. Also, this additional information is seldom used to improve decision-making in model-based fault diagnosis [START_REF] Anrig | Model-based reliability and diagnostic: a common framework for reliability and diagnostics[END_REF].

The aim of the paper is to propose a new approach which increases the performance of the decision-making in fault diagnosis by taking into account a priori knowledge on the system state through a dynamic Bayesian network.

The paper is organized as follows: Section 2 presents the principle of the decision-making in model-based fault diagnosis and fault isolation problem is stated. Section 3 describes the method which achieves the decisionmaking based on the Bayesian network (BN) inference. Then, the proposed approach used to merge both the fault diagnosis and the dynamic reliability is presented. The proposed approach is illustrated through a simulation example in Section 4. Finally, conclusions and perspectives are given in last Section.

PROBLEM STATEMENT

Symptoms generation

The first step produces indicators sensitive to the faults. These indicators are defined as Residuals, noted "r". A residual is defined as the difference between a measurement and the corresponding reference value estimated with a model of the fault-free system (Figure 1). ). The objective is to decouple the faulty effects from each residual. Each residual is designed so that it is sensitive to different faults or subsets of faults.

Usually, the second step of the diagnostic procedure, Residuals evaluation, is based on the assumption that if a fault occurs, the statistical characteristic of the residuals is modified. The residuals evaluation involves statistical testing such as limit checking test, generalized likelihood ratio test, trend analysis test. The output vector of the statistical test, called coherence vector U, can be built according to a test applied to the set of J residuals: U = [u 1 ,…u j ,… u J ] T where u j represents the status of the residuals: u j is equal to "0" when the residual signal is closed to zero in some sense and equal to "1" otherwise. u j is called the symptom associated to the residual r j .

Unfortunately, the residual is corrupted by noise, which affects the decision-making. The efficiency of the detection is related to the false alarm and missing detection probability. In the residual evaluation, the problem may be formulated as a hypotheses testing problem. Let us recall the main definition:

• 0 H : the residual is not affected by a fault; Hits 0

H Reject 0 H ) ), ( ( 0 H k r f ) ), ( ( 1 H k r f ) (k r
( ) true is H H reject p 0 0 = α (1)
and the probability of missing detection correspond to the surface risk β defining the conditional probability:

( ) true is H H reject not p 1 0 = β (2)
More details can be found in the reference book of Basseville and Nikiforov [15] or in the second chapter of the book of Young and Calvert [16].

The risk values α and β are calculated a priori if the density functions ) ), ( ( i H k r f are known. However, this is rarely the case in practice. If the density functions are not available, then an estimation of the probabilities is computed by mean of a frequentist approach based on data including a posteriori detection results.

Incidence matrix

Several approaches have been proposed to generate structured residuals and consequently to generate the incidence matrix [START_REF] Frank | Fault diagnosis in dynamic systems using analytical and knowledge-based redundancy -A survey and some results[END_REF]. Let us consider the following example where three different faults (F 1 , F 2 and F 3 ) can be isolated by designing three symptoms u 1 , u 2 and u 3 (Table 1).

Table 1: Incidence matrix example.

F 1 F 2 F 3 u 1 0 1 0 u 2 1 0 0 u 3 1 0 1
In Table 1, a "1" denotes that a symptom u j is sensitive to a fault (F 1 , F 2 or F 3 ), while a "0" denotes insensitivity to a fault. This table is called an incidence matrix and can be considered as an a priori knowledge. Each column of the incidence matrix represents a fault signature: the vector [0 1 1] T corresponds to the signature of the faulty element F 1 . In this paper, incidence matrix is annotated D with different elements D(n,j), where n is the number of elements suspected to be faulty (n=1...N) and j is the number of residuals (j=1...J).

Fault Isolation

Usually, a very simple logic analysis between each fault signature F n and each coherence vector U is used to isolate the faulty component. In practical cases, false alarms occur and corrupt the decision logic. The coherence vector can be different from all signatures. Therefore, the goal of the decision-making is to minimize the false alarms and missing detection rates due to the effects of modelling errors and unknown disturbances that affect residuals.

Moreover, in spite of the residuals generation and evaluation robustness, a simple logic rule is not efficient enough to isolate faults when simultaneous multiple faults occur [START_REF] Weber | Multiple fault detection and isolation[END_REF]. This is justified by the fact that if D(n,j)=0, then u j cannot bring any information about the occurrence of fault F n , since the residual r j might be different from 0 due to noise or modelling errors or another fault F k (with D(k,j)=1) affecting the system. Thus, new decision-making method is necessary.

Moreover, a new source of information should be integrated in fault diagnosis.

[18] has recently proposed to introduce reliability analysis for fault diagnosis purpose. System reliability analysis allows determining the degradation degree of the system components. The paper aims at developing a method which integrates a dynamic reliability estimation of the system component as presented in the next section. In this new approach the reliability of the components is computed according to the states of the residuals (with the bayesian approach)

and not only an a priori estimation of the reliability. This is the reason why a dynamic computation of the reliability is used.

DECISION-MAKING PROCESS DESIGN

Bayesian network equivalent to incidence matrix a) Bayesian network (BN):

Bayesian network are probabilistic networks based on graph theory.

S

. A probability distribution over these states is defined as a vector

) (n p with ) ( n m s n p = [ ] ( ) M 1,..., m ∈ ∀
is the marginal probability of n being in states n m s . As illustrated in the graph depicted at Figure

n j n i j n s 1 … j n L s i n s 1       = = i j n i n j s n s n p 1 1 … … … … i n M s …       = = i j n M i n L j s n s n p
Concerning the root nodes, i.e. those without parent, the CPT contains only a row including the a priori probability of each state.

Various inference algorithms can be used to compute marginal probabilities for each unobserved node given information on the states of a set of observed nodes. The most classical one relies on the use of a junction tree.

Inference in BN then allows to take into account any state variable observation (an event) such that it updates the probabilities of the other variables. Without observation, the computation is based on a priori probabilities.

When observations are given, this knowledge is integrated into the network and all the probabilities are updated.

Knowledge is formalised as evidence. A hard evidence of the random variable X indicates that the state of the node n is one of the states

{ } n M n n s s ,... : 1 S . For instance, X is in state n s 1 : 1 ) ( 1 = = n s n p and 0 ) ( 1 = = ≠ n m s n p .
Moreover, when this knowledge is uncertain, soft evidence can be used to define the distribution over n.

b) Bayesian network model as an incidence matrix:

The relationship between symptoms and faults are represented by a graph. Obviously a fault can be considered as the cause of the residual deviation. Therefore, some connections can be established from the fault to the symptoms in order to define the relation of causality between fault occurrence and the symptom states.

Whereupon, a Bayesian network can define directly an incidence matrix D(n,j). Let us consider two incidence matrices that represent two different cases of possible isolability conditions (Table 3).

Table 3: Two incidence matrices.

a) b) D(n,j) F 1 F n F N D(n,j) F 1 F n F N u 1 1 0 0 u 1 0 1 1 u j 0 1 0 u j 1 0 1 u J 0 0 1 u J 1 1 0
As depicted in Figure 4, each above incidence matrix is represented as elementary graph with their appropriate arc which corresponds to the link between fault signature and coherence vectors. The fault occurrence probability is modelled as a random variable F n associated to each fault. F n is described by two states {not occurred, occurred}. Moreover, the symptoms are represented also as random variable u j defined over the set of two states: {not detected, detected} with p(u j =detected), if the fault affects the system and the residual r j is detected different from 0.

The probability distribution over the symptom states depends on the false alarms and missing detections. Using Bayesian network model, a CPT is determined to model the relation between variables. In order to compute the probability distribution of symptoms u j , a CPT is defined according to the fault F n parent of u j . For instance, when only one symptom is associated to one fault, as presented at Figure 4a, then the CPT has the structure presented in Table 4:

Table 4: CPT for node u j in Figure 3a.

u j F n not detected detected not occurred n - 1 α n α occurred n β n - 1 β
where the probability n β of missing detection and the probability n α of false alarms both for the fault F n are defined such as:

( ) ( ) occurred F detected not u p occurred not F detected u p n j n n j n = = = = = = β α (3)
Therefore, the probability distribution over the states of the causes (fault occurrence) is based on the residual evaluation result. The node u j is defined as hard evidence. If changes in the residual are detected, then: p(u j = detected) = 1 and p(u j = not detected) = 0 (4) Otherwise:

p(u j = detected) = 0 and p(u j = not detected) = 1 (5)
The Bayes theorem is used to compute ( )

j n u F p
, for instance, relatively to Figure 4a,( )

j n u F p is equal to: ) ( ) ( ) ( ) ( j n j n j n u p F u p F p u F p = (6)
where ( ) Generally, several faults are associated to one symptom, as illustrated in Figure 4b. In this case, the CPT is more difficult to determine and it possesses a different structure from Table 4. For instance, when one symptom is associated to few faults, as presented at Figure 4b where few are equal to two, then the CPT is defined as the Table 5.

Table 5: CPT for node u j in Figure 3b.

u j F 1 F N not detected detected not occurred 1-α j α j not occurred occurred β N 1-β N not occurred β 1 1-β 1 occurred occurred β 1,N 1-β 1,N
The probability that a false alarm exists, when few faults are the cause of the residual deviation, is estimated independently to the faults. On the contrary, the missing detection depends on the fault, as a fault 1 F has not the same impact on the residual as a fault N F . Let us consider that false alarms and missing detections estimated as follows: As presented in Table 5, ( ) The Bayes' theorem is applied in the BN inference algorithm to determine ( )

( ) ( ) ( ) ( ) N N j N N j N j j N j occurred F occurred F not u occurred F not F not u o n F occurred F not u not F not F detected u
H Reject 0 H ) ), ( ( 0 H k r f ) (k r -0.05 0 0.05 0.1 ) (k r risk β 1 ) ), ( ( N N H k r f -0.05 0 0.05 0.1 risk β N ) ), ( ( 1 1 H k r f ) (k r -0.05 0 0.05 0.1 risk β 1,N ) ), ( ( , 1 , 1 N N H k r f ) (k r
N j ,F F u p
j N u u F p , 1
from the states of the symptom u 1 and u j such as:

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ∑ = = i F i J i i i i J i i i i J i J i i J i F u p F u p k n F p F u p F u p k n F p k n u u p k n u u F p k n u u F p 1 1 1 1 1 ) ( ) ( ) ( , , ) ( , , , ) ( , , (10)
In this formulation the probability of failure occurrence is deduced from the reliability of the component ) (k n i and the residuals states ) (k u i and ) (k u J

. Moreover, the back propagation (from ) (k

F i to ) (k n i
) allows to verify the coherence when computing the reliability of the component ) (k n i according to the states of the residuals.

ILLUSTRATION EXAMPLE

Process description and fault diagnosis

To illustrate our approach, a simulation example is considered: a heating water process. The process, presented in Figure 8, is composed of a tank with a section S equipped with two heating resistors R 1 and R 2 . The inputs are the water flow rate Qi, the water temperature Ti and the heater electric power P. The outputs are the water flow rate Qo and the temperature T which is maintained around an operating point. The temperature of the water Ti is assumed to be constant. The objective of the thermal process is to provide a constant water flow rate at a given temperature. Using the hydraulic and thermal equations, the system can be described by the following equations:

V

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )        × - - = - = t Sh t q T t T t CSh t P dt t dT t q t q dt t dh S i i i ρ (13)
where C ρ represents a constant thermal variable and i T is equal to 20°C.

Based on previous equation, a discrete state space representation of the system around an operating condition ( op h =0.6m, op T =50°C) is determined as follows:

( ) ( ) ( ) ( ) ( )      = + = + k x k y k u B k x A k x d d 1 ( 16
)
where the output vector y is equal to [ ] T h T and the input vector u defines [ ] T i P q

. It should be noted that the sampling period is fixed to 360s in order to respect the closed-loop time constants.

Moreover, the measured output flow rate can be determined by using the Torriceli-rule as:

( ) ( )

k h k q o η = ( 19
)
where η defines the outflow coefficient.

In this study, only sensor failures are considered: level sensor H, output temperature sensor T and output flow rate sensor Qo. As indicated previously, several methods have been proposed in the literature to generate structured residuals and to perform the fault diagnosis [START_REF] Isermann | Trends in the application of model based fault detection and diagnosis of technical processes[END_REF]. However in this example, due to the property of matrix A (which is diagonal), structured residuals can be generated directly with a conventional observer: each residual is sensitive to one fault. Then based on the state space representation, a conventional Luenberger observer is considered to generate the residuals vector ( ) ( )

k y k y - , such that: ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )      = - + + = + k x k y ˆk y k y K k u B k x  k x ˆd d 1 ( 24
)
Based on the residuals vector

( ) ( ) [ ] ( ) ( ) ( ) ( ) [ ] T T k h k h k T k T k r k r - - = 2 1
evaluation, the symptoms vector

( ) ( ) [ ] T k u k u 2 1
is performed in order to detect fault occurrence on H level sensor or T temperature sensor.

Moreover, according to the physical equation between output flow rate Qo and liquid level H, a residual

( ) ( ) ( ) [ ] k q k q k r o o - = 3 is calculated.
Based on the residual evaluation of

( ) ( ) ( ) [ ] T k r k r k r ) k ( r 3 2 1

=

, the associated fault incidence matrix is defined in the following table :   Table 10: Incidence matrix.

H Qo T u 1 0 0 1 u 2 1 0 0 u 3 1 1 0 4.

Dynamic BN design for decision-making

The proposed approach has been designed with the help of the software BayesiaLab (www.bayesia.com). The incidence matrix, defined in Table 10, leads to a DBN model presented in Figure 9. For all faults of the system, the probability of missing detection is assumed to be fixed to 0.02 and the probability of false alarms is equal to 0.05. The probability of missing detection with two simultaneous faults is fixed to 0.01. Consequently as presented in §3.1 (see Table 4) the CPT of u 1 , is defined Table 11, and also the CPT of u 3 (see Table 5) is defined Table 12. In order to define the dynamic reliability model, Figure 10 to Figure 12 present the MC and the Mean Time To Failure (MTTF) which is used to determine the failure rates λ. This failure rate quantifies the transition between the states of 3 considering faulty components and associated probabilistic state matrix X P defined in eq. ( 9). The Markov Chains of the components are supposed to be independent. It should be noted that two states {up, down} are considered for sensors Qo and H, but one more state (dgd) is considered for sensor T which corresponds to a degraded state of the component. As defined in §3.2., the CPT used to simulate the MC for the sensor H (resp. T) reliability described in Figure 11 (resp. Figure 12) is presented in Table 13 (resp. Table 14). The state (dgd) is a functioning state so no fault occurred in this state as implemented in the CPT presented in Table 15.

MTTF 1 =45 000 h λ 1 =0. 22 10 -4 up down λ 1

Results and comments

Based on the incidence matrix (see Table 10) and under any assumptions of the number of faults, then if the coherence vector issued from the residual evaluation at sample k is equal to [0 0 1] T or to [1 1 1] T , for example, the fault indicators I generated by a logic test is as in Table 16.

Table 16: Fault indicators.

U I H I T I Qo [0 0 1] T 0 0 1 [1 1 1] T 1 1 1
Because H and T fault signatures are different, and Qo fault signature is included in H fault signature, the fault isolation is not easy to perform:

when coherence vector is equal to [0 0 1] T , the decision-making provides a fault isolation on sensor Qo (I Qo =1) then a maintenance action is performed to repair this sensor; -If the coherence vector is equal to [1 1 1] T , the three sensors are suspected to be down with the same possibility ((I H = I T =I Qo =1).

However, based on our approach, it could be possible to optimize the maintenance action.

In order to illustrate the performance and the limitation of the proposed method, various faults scenarios have been considered: Scenario A) A false alarm occurs at sample k=6 which appears as an outlier on sensor T.

Scenario B)

A bias on the sensor Qo is assumed to occur at sample k=9 and to repair at sample k=14.

Scenario C)

The system is in a fault-free case.

Scenario D)

T and H sensors faults are supposed to occur simultaneously at sample k=26.

The dynamic behavior of the structured residuals vector has been illustrated in the presence of the various faults scenarios Figure 13.

Scenario A) and Scenario B)

The residual sensitive to a fault is affected. The other residuals are close to zero.

Scenario C)

The system is in a fault-free case. All residuals are close to zero.

Scenario D)

According to the fault incidence matrix issued from the structured residuals (see Table 10), all residuals are different from zero at sample k=26 when the sensor faults occur simultaneously. These residuals are evaluated using the statistical test and are detected isolated correctly as shown in Figure 14: Scenario A) A false alarm occurs at sample k=6 which appears as an outlier on the first symptom which switch to "1" during one sample.

Scenario B)

According to the structured residuals defined in the incidence matrix (see Table 10), only the third symptom switched to "1". Few samples after the third symptom switch to "0" due to a maintenance action. The failure probabilities for the three sensors are presented at Figure 15 without taking into account the dynamic reliability of components. Figure 16 is devoted to the illustration method through the failure probabilities evolution including the dynamic reliability of components.

Scenario A)

The outlier generates a false alarm, the CPT for T can only reduce the value of failure probability to 0.9515 computed based on the Baye's theorem according to the false alarm probability α T =0.05 and the missing detection β T =0.02 with a priori probability distribution: p(Fault T= occured) = 0.5 and p(Fault T = not occured) = 0.5 (see Figure 15). However, as shown at Figure 16, the decision-making is based on reliability of the component. Thanks to the Bayes' theorem, the inference algorithm in the BN computes the reliability of the component taking into account the Markov Chain model with failure rate parameters and the online information based on the residual evaluation (symptoms). Therefore, if all the symptoms are known and close to 0 then the reliability of the component is close to 1. However, the fault is suspected to have occurred if the residuals deviate from zero. Another advantage, of taking reliability into account, lies in the fact that if the residuals states are not computed or if no ambiguity appears in the fault signature, and then the knowledge relative to the reliability of the component is dominating. In this scenario, the integration of reliability and the symptoms result in similar way a sliding window and therefore annihilate the false alarm.

Scenario B)

The two BN methods isolate the fault. It could be noted that a time delay is observed for the second one due to reliability consideration. After tree sampling periods, the probability of fault occurrence I Q is equal to one. The probability of the I H fault occurrence is increased. This fault generates deviation on residual u 3 but according to the inference in the BN, this fault is considered as not coherent with the states of the other symptoms. Therefore, the probability of the I H fault occurrence stays close to 0.

Scenario C)

When a maintenance action is realized, the decision-making is back to a fault free case. The probability that the fault I Q occurred is re-initialized to zero. is explained by the failure on sensor H, then the Qo failure probability is computed based on the Baye's theorem, equal to "0.5156" this value translate the uncertainty (see Figure 15). However, according to the DBN, then the Qo failure probability increases by taking into account the reliability of components (Figure 16). A priori the fault I Q has not occurred. Nevertheless, after long delay, the fault is suspected due to the degradation of the component modeled in the DBN from the failure rate of the component. With the proposed method, it is possible to plan a maintenance action without visiting the Q sensor at the first place due the low level of failure probability. Then, the maintenance action can be focused to the others, T and H sensors, showing a higher level of failure probability. Therefore, if all the symptoms are known and close to 0 then the reliability of the component is close to 1. For complex systems, the problem of the decision-making when various fault signatures vectors are identical or similar can be allayed by using a suitable dynamic Bayesian network. The simulation example, a heating water process, has highlighted the performances of the method.

The design of the dynamic Bayesian network requires the false alarms and miss-detection probabilities of the residual evaluation parameters that are not always possible to assess. Nevertheless, the results obtained in this paper are encouraging and allow us to advocate the method in order to optimize the maintenance actions. Therefore, for a system which is liable to various occurring faults simultaneously or which is defined through an incidence matrix with similar fault signatures, the fault probabilities, provided by the method will enable to plan the maintenance actions. In future works we are interesting to analyse the sensitivity of the decision making to the parameter as false alarms and missing detection.

Figure 1 :

 1 Figure 1: Illustration of residual generator. While a single residual is sufficient to detect a fault, a set of residuals is required for fault isolation. Several methods have been proposed in the literature to generate structured residuals and to perform the fault diagnosis [1]. One of the most popular is the observer-based design [12], [13]. The causal knowledge-based addressed in [14] deals with the complexity of large scale systems. Another way to generate structured residuals is to develop models having an appropriate structure. The structured residuals may be generated by parity equations from the input-output model, or balance equations (chap. 6 and 9, [4]). The objective is to decouple the faulty effects

• 1 H

 1 : the residual is corrupted by fault.

Figure 2 :

 2 Figure 2: Illustration of the error of the first kind α and error of the second kind β .

Figure 2 or 1 H

 21 Figure 2 gives an illustration of hypotheses testing problem with Gaussian density functions. The density functions ) ), ( ( 0 H k r f and ) ), ( ( 1 H k r f characterize the probability that the residual ) (k r respect the hypothesis

Figure 4 :

 4 Figure 4: Two structured Bayesian networks.

  a posteriori distribution of probability over the fault states according to the states of the symptoms. ( ) n F p is the a priori distribution of probability over the fault states and ( )

  not d. (resp. not o.) denotes "not detected" (resp. "not occurred"). The definition of false alarms and missing detections are represented Figure5.

Figure 5 :

 5 Figure 5: Illustration of the error of the first kind j α and error of the second kind 1 β and N β .

1 is

 1 generated according to the incidence matrix defined in Figure4b.

  to the missing detection rate β 1 (resp. β N ) and the false alarm rate α j for the fault F 1 (resp. F N ).

Figure 8 :

 8 Figure 8: Heating water process.

Figure 9 :

 9 Figure 9: Graphical model of the decision-making with DBN.

Figure 10 : 1 Figure 11 : 3 Figure 12 :

 10111312 Figure 10: Reliability MC model for sensor Qo.

Figure 13 :

 13 Figure 13: Residuals behavior.

Figure 14 :

 14 Figure 14: Symptoms scenario.

  During this period, no fault occurs. Symptoms are equal to "0".Scenario D)T and H sensors faults are supposed to occur simultaneously. Based on their fault signatures, all symptoms switched simultaneously to "1".

Figure 15 :

 15 Figure 15: Failure probabilities for the three sensors with BN.

Figure 16 :

 16 Figure 16: Failure probabilities for the three sensors with DBN and reliability.

  5. CONCLUSIONThis paper presents a new strategy to increase the performance of the decision-making in model-based fault diagnosis. The developed approach consists in taking into account in fault diagnosis scheme a priori knowledge on the faulty or non faulty system by a Markov chains modelling. Thanks to the Bayes' theorem, the inference algorithm in the BN computes the reliability of the component taking into account the a priori Markov Chain model with failure rate parameters and the on-line information based on the residual evaluation (symptoms).

  and the arcs indicate direct probabilistic relations between the nodes. A discrete random variable X is represented by a node n with a finite number of mutually exclusive states n

						m s . A set of states is
	defined as	n	{ s ,... n : 1	s	n M	}

They are directed acyclic graphs used to represent uncertain knowledge in Artificial Intelligence [19]. Each node represents a discrete variable defined over several states

Table 2 :

 2 Conditional Probability Table for n j .

	3, nodes i n and j n are linked by an arc where i n is considered as a parent of j n .
				S	i n	i n s ,... : 1	s	n i n i n M and { } j j j n L n n s s ,... : 1 S . The CPT of j n is then calculated by the conditional probabilities
	p	(	n	i j n	)	over each j n state according to its parents i n states as presented in Table 2.

j Figure 3: Elementary Bayesian network.

A conditional probability distribution quantifies the probabilistic dependency between j n and its parent i n and is defined through a Conditional Probability Table (CPT). Therefore, the nodes i n and j n are defined over the sets { }

Table 11 :

 11 CPT of the node u 1 .

		u 1
	Fault T	not detected detected
	not occurred	95	5
	occurred	2	98

Table 12 :

 12 CPT of the node u 3 .

	u 3

Table 13 :

 13 CPT of the node H_sensor(k).

		H_sensor(k)
	H_sensor(k-1)	up	down
	up	99.98	0.02
	down	0	100

Table 14 :

 14 CPT of the node T_sensor(k).

			T_sensor(k)
	T_sensor(k-1)	up	dgd	down
	up	99.985	0.012	0.002
	dgd	0	99.967	0.033
	down	0	0	100

Table 15 :

 15 CPT of the node Fault_T.

		T_sensor(k)
	T_sensor(k)	up	down
	up	100	0
	dgd	100	0
	down	0	100

Dynamic model of reliability: a dynamic Bayesian network solution [20]

In order to model dynamic behaviour of the system performances degradation, dynamic Bayesian network has been considered. Let us recall some fundamental Markov Chain (MC) model.

In the framework of decision-making, a discrete random variable X with two states {up, down} is considered.

These states represent respectively the operational and failure state of the component. Associated to a discrete random variable X , a matrix X P defines the probabilistic state transitions between (up) and (down):

The transient analysis of the MC based on the Chapman-Kolmogorov equation [START_REF] Cassandras | Introduction to Discrete Event Systems[END_REF] provides an expression for

Under dynamic consideration in a Bayesian network, the state of the variable i X (to model the i th component of a system) is represented at sample k by a node n i (k) with a finite number of states.

)) ( ( k n p i denotes the probability distribution over these states at time step k. The dynamic Bayesian networks allow to represent random variables and their impacts on the future distribution of other variables [START_REF] Murphy | Dynamic Bayesian Networks: Representation, Inference and Learning[END_REF]. Starting from an observed situation at sample k=0, the probability distribution

)) ( ( k n p i over M states for the component i X associated to the node i n is computed by the dynamic Bayesian network inference.

Indeed, it is possible to compute the probability distribution of any variable i X at sample k based on the probabilities defined at sample k-1 as shown in the elementary network presented in Figure 6. The first slice contains the nodes corresponding to the current time step k-1, the second one contains the nodes of the following time step k. Observations, introduced as hard evidence or probability distributions, are only realised in the current time slice. The time increment is carried out by setting the computed marginal probabilities of the node at sample k as observations for its corresponding node in the previous time slice. The CPT in dynamic Bayesian network is equivalent to X P (Table 6).

Fusion of the Incidence Matrix and the dynamic Bayesian network model of reliability: a solution for increasing effectiveness of decision-making

As presented previously, the dynamic Bayesian network models the component reliability which takes into account the time degradation of components. The representation of incidence matrix as a graph, defined at the beginning of the paper, provides a formalism, which is used to realize the fusion between fault diagnosis and reliability model. The decision-making is established after fusion of information obtained from both residual analyses and reliability estimation. Therefore, based on the Bayesian network representation, image of the incidence matrix, the dynamic evolution of the component reliability is taken into account on each node F i as follows: Figure 7. This relationship involves the definition of a CPT. The CPT of F i is very simple to define if the component reliability is modelled with the states (up) or (down) which is a common case in fault diagnosis. Then identity matrix is used to compute the distribution on F i as presented in the Table 7. Based on the elementary example, illustrated in Figure 6, the computation of ( )

is performed thanks to the inference algorithm in dynamic Bayesian network a simplification of the formulation is presented by:


Gareth M Jones
email: g.m.jones@warwick.ac.uk

Mattia Biesuz
email: mattia.biesuz@unitn.it

Wei Ji
email: jiwei@whut.edu.cn

Sandra Fisher
email: sandra.fisher-john@lucideon.com

Carolyn Grimley
email: carolyn.grimley@lucideon.com

Charles Manière
email: charles.maniere@ensicaen.fr

Claire E J Dancer
email: c.dancer@warwick.ac.uk

Promoting microstructural homogeneity during flash sintering of ceramics through thermal management

Flash sintering (FS) is a novel field-assisted sintering technology, where the ceramic is heated internally by the Joule effect. While FS promises a tremendous reduction of ceramic firing time and furnace temperature, it has been applied only at the laboratory scale to date. The key limitation of scaling up the technique to the industrial manufacturing level is the intrinsic difficulty managing the heat generation and obtaining homogenous microstructures in components of industrial interest. Heterogeneous regions primarily originate from the different types of thermal gradients that develop during FS; therefore, the management of heat generation is crucial to achieve uniformity. In this article, we discuss the advantages of controlling the microstructural homogeneity of ceramics during FS, and the technical routes to achieve this. The origin and formation mechanisms of thermal gradients upon flash sintering are outlined. Possible approaches to reduce thermal and microstructural gradients are identified. The opportunities and challenges in scale-up of FS are discussed from both industrial and scientific perspectives.

Introduction

The development of energy-efficient firing technologies has driven research on sintering throughout the twentieth century, with both economic and environmental motivations. Several novel low-energy consolidation techniques have been developed in the past decade, of which flash sintering (FS) is one of the most promising.1-3 FS is a field-assisted sintering technology where the simultaneous application of electric fields and heat to a ceramic green component leads to a sudden drop of its electrical resis-tivity, an internal heat generation by the Joule effect,4,5 and densification within a few seconds. Attempts have been made to scale-up the flash process using other well-established sin-tering equipment such as flash-spark plasma sintering (flash-SPS)6-8 and flash-microwave sintering,9,10 though to date, the technique primarily remains at the laboratory scale. In addition to reducing energy consumption, FS ena-bles other enhancements to glass and ceramic manufactur-ing. 11 The flash event can be used to facilitate the synthesis of complex oxides,12 formation of ceramic-ceramic13 and metal-ceramic14,15 joints, and viscous16 and plastic17 defor-mations. The flash event has also allowed the consolidation of microstructures18 and phases19 that cannot be achieved via conventional heating. Additionally, it has been shown to introduce a significant perturbation of the crystallographic and electronic defect chemistry.11,20,21

Several types of thermal gradient develop during FS as heat is internally generated by the Joule effect.4,5 These gradients can occur over all length scales characteristic of the polycrystalline ceramics, from the macroscale bulk, to the microstructural features such as the grain size scale and the porosity scale, and the finest scales such as the grainboundary and other interfaces.2 Such gradients can lead to undesired microstructural inhomogeneities in the sintered component. For the industrial adoption of FS in ceramic manufacturing, a key technological challenge is the fine control of current flow and hence heat transfer in the ceramic component dur-ing FS. In this article, different types of thermal gradients are described and their origins clarified. This is followed by a comprehensive discussion of solutions to avoid thermal and microstructural gradients.

Temperature and microstructural heterogeneity

Origin of microstructural inhomogeneities

FS experiments are often carried out on small specimens to minimize thermal heterogeneity. However, the use of samples of even a few millimeters in thickness has been reported to generate temperature heterogeneities3 ranging over a few hun-dred degrees Celsius in some cases. Evidence of this comes from thermal images,22-24 microstructural variations,25 and X-ray diffraction peak width broadening. 26 The origins of the thermal gradients and their consequent microstructural inho-mogeneities are indicated in Figure 1. We can outline two main categories of thermal gradient: in the sample cross sec-tion (i.e., orthogonal to the current flow) and along the gauge length (i.e., in the direction of the current flow). In the fol-lowing sections, we expand on this categorization and discuss how they enable the development of different microstructural variations.

Thermal inhomogeneities and gradients

Thermal gradients in the sample cross section upon FS arise cooling fluxes. The negative temperature coefficient (NTC) resistivity of most ceramics implies a higher current concentration and thermal dissipation in the hot zones. 4 The developed surface-to-core temperature differences depend on the applied electric power, the thermal conductivity, the sample geometry, and the activation energy for electrical conductivity. This NTC phenomenon, combined with the specimen cooling at the boundaries, tends to amplify the temperature gradients and may generate highly unstable heating regions.27,28 These so-called hot spots can be detected in samples after sintering through microstructural features, including local melting,3 and abnormal grain growth.29 Grain core-boundary temperature differences are minimal and unlikely to contribute significantly to the temperature gradients.30

Several researchers have modeled the heat distribution in the samples during FS.31-35 Homogenous heat capacity and transfer are commonly assumed for the bulk material. These studies focus on macroscale effects over the whole sample or on nanoscale effects at particle-particle interfaces, with intermediate scales neglected. A simulation study by Dong based on perturbation methods tests the validity of this assumption and shows the hot spot phenomenon may be active for specimens above a certain critical size ( Zcr ) as low as a few millimeters for yttria-stabilized zirconia (YSZ):36 from a disparity between the rate of Joule heating in the sample and the rate of heat transfer and are initiated by surface

Z cr = 2κ(T )ρ(T )RT 2 = QE 2 2kRT 2 , pQ (1)
where κ(T ) represents the thermal conductivity, ρ(T ) represents the electrical resistivity, Q represents the corresponding activation energy, p represents the electric power density, R represents the gas constant, T represents the absolute temperature, and E represents the electrical field.

The existence of a critical hot spot size is due to the fact that the volumetric heat produced by a temperature perturbation increases with the square of the perturbation radius (λ), whereas the heat dissipated toward the neighboring regions of the specimen is linear with λ. Thus, when λ > λ cr, the volumetric heat is dominant and the perturbation is stable, on the other hand, if λ < λ cr heat is exchanged with the surrounding and the perturbation vanishes. Additional thermal gradients could develop close to microstructural features such as grain boundaries2 and pores.37

Geometry and electrode contact effects

The geometry of the sample and attached electrodes has a significant effect on the development of thermal gradients during FS. The magnitude of a thermal gradient from surface to core depends primarily on the exposed surface area of the sample, power dissipated (both dependent on geometry), and thermal conductivity (material property). Small non-industrially relevant sample geometries such as dogbones have been extensively used in FS studies. Their minimized cross-sectional area reduces the severity of the gradient between the surface and core. For long gauge lengths, anode/cathode and electrode/core gradients are nonnegligible. Field enhancement at the electrode hole can further enhance the thermal gradient.38 Dogbones have typical gauge length of 15-20 mm with a cross-sectional area of 3-6 mm2. 39-42 In the case of disk-shaped samples, which have a larger diameter compared to their height, the ratio between the area in contact with the electrodes and the curved side surface depends on the height-to-diameter ratio with heat dissipation through the electrodes becoming dominant in the case of thin samples. The small height minimizes any gradient between the sample center and the electrode region, but also facilitates the formation of hot spots and radial gradients. 43 Common sizes of disk vary, but most used in FS experiments are less than 10 mm in diameter,44-46 the exception being flash-SPS where larger disk diameters are more common47-51 When "pellet" samples are used (i.e., cylindrical-shaped samples with similar diameter to height),27 geometry effects lie between the two extremes of the dogbone (low diameter-toheight ratio) and the disk (high diameter-to-height ratio).

Poor contact of the sample to the electrode introduces high resistance in the electrode region over which power is dissi-pated. Overheating of the electrode region can produce exces-sive grain growth at both the anode and cathode for direct current (DC) and alternating current (AC).38 Good thermal contact between the sample and electrodes allows heatsinking from the sample to cooler electrodes.27

Material-dependent effects

The bulk material properties affect thermal homogeneity through the thermal conductivity and the electrical resistivity (Equation 1). 36 The activation energy (Q) for electrical conductivity in ceramics varies significantly and this affects the probability of hot spot formation. Hot spots become more stable in materials possessing high Q and low k , or when high power density is applied (Equation 1). The behavior of materi-als undergoing flash sintering depends on the relative values of their electrical and thermal conductivity. The presence of impurities, dopants, or secondary phases in ceramic compos-ites can alter the electrical and thermal conductivity and hence the FS behavior. Composite materials can have complicating effects if the conductivity of the constituents varies as in alu-mina-titania52 and aluminazirconia composites53 Similar effects may occur when using reactive flash. 54 The presence or formation of thermally conductive phases during FS can act as a heatsink or thermal buffer to mitigate hot spot forma-tion, such as the glass phase in porcelain. 55 The packing of the ceramic compact also affects the overall contact resistance and ease of FS, smaller particles pack better with a higher ratio of particle-particle contact to surface area compared to larger particles. 56 Electrochemical effects associated with the application of a DC electric potential have been observed during FS in several ceramic systems. These phenomena are particularly relevant in ionic conductors where the displacement of ions under the effect of the external field leads to variations in the sample composition between the anode (+) and cathode (-). In alkali-ion conductors, such as soda lime silicate glass, FS causes a substantial accumulation of the positive charge carriers at the cathode and a highly resistive alkali depletion layer forms at the anode.57,58 The local increase in the electri-cal resistivity causes a substantial overheating of the anodic region (even > 1000°C), observed using thermal imaging.57,58 Electrochemical phenomena have also been reported in oxygen ion conductors, including YSZ23,59 and gadolinium-doped ceria.20,60,61 Here, the mobile charge carriers are oxygen vacancies, vO , which move and accumulate at the cathode. vO accumulation is somewhat counterbalanced by molecular oxygen reduction from the atmosphere at the cathode, which acts as a partially blocking electrode.

The molecular oxygen reduction rate is highly dependent on the electrode configuration59 and is often not sufficient to sustain the current flow. This leads to a strong electrochemical reduction of the cathodic region, which propagates to most of the gauge length during FS. This phenomenon is associated with visible electrochemical blackening62 due to the formation of donor levels in the bandgap. 23 This anode/cathode asymmetry causes thermal gradients due to imbalance in the concentrations of positive and negative charge carriers, leading to overheating of the anode.23,24 In addition to these thermal effects, cathodic electrochemical reduction induces microstructural gradients due to localized changes in the diffusivity of cations.63,64 While the most extreme electrochemical effects have been observed in ionic conductors, anode/cathode asym-metry, reduction, or blackening has been reported during FS in several ceramic systems.11

Microscale microstructural development

While FS can cause the formation of large inhomogeneous regions on the scale of the sample that are clearly detrimen-tal, on the microscale, it promotes significant homogeneity compared to other sintering techniques. Recently, Ji et al.65 compared the microstructure of 3YSZ produced by different densification technologies: conventional sintering, fast firing, self-propagating high-temperature synthesis, and FS (Figure 2). By heating to different temperatures, both partially and fully dense materials were obtained using each technique. Conventionally sintered (Figure 2a-b) samples showed heterogeneous microstructures with pores of high coordination num-ber originating from the green compact. Pores and abnormal grain growth were also observed in samples produced by fast firing (Figure 2c-d) and self-propagating high-temperature synthesis (Figure 2e-f). Large pores shrink proportionally to the macroscopic shrinkage of the body and are difficult to fully remove. In contrast, the partially dense FS specimen had a fine and uniform porosity, and the fully dense sample was homogeneous (Figure 2g-h). In the FS sample, only small pores with low pore coordination number were observed, demonstrating that FS accelerates the homogenization of the morphology in the microscale. This difference in pore size contributes signifi-cantly to the greater densification rate following rapid heating with or without an electric field.66

Strategies to reduce temperature and microstructural gradients

Modeling thermal gradients

Thermal gradient development is highly influenced by the cavity configuration in FS and can be most easily explored using finite element simulation. 67 "Hybrid" heating uses an additional external heat source to homogenize the thermal field. For "direct" heating of thick specimens, hightemperature differences are expected, even assum-ing a uniform volumetric dissipation (Figure 3a). 36,67 The relative cooling fluxes at the specimen boundaries are cal-culated as ~ 70% radiation, ~ 20% convection, and ~ 10% conduction for the configuration shown in Figure 3a at high temperature. The dominant radiative contribution is easily explained by the T 4 Stefan-Boltzmann Law. Egorov et al. 68 state the convective relative cooling proportion may be higher if the heated specimen is immersed in a large air cavity.

Thermal gradients originating from cooling fluxes due to boundaries can be reduced by three main approaches: 67 first, by lowering the heating cavity air volume using thermal insulation; second, by delaying the cooling core/edge gradient development by imposing ultrarapid heating from the start (or using double-step or multistep FS); and third, by employing hybrid heating using a susceptor or an external heating element. As indicated by Figure 3c, the most efficient way to reduce the thermal heterogeneities in direct heating is a configuration combining thermal insulation and high heating rate. 67 This conclusion is valid in both flash-microwave sintering and "standard" FS. A similar approach was applied to FS by thermally insulating the specimen using coarse zirconia powder or an alumina foil/ wool. This approach increases the homogeneity of the microstructure while significantly reducing the energy consumption. 69 In microwave sintering, a susceptor is used to homogenize and stabilize the heating. 70 Hybrid heating can also be applied to "traditional flash" and microwave flash experiments where the susceptor is electrically heated to follow the heating of the sample. This approach is commonly employed in flash-SPS, where a graphite foil/felt is used to activate the flash transition. This allows the production of specimens with large diameters (30-60 mm) 7,49 and homogeneous microstructures can be obtained. Simulation explains this stability by the lateral graphite elements heating (electrical susceptor) and thermal confinement originating from thermal contact resistance at each interface. 8

Experimental solutions

Power delivery and heating rate

In this section, we discuss possible technical solutions to avoid or mitigate the effects of the formation of thermal gradients during FS.

As the hot spot critical size (Equation 1) decreases when increasing the applied power, the power peak at the voltageto-current control switching facilitates hot spot formation and so has a detrimental effect on the final microstructure. One simple solution is to operate under a low current limit from the start. This "power-controlled" method results in a much slower heating rate, 71 but limits the likelihood of hot spot formation. An alternative is to adapt the two-step conventional sintering process to FS. 25,46 The interruption in current flow reduces instances of current concentration, as when the flash is reinitiated the current will find a different unblocked path. Proportional integral derivative (PID) control produces current ramps with rates varying from 1 A/cm 2 /s to 0.001 A/cm 2 /s. At slow ramp rates, increased grain growth was observed. 72 Using a similar method, different ramp profiles (linear, squareroot, and parabolic) were trialed. 73 Current limit stepping via PID has been built into FS protocols 74 and resulted in superior homogeneity between the surface and core. 75 The use of the current ramp FS approach also resulted in a smaller anode/ cathode grain size difference in 3YSZ. 41 To provide reproducibility and control required to adapt FS to industrial applications, Lucideon Ltd. has developed a control system based on a nonlinear algorithm. In conjunction with a real-time computer, the software controls current, voltage, or power outputs on a millisecond basis, preventing any value from overshooting its set point. Controlled ramps of any the previously discussed parameters and pulsed power dissipation can be achieved, and the system functions at frequencies ranging from DC up to the megahertz range. [START_REF] Lucideon | Flash Sintering -Technology Scale-up for Advanced Manufacturing[END_REF]

Electrode configurations

The electrode material and interface paste composition (if used) can greatly affect the contact resistance in FS. Several different electrode configurations have been proposed to solve or alleviate electrode effects. The effect of electrochemical redox under DC can be mitigated through the use of porous mesh electrodes and pastes 59 (i.e., electrodes that facilitate O 2 diffusion and its incorporation at the cathode). Another solution is the use of molten salt electrodes, which can supply ionic species to the depleted electrode region in alkaliion conductors, 28 though the use of AC fields is more generally applicable. 23,[START_REF] Mclaren | [END_REF] The adoption of an optimized electrode shape 78,79 can also alter the thermal gradient. A solution to enable continuous production has been developed by a collaboration between the University of Colorado and Lucideon Ltd. using sliding electrodes above and below a whiteware sample within a rolling kiln. 80 Contactless FS opens up opportunities for components with unfavorable geometries, including thin films for solid-oxide fuel cells and solid-state batteries, environmental and thermal barrier coatings, as well as complex geometries with nonuniform thickness. Contactless FS may avoid the formation of thermal gradients caused by electrode-specimen contact as it delivers current to the sample via a conductive gaseous medium between the electrode and the sample. Successful trials of contactless FS have used electric arcs, cold plasma, and local heating using an electrically conductive flame. [81][82][83] Lucideon Ltd. has developed a demonstration scale system for contactless flash sintering. [START_REF] Lucideon | Flash Sintering -Technology Scale-up for Advanced Manufacturing[END_REF] The system consists of a multiaxis robot capable of full dimensional rastering allowing homogeneous flash sintering of large areas and thin (< 1 mm) geometries. The system has also been used for localized roomtemperature repair of glaze defects on sanitary ware.

Additives

From Dong's equation (Equation 1), at a given temperature, the critical size of a hot spot rises as the activation energy for electrical conductivity is reduced. This can be achieved by adding electrically conductive phases to the ceramic, including carbon nanotubes, 84 graphene, [85][86][87] or graphite. 88 The addition of these compounds also increases the overall thermal conductivity of the specimen, which further decreases the hot spot formation probability and reduces the surface/core thermal gradients. Preliminary results have also shown that the addition of water ("cold FS") allows flash of thin pellet-like samples at room temperature without hot spot formation. 89

Summary

We have described how a variety of thermal gradients form during flash sintering, leading to microstructural inhomogeneities. A range of strategies have been proposed to mitigate and manage this formation; these are summarized in Table I.

It is clear that the strategies chosen to implement FS on an industrial scale will be shaped by the requirements of each product and the benefits pursued by its manufacturer. For instance, small additions of a more conductive phase will lead to homogeneity improvement by mitigating the thermal gradient between the sample core and edges. At the same time, this tunes the optical and electrical properties of the final product, which may be detrimental for some applications. While it is unlikely that any single solution would be universally appropriate, some inhomogeneity mitigation strategies are more amenable to large-scale production.

Contactless FS offers the most extensive set of benefits, by removing the need for expensive electrode materials, mechanisms to ensure contact is maintained (conductive paste and/ or uniaxial pressure), and the limitation to geometries, which possesses a regular cross section and provides two flat surfaces for contact. Currently, contactless has been successfully used for applications requiring thin ceramic layers, such as thermal barrier coatings. 82 The next most practical solutions are the use of AC and electrical control systems. In the laboratory, uncontrolled DC is preferred due to its simplicity. For a manufacturer, any additional difficulty introduced by controlled AC is small compared to the overall cost of retrofitting a line for FS and is beneficial for stability and homogeneity. These two strategies can greatly reduce hot spot formation and electrochemical reduction effects.

Reducing the macroscopic thermal gradients induced by heat loss at free surfaces and electrodes is key to the adoption of FS as an industrial-scale technology. Continued work on the effect of alternate paste compositions and morphologies would aid in the development of flash for more specialized batch processes. However, continuous throughput processes are not readily amenable to the use of pastes or physical insulation. Thermal homogeneity through the use of insulation has been attempted by researchers developing a flash-SPS hybrid. 90 Similar effects may be applicable to microwave sintering or microwave flash, which can display a similar thermal runaway mechanism. 67 Future efforts should explore both engineering solutions for insulation and alternate creative technical approaches.

Figure 1 .

 1 Figure 1. Different types of thermal and microstructural gradients that can develop during flash sintering. In the cross section, both hot spots (small regions within the bulk) and surface/core gradients can develop, resulting in different regions of microstructural inhomogeneity. Differences in microstructure along the gauge length are due to the electrode composition and the electric contact between electrode and sample. The microscope images show examples of these effects for (a) cross-sectional differences in regions in a flash sintered sample of 8 mol% yttria-stabilized zirconia25 and (b) differences between core, surface, and electrode regions of a flash sintered sample of 3 mol% yttria -stabilized zirconia.27 On the right (c) is an example of a thermal imaging measurement of the thermal gradient in soda lime silicate glass during flash sintering.28

Figure 2 .Figure 3 .

 23 Figure 2. Morphology of partially and fully dense samples upon (a, b) conventional sintering, (c, d) fast firing, (e, f) self-propagating high-temperature synthesis, and (g, h) flash sintering65.

Table I . The origins of thermal gradients and strategies to mitigate their formation.

 I

	Type of gradient	Origin	Materials	Strategies
	Surface/core (core overheated with	Heat radiation from the hot surface	• Materials with low thermal conductiv-	• Addition of a thermally conductive phase
	respect to surface)		ity	• Ultrafast (very small time scale)/pulsed
				• Decrease heating rate
				• Use of susceptor (hybrid heating)
				• Use of thermal insulation
	Hot spot	Thermal instability related to the	• Materials with low thermal conductiv-	• Addition of a thermally conductive phase
		negative temperature coefficient for	ity	• Addition of a conductive phase with
		electrical resistivity	• Materials with negative temperature	low activation energy for electrical
			coefficient for electrical resistivity	conductivity
			• Materials with heterogenous con-	• Decrease of the power peak (step or
			ductivity	current rate flash sintering)
				• Increase furnace temp (tune to material
				conductivity)
				• Decrease heating rate or ultrafast (very
				small time scale)/pulsed
				• Continuous flash processes (contact or
				contactless) (variance of current path-
				way/limited duration of heating)
				• AC flash (variance of current pathway)
	Electrode/center	High contact resistance between	• Any	• Use of conductive pastes or pressure to
	(Electrode region overheated with	electrodes and sample		reduce the contact resistance
	respect to the center)			• Electrodes with higher heat capacity
				• Matching electrode positive temperature
				coefficient behavior to temperature
				• Matching AC frequency to charge
				transfer
				• Contactless flash sintering
	Center/Electrode	Electrodes act as heatsinks	• Any	• Electrodes with lower heat capacity
	(Core overheated with respect to the			• Contactless flash sintering
	electrode region)
	Anode/cathode gradient	Electrochemical redox at DC/low	• More evident in ionic conductors	• AC flash sintering
		frequency AC		• Use of nonblocking electrodes
				• Reintroduce lost ionic species (atmos-
				phere/molten salt)

Acknowledgments

The authors thank D. Pearmain of Lucideon Ltd. for helpful discussions and support. Funding for the PhD RD&I Associate for G.M.J. was provided from Lucideon Ltd. and the European Regional Development Fund as part of the "Smart Energy Network Demonstrator" based at Keele University (Grant Reference 32R16P00706). C.E.J.D. acknowledges funding from the Engineering and Physical Sciences Research Council (EP/ R029873/1).


Etienne Genier

Sacha Grelet

Rasmus D Engelsholm

Patrick Bowen

Peter M Moselund

Ole Bang

John M Dudley

Thibaut Sylvestre

An ultra-flat, low-noise and linearly polarized fiber supercontinuum source covering 670 nm-1390 nm

An ultra-flat, low-noise and linearly polarized fiber An ultra-flat, low-noise and linearly polarized fiber supercontinuum source

Supercontinuum (SC) generation in all-normal dispersion (ANDi) photonic crystal fibers (PCF) has recently been investigated due to its ability to generate a low-noise broadband coherent spectrum [START_REF] Hori | Flatly broadened, wideband and low noise supercontinuum generation in highly nonlinear hybrid fiber[END_REF][START_REF] Hooper | Coherent supercontinuum generation in photonic crystal fiber with all-normal group velocity dispersion[END_REF][START_REF] Tarnowski | Polarized all-normal dispersion supercontinuum reaching 2.5 m generated in a birefringent microstructured silica fiber[END_REF][START_REF] Hartung | Design of all-normal dispersion microstructured optical fibers for pulse-preserving supercontinuum generation[END_REF][START_REF] Gonzalo | Polarization noise places severe constraints on coherence of all-normal dispersion femtosecond supercontinuum generation[END_REF][START_REF] Genier | Amplitude noise and coherence degradation of femtosecond supercontinuum generation in all-normal-dispersion fibers[END_REF], with a noise level much lower than typical SC sources using anomalous dispersion. The main issues with typical soliton-based SC generation is the poor flatness and the coherence degradation for a soliton number higher than 16, associated with large pulse-to-pulse fluctuations due to noise amplification by modulation instability (MI), soliton fission and rogue wave generation [START_REF] Dudley | Supercontinuum generation in photonic crystal fiber[END_REF]. On the contrary, ANDi SC generation, which is essentially based on nonlinear coherent effects such as self-phase modulation (SPM) and subsequently optical wave breaking (OWB), has a high degree of coherence and pulse to pulse stability [START_REF] Finot | Beneficial impact of wave-breaking for coherent continuum formation in normally dispersive nonlinear fibers[END_REF][START_REF] Heidt | Limits of coherent supercontinuum generation in normal dispersion fibers[END_REF]. ANDi SC generation is thus a topic being extensively investigated for several applications such as optical coherent tomography [START_REF] Jensen | Noise of supercontinuum sources in spectral domain optical coherence tomography[END_REF], coherent amplification [START_REF] Heidt | Low noise all-fiber amplification of a coherent supercontinuum at 2 m and its limits imposed by polarization noise[END_REF], nonlinear microscopy [START_REF] Dwapanyin | Generalized spectral phase-only time-domain ptychographic phase reconstruction applied in nonlinear microscopy[END_REF], hyperspectral stimulated Raman scattering microscopy [START_REF] Abdolghader | All normal dispersion nonlinear fibre supercontinuum source characterization and application in hyperspectral stimulated raman scattering microscopy[END_REF] , and coherent anti-Stokes Raman scattering (CARS) spectroscopy [START_REF] Viljoen | Implementation of temporal ptychography algorithm, i 2 pie, for improved single-beam coherent anti-stokes raman scattering measurements[END_REF]. Several recent studies have however shown some limitations on the input pulse duration due to stimulated Raman scattering (SRS) and polarization modulation instability (PMI) [START_REF] Gonzalo | Polarization noise places severe constraints on coherence of all-normal dispersion femtosecond supercontinuum generation[END_REF][START_REF] Genier | Amplitude noise and coherence degradation of femtosecond supercontinuum generation in all-normal-dispersion fibers[END_REF][START_REF] Heidt | Limits of coherent supercontinuum generation in normal dispersion fibers[END_REF][START_REF] Møller | Intensity noise in normal-pumped picosecond supercontinuum generation, where higher-order raman lines cross into anomalous dispersion regime[END_REF][START_REF] Feehan | Decoherence due to xpm-assisted raman amplification for polarization or wavelength offset pulses in all-normal dispersion supercontinuum generation[END_REF]. Specifically, it has been shown that for input pulse duration longer than ∼ 500 fs, SRS-induced noise amplification impacts on the SC pulse coherence and spectral flatness [START_REF] Heidt | Limits of coherent supercontinuum generation in normal dispersion fibers[END_REF]. Furthermore, for input pulse duration longer than ∼ 120 fs, PMI can take place between the two orthogonal eigenmodes of the fiber, thus degrading the relative intensity noise (RIN) and the output polarization for non-polarization maintaining (PM) fibers [START_REF] Gonzalo | Polarization noise places severe constraints on coherence of all-normal dispersion femtosecond supercontinuum generation[END_REF].

In this work, we overcome these limitations by using a silica-based PM ANDi PCF pumped with a low-noise, short pulse duration (180 fs) compact ytterbium mode-locked laser at 1049 nm. Using this system, symmetrically flat and ripple-free SC generation has been generated from 670 nm to 1390 nm with a maximum output power of 720 mW. It is shown in particular that pumping on the fast axis of the fiber provides the flattest and most stable SC while removing the cross-phase modulation instability (XPMI) that has recently been evidenced in PM ANDi fibers [START_REF] Genier | Cross-phase modulation instability in pm andi fiber-based supercontinuum generation[END_REF]. A polarization extinction ratio (PER) of 17 dB and an average RIN of 0.54 % from 700 nm to 1100 nm were further measured. We finally compare the spectrally-resolved RIN of our ANDi-fiber-based SC source with those recently published in the literature and we show excellent noise and polarization characteristics. The experimental setup used to observe and analyze coherent SC generation in the PM-ANDi PCF is shown in Fig. 1. As a pump laser, we used a femtosecond ytterbium solid-state mode-locked laser (Origami 10 HP, NKT Photonics) centered at 1049 nm, delivering a 180 fs pulse train at an 80 MHz repetition rate with a maximum average power of 4 W, a low RIN of 0.5%, and a PER of 31 dB. The laser output power was controlled using a variable neutral density (ND) filter. A half-wave plate was used to rotate the field polarization state at the fiber input. An aspheric lens with a focal length of 6 mm was used to couple the light into the fiber core. The fiber is a PM ANDi PCF from NKT Photonics (Model NL-1050-NEG-PM). This fiber has a core diameter of 2.4 m, a relative hole diameter of d/Λ = 0.45, a small hole-to-hole pitch of Λ = 1.44 m, and a nonlinear coefficient of =26.8 W -1 km -1 at 1040 nm. The fiber's dispersion and loss profiles are shown as an inset in Fig. 1 [START_REF] Genier | Cross-phase modulation instability in pm andi fiber-based supercontinuum generation[END_REF]. The dispersion has an all-normal parabolic profile with a peak at -13 ps/nm/km at 1040 nm. The PM property of the PCF is mainly stress-rod induced and it has a slightly elliptic core, which together gives it a high birefringence of 4.7 10 -4 at the pump wavelength. The output spectrum is recorded using an optical spectrum analyser (OSA) and analysed as a function of input polarization and pump power, while the RIN was measured in another path using a monochromator as filter passband, two different fast photodetectors to cover the whole SC bandwidth and an oscilloscope. Figure 2 shows the spectral intensity profile for 3 selected polarizations and an average output power from the fiber of 720 mW (estimated input peak power of 48 kW, including losses and Fresnel reflections). We obtain the broadest bandwidth while pumping along the slow axis (red curve), with a spectrum going from 685 nm to 1395 nm at -3 dB if we disregard the oscillations in the central part of the spectrum. We find the spectrum becomes more noticeably modulated as we move closer to the center of the spectrum. Furthermore, the modulations appear together with the observation of an XPMI peak [START_REF] Genier | Cross-phase modulation instability in pm andi fiber-based supercontinuum generation[END_REF] when pumping along the slow axis. This clarifies the importance of avoiding XPMI, since the modulation can be detrimental to interference based imaging. The spectrum is the narrowest when pumping at 45 • from the principal axes (pink curve) spanning from 700 nm to 1370 nm at -3 dB, which is expected because only half of the power is available on each axis. We observe that the spectrum for 45 • is highly modulated from 870 nm to 1120 nm. The XPMI sideband is again present and slightly stronger than the one obtained on the red curve, as expected from the theory [START_REF] Genier | Cross-phase modulation instability in pm andi fiber-based supercontinuum generation[END_REF][START_REF] Kudlinski | Simultaneous scalar and cross-phase modulation instabilities in highly birefringent photonic crystal fiber[END_REF][START_REF] Drummond | Cross-phase modulational instability in high-birefringence fibers[END_REF]. Then, we obtain an extremely flat spectrum when pumping on the fast axis (blue curve) with a spectrum spanning from 670 nm to 1390 nm at -3 dB without any trace of an XPMI sideband around 1400 nm. Fig. 2: Experimental SC spectra at the output of the PM-ANDi PCF for output average power of 720 mW with an input pulse duration of 180 fs while pumping along the slow axis (red curve), fast axis (blue curve) and at 45 • from the axes (pink curve). The blue and red curves are offset by 15 dB and 30 dB from the pink curve, respectively.

We now investigate the case with pumping along fast axis, for which we observed the best spectrum with the least ripples. Thus, Fig. 3 shows the evolution of the SC generation for a fiber output average power increasing from 1 mW to 720 mW (peak power going from 8 W to 48 kW). We observe two SPM fringes on the pink curve (for a power of 6 mW). Then, by increasing further the power up to 30 mW (brown spectrum), we can clearly see the appearance of OWB on the edges. Furthermore, the SC bandwidth goes from 9 nm (in the linear regime) to 720 nm (going from 670 nm to 1390 nm) at -3 dB, which is more than one octave with a power spectral density (PSD) higher than 0.4 mW/nm (-4 dBm/nm). Next, as expected from coherent ANDi SC generation, the spectrum broadens symmetrically and is extremely flat (the maximum amplitude of the ripples in the SPM area is 1.65 dB). Finally, we do not observe any SRS components in the final spectra at high power (red and black spectra) across the whole bandwidth [START_REF] Gonzalo | Polarization noise places severe constraints on coherence of all-normal dispersion femtosecond supercontinuum generation[END_REF][START_REF] Heidt | Limits of coherent supercontinuum generation in normal dispersion fibers[END_REF]. Thus, we can expect this SC generation to be extremely low-noise and highly stable. Fig. 3: Experimental SC spectra at the output of the PM-ANDi PCF for output average power of 1 mW to 720 mW with an input pulse duration of 180 fs and an input beam polarized along the fast axis.

Using a polarizer at the output of the ANDi-PCF, we measured the PER of the generated SC using the OSA. Fig. 4 depicts the SC spectra when the input beam is polarized along the fast axis and the output polarizer is aligned either on the fast axis (red curve) or on the slow axis (blue curve). Specifically, the polarizer is rotated to measure the minimum power including the whole spectrum and then the spectrum is measured at this angle and at a 90 degree rotation. This way, we can measure an average PER of 17 dB over the full SC bandwidth, which means that the SC output is nearly linearly polarized. We noticed however a PER degradation of 14 dB compared to the input pump laser (PER of 31 dB). This degradation of 14 dB is also measured at low power (no spectral broadening) and may be explained by the collapse of the air holes at the entrance of the fiber required to remove back-reflection in the laser cavity. Fig. 4: Experimental SC spectra when pumping on the fiber's fast axis and using a polarizer aligned either on the fast axis (red) or on the slow axis (blue). The difference between the flat parts of the spectra gives a PER of about 17 dB.

Let us now turn our attention to the noise properties of the pump and the generated supercontinuum. For the pump noise measurement, we coupled the output beam onto a silicon photodiode connected to a fast oscilloscope (DSA91204A from Agilent). After recording 10000 pulses, we obtained the histogram depicted in Fig. 5. Using a Gaussian fit of this distribution, we extract a mean value of 2.5 and a standard deviation of 0.011. Thus, calculating the ratio of the standard deviation to the mean, we obtain a RIN value of 0.44 % using this method, in rather good agreement with the laser manufacturer data (0.5 %). After measuring the pump laser noise, we further investigated the RIN over the whole SC bandwidth. To this end, we used a monochromator (from Princeton instrument) as a tunable spectral filter and the spectrally-sliced SC was sent to the photodetector and oscilloscope. We specifically used two photodetectors to cover the SC bandwidth: a silicon photodetector from 650 to 1100 nm and an InGaAs photodetector for longer wavelengths. The SC RIN was still measured using the Gaussian fit method by saving 10000 pulses onto an oscilloscope (PicoScope) every 10 nm, with a 3 nm bandwidth for the silicon photodetector and a 1 nm bandwidth for the InGaAs one. Fig. 6 depicts the two spectrally-resolved RIN measurements as red dots (silicon) and as green dots (InGaAs), while the blue curve shows the related SC spectrum while pumping the fiber at 45 • off the axes. We find an average RIN from 700 nm to 1100 nm down to 0.54 %, when we omit the high edge RIN values. The results also show that the RIN increases drastically on the SC edges due to both the reduction of the SC power and the strong effect of laser peak power fluctuations on the supercontinuum bandwidth [START_REF] Genier | Amplitude noise and coherence degradation of femtosecond supercontinuum generation in all-normal-dispersion fibers[END_REF]. Furthermore, the pink curve in Fig. 6 depicts the estimated minimum RIN values measured by our system as a function of wavelength. Indeed, our scope has an estimated dark noise standard deviation of 1 mV. Thus, if we divide this number by the mean of our Gaussian distribution, for each spectrally resolved bin, we obtain a minimum estimated spectrally resolved RIN value of 0.14 %. Finally, we observe that our measured spectrally RIN values are all slighty above the estimated minimum.

Mininum RIN

Fig. 6: Measured spectrally-resolved RIN using the oscilloscope method with silicon (red dots) and InGaAs (green dots) photodetector and corresponding average normalized SC spectrum (solid blue), pumping at 45 • from the principal axes plus simulated SC (black dots) and RIN spectra (dashed brown). The pink curve shows the estimated minimum RIN values from our measurement system. These experimental measurements have been further compared to numerical simulations of the two-coupled generalized nonlinear Schrödinger equations (CGNLSEs), as those described in Ref. [START_REF] Ghosh | Chalcogenide-glass polarization-maintaining photonic crystal fiber for mid-infrared supercontinuum generation[END_REF] for highly-birefringent PM fibers. We performed SC and RIN simulations solving the CGNLSEs using the split-step method and the initial conditions as described in [START_REF] Genier | Amplitude noise and coherence degradation of femtosecond supercontinuum generation in all-normal-dispersion fibers[END_REF] with an input pulse duration of 180 fs (full width half maximum), an input peak power of 48 kW, an input amplitude noise value of 0.22 % (amplitude noise ≈ RIN/2) and associated anti-correlated pulse duration fluctuations to keep the energy of the mode-locked laser constant. To plot the simulated RIN, we used the following equation with an ensemble size of 50 independent pulses, which is sufficient to minimise modelling error:

RIN( ) = ˜ ( ) 2 -˜ ( ) 2 2
/ ˜ ( )

2

.

(

) 1
The simulation results are shown in Fig. 6 (black dots for the spectrum and dashed brown for the RIN). As can be seen, the agreement is excellent with the experimental spectrum, both for the shape and bandwidth, with a highly modulated spectrum and the generation of an XPMI sideband near 1400 nm. This agreement is also rather good between the measured and simulated RIN. Table 1 compares the noise performances of our SC source to those of the literature [START_REF] Gonzalo | Polarization noise places severe constraints on coherence of all-normal dispersion femtosecond supercontinuum generation[END_REF][START_REF] Rao | Ultra-low-noise supercontinuum generation with a flat near-zero normal dispersion fiber[END_REF]. The second column indicates the SC bandwidth. The third one explained the method used to measure the RIN of the SC system, either based on an oscilloscope in the time domain or on an Electric spectrum analyser (ESA) in the frequency domain. The fourth and fifth ones stipulate the measurement and filter bandwidth to obtain the average RIN specified in the last column. Few methods have indeed been implemented to measure the coherence or the RIN of ANDi SC generation [START_REF] Gonzalo | Polarization noise places severe constraints on coherence of all-normal dispersion femtosecond supercontinuum generation[END_REF][START_REF] Rao | Ultra-low-noise supercontinuum generation with a flat near-zero normal dispersion fiber[END_REF][START_REF] Tarnowski | Compact all-fiber source of coherent linearly polarized octave-spanning supercontinuum based on normal dispersion silica fiber[END_REF]. For the coherence, in [START_REF] Tarnowski | Compact all-fiber source of coherent linearly polarized octave-spanning supercontinuum based on normal dispersion silica fiber[END_REF], the authors used a dispersive Fourier transform technique to highlight the stability of ANDi SC generation. Then, using a fiber-based unequal-path Michelson interferometer, they demonstrate a high fringe visibility meaning a high degree of coherence for the ANDi SC generation. Also, in Refs. [START_REF] Gonzalo | Polarization noise places severe constraints on coherence of all-normal dispersion femtosecond supercontinuum generation[END_REF] and [START_REF] Rao | Ultra-low-noise supercontinuum generation with a flat near-zero normal dispersion fiber[END_REF], the authors used a setup composed of bandpass filters (10/12 nm), a fast photodiode and an oscilloscope to measure the spectrally-resolved RIN of the SC spectrum. In [START_REF] Gonzalo | Polarization noise places severe constraints on coherence of all-normal dispersion femtosecond supercontinuum generation[END_REF], the authors reported an average RIN of 33 %. This value can be considered high but, in this study, the coherence of ANDi SC generation was degraded by the PMI effect. In [START_REF] Rao | Ultra-low-noise supercontinuum generation with a flat near-zero normal dispersion fiber[END_REF], the authors reported a very low average RIN of 1.2 %, due to the low input peak power used to obtain their ANDI SC and the large bandwidth of their filters.

From Table 1, we see that we obtain the lowest average RIN value using the oscilloscope measurement method. About the value reported by [START_REF] Heidt | Low noise all-fiber amplification of a coherent supercontinuum at 2 m and its limits imposed by polarization noise[END_REF], it is not possible to compare directly our RIN value to the 0.07 % announced because they don't use any passband filters and thus measure a single RIN value for the whole SC spectrum. Furthermore, they limit their ESA measurement bandwidth from 10 Hz to 20 MHz while the repetition rate of their laser is up to 80 MHz. Therefore, it reduces the overall RIN value. Finally, we underline the setup composed by the Origami 10 HP laser and the PM-ANDi silica PCF is the system with the lowest average RIN to our knowledge and the highest PER (measured PER of 10 dB in [START_REF] Heidt | Low noise all-fiber amplification of a coherent supercontinuum at 2 m and its limits imposed by polarization noise[END_REF]), and thus it can be also considered as one of the state-of-the-art systems.

In conclusion, using a compact ytterbium mode-locked laser producing 180 fs pulses and a highly-nonlinear polarization-maintaining ANDi PCF, we were able to generate an ultra-flat octave-spanning SC spectrum with a bandwidth up to 720 nm (from 670 nm to 1390 nm) and associated PSD higher than 0.4 mW/nm. We have shown that the SC is linearly polarized with a PER of 17 dB across the whole spectrum bandwidth and with a spectrally-resolved RIN below 0.54 % from 700 nm to 1100 nm, which is, to the best of our knowledge, the lowest average SC ANDi RIN value ever reported with this measurement method. This system could find many applications in OCT, metrology, or as a wavelength-tunable fs laser source thanks to its octave-spanning bandwidth, flatness, low-noise, high coherence, polarization and near-linear chirp properties.

Fig. 1 :

 1 Fig. 1: Schematic of the setup, including ytterbium femtosecond mode-locked laser (ORIGAMI 10 HP), variable neutral density filter (ND), half wave plate (HWP), aspheric lenses (AL), 3D translation stages (S), 2 m of PM-ANDi PCF, flip-mirror (FM), integrating sphere (IS), 2 m of multimode pick-up fiber (MM), optical spectrum analyzer (OSA), monochromator (MONO), photodiode (PD), and oscilloscope (OSCI).

Fig. 5 :

 5 Fig. 5: Histogram of the variation of the energy of 10000 pulses measured with the oscilloscope. The red curve represents a Gaussian fit of the distribution.

Table 1 :

 1 Average RIN values in ANDi SC generation from the current literature and our system.

	Ref	Δ	Method Δ	Δ	RIN
	[5]	620 nm, -10 dB	Osci.	600 nm	10 nm	33 %
	[21] 480 nm, -30 dB	Osci.	200 nm	12 nm	1.2 %
	[11] 300 nm, -3 dB,	ESA	300 nm		0.07 %
	Ours 720 nm, -3 dB	Osci.	400 nm	3 nm	0.54 %

Funding. Horizon 2020 Framework Programme (722380); Agence Nationale de la Recherche (ANR) (ANR-15-IDEX-0003, ANR-17-EURE-0002).

Disclosures. The authors declare no conflicts of interest.