C Bordreuil

A Niel

Modelling of hot cracking in welding with a cellular automaton combined with an intergranular fluid flow model

Keywords: Hot cracking, Welding, Microstructure Prediction, Grain boundaries, Fluid network, Solidification

Hot cracking phenomena in welds are still a major issue to choose suitable welding conditions. Solidification cracking in welds depends on microstructure morphology, thermal and mechanical loadings. It is still difficult to discriminate the influence of the different phenomena because they are all connected. To investigate solidification cracking, numerical modelling that integrates most of these phenomena is developed in this paper.

Thermal and mechanical fields are computed with a finite element model at the process scale. The thermal field is combined with a cellular automata (CA) to predict microstructure at a meso scale. Based, on this prediction, the intergranular network is extracted.

Finally, the pressure in the intergranular fluid network is computed during solidification.

The fluid flow model along grain boundaries predicts an eventual nucleation of pore/void in the remaining liquid by investigating the maximum pressure drop in the solidifying microstructure. The numerical model is able to describe some hot cracking phenomena observed during experimental tests.

Introduction

Ensuring crystal cohesion during solidification in welding by avoiding cracks is a major challenge in material science and engineering. Factors affecting weld metal solidification cracking can be considered as either metallurgical or mechanical and often Email address: cyril.bordreuil@univ-montp2.fr (C.Bordreuil, A.Niel) June 11, 2013 interact together. The weld solidification cracking of alloys often appears in the mushy zone. This zone is a transition region between the weld pool and the base material. The behaviour of this region is complex because it mixes liquid and solid phases that interact under solidification. Different regions depending on the solid fraction, can be defined at the microstructure scale [START_REF] Dantzig | Soldification[END_REF]. The first region correspond to a high liquid fraction around dendrites where liquid can easily flow around solid skeleton. The second region is for high solid fraction but still with a continuous liquid network around dendrites or along grain boundaries. The last region is at a very high solid fraction with isolated liquid films and solid bridges are created between grains. In general, the coherency solid fraction f coh is defined around the transition between the two last regions. It is known from casting that defects as hot tearing and porosities can be due to solidification shrinkage and/or to tensile strain rate. Solidification shrinkage and plastic strain rate opening has to be compensated by fluid flow. If less flow is available, a depression can be created and a pore can nucleate. One important point that is rarely taken into account is the morphology of the intergranular fluid network and the flow along it. By using a microstructure prediction, the model is able to take into account this feeding along grain boundaries.

In welding, the hot cracking phenomenon is difficult to predict due to high gradient arising in the vicinity of the weld pool. These gradients load the solidifying microstructure. In particular, tensile stresses can be generated by differential thermal shrinkage upon cooling or by difference of plastic behaviour [START_REF] Chihoski | Understanding weld cracking in aluminium sheet[END_REF]. It is also often emphasized that grain morphology has an influence on hot cracking [START_REF] Wang | Solidification cracking of superalloy single and bicrystals[END_REF]. Due to fast cooling rate induced by the process, the microstructure morphology is heterogeneous at micro, meso and macroscale. At the microscopic scale, the grain demonstrates a dendritic substructure [START_REF] Prangnell | Effect of welding parameters on the solidification microstructure of autogeneous TIG welds in an Al-Cu-Mg-Man alloy[END_REF]. At the weld bead scale, the process generates a morphology with columnar and equiaxed grains. There is a transition from a columnar to equiaxed dendritic morphology called CET (columnar equiaxed transition). For the 6061 Al Alloy, the crack onset was observed in the columnar zone, far from the the weld pool in region of high solid fraction [START_REF] Niel | Modeling hot cracking phenomena in 6061 Aluminum alloy weld with a microstructural based criteria[END_REF] and at grain boundaries. The solidification shrinkage and/or tensile stresses can create a depression in the interdendritic spacing or at grain boundaries. Tensile stresses perpendicular to the grain boundaries tend to create an opening at grain boundaries.

The literature on hot cracking phenomena in welding is abundant [START_REF]Hot cracking phenomena in welds[END_REF][START_REF]Hot cracking phenomena in welds[END_REF][START_REF]Lippold and Bollinghaus and Cross[END_REF]. Chihoski [START_REF] Chihoski | Understanding weld cracking in aluminium sheet[END_REF] demonstrates, first, the influence of structural effects induced by the localized heat input. Matsuda [START_REF] Matsuda | Dynamic observation of solidification and solidification cracking during welding with optical microscope (i) -solidification front and behavior of cracking[END_REF] introduced the Brittle Temperature Range (BTR) based on observation. The BTR defined a temperature interval where the microstructure is in a critical configuration. In this model, strain rate and temperature are used to predict crack onset and doesn't introduce any metallurgical considerations (or variables). On superalloy, Wang [START_REF] Wang | Solidification cracking of superalloy single and bicrystals[END_REF] adapted a criterion of Rappaz [START_REF] Rappaz | M "A new hot tearing criterion[END_REF] for hot tearing to the case of hot cracking in welding of single and bi crystals. The Rappaz -Drezet -Gremaud [START_REF] Rappaz | M "A new hot tearing criterion[END_REF] (noted RDG in the rest of the text) criterion introduced some microstructural variables as well as macroscopic variables. This model can take into account the influence of lack of liquid feeding as well as some plastic strain localisation in the formulation. The variables are linked with the theory of multiphasic media and was able to predict crack onset in some configuration of welding [START_REF] Wang | Solidification cracking of superalloy single and bicrystals[END_REF]. Because, welding induces characteristic length (grain size, weld pool width) of the same order of magnitude, hot cracking susceptibility modelling needs an integrated approach to take into account process, microstructural morphology and liquid-solid phases under solidification. Due to the difficulty to measure mechanical and thermal fields in the vicinity of the moving weld pool, computational modelling must be performed to approximate stress/strain distribution in the mushy zone and to understand the influence of the different mechanisms on the hot cracking [START_REF] Safari | Hot cracking in stainless steel 310s, numerical study and experimental verification[END_REF]. The modelling needs to integrate process simulation, microstructural prediction and depression computation in the intergranular fluid network.

The purpose of this paper is to propose a methodology to investigate the couplings that can coexist between process and microstructure. It describes the modelling to approximate pressure drop in the continuous intergranular network. The computational model that combines thermo-mechanical, cellular automaton and pressure computation during the solidification is explained. The process simulation is done with a finite element thermo-mechanical model. The microstructural prediction is based on a cellular automaton based on solidification mechanisms driven by the thermal field. The pressure in the liquid network is then computed by adapting the RDG model to the microstructure predicted by the cellular automaton algorithm and loaded with thermal and mechanical fields. Some comments are also given on the influence of the coalescence of grain boundaries. Then, some results given by this approach are shown. The methodology is depicted with some two dimensional simulations done on 6061 Aluminium Alloy that was compared with some experiments [START_REF] Niel | Modeling hot cracking phenomena in 6061 Aluminum alloy weld with a microstructural based criteria[END_REF]. The relationship between microstructural prediction and fluid flow at grain boundaries is new an try to contribute to the explanation of hot cracking in welding.

Computational model

Different phenomena, acting at different scales, are involved in weld solidification. To model this phenomenon, different scales must be integrated:

• process scale : thermal and mechanical fields (strain rate).

• microstructure : fluid pressure computation around the solid skeleton in the intergranular fluid network.

To determine the temperature and stress/strain fields in the mushy zone, a thermomechanical finite element model is performed. This model allows to estimate fields during welding process. These fields are determined at the (macroscopic) process scale and drive the microstructure generation. To predict cracking, the microstructure morphology has to be determined. A cellular automaton (CA) is implemented to directly track grain morphology. The CA can take into account different physical mechanisms as competitive grain growth, grain selection and CET. In welding, the microstructure is assumed to be driven by the process and the coupling between microstructure evolution and thermal field can be considered as weak [START_REF] Rappaz | Prediction of grain structures in various solidification processes[END_REF]. In this model, micro and macro segregation are not taken into account. The CA results give grains boundaries but also crystallographic orientation and solid fraction can be estimated with thermal fields. Once microstructure is known, the pressure at grain boundaries can be determined by adapting the finite element RDG model developed in [START_REF] Niel | Modeling hot cracking phenomena in 6061 Aluminum alloy weld with a microstructural based criteria[END_REF].

The integrated algortihm considering different modelling is organized as follows:

1. Thermo-mechanical finite element model of the process to determine macroscopic fields.

2. Microstructural prediction with CA.

3. Conversion of CA data structure to multiphasic finite element to predict pressure in the intergranular network.

Then by comparing the pressure the intergranular network to cavitation pressure in as in RDG criterion [START_REF] Rappaz | M "A new hot tearing criterion[END_REF], the crack susceptibility can be determined. In this model, the important effect of hydrogen content [START_REF] Atwood | A three-phase model of hydrogen pore formation during the equiaced dendritic solidification of aluminium silicon alloys[END_REF] can be taken into account by a dependence of the cavitation pressure to the hydrogen content as explained in [START_REF] Coniglio | Mechanisms for solidification crack initiation and growth in aluminum welding[END_REF].

In the next sections, the different modeling are briefly discussed.

Thermo-mechanical finite element model

The model is conducted to simulate the behaviour of thin plate for full penetrated weld pool. For the model under consideration, it is not necessary to solve temperature and stress/strain fields simultaneously due to the weak couplings between the two problems. An uncoupled thermal and thermo-mechanical analysis is adopted. The thermal transient analysis is first performed followed by a quasi static mechanical analysis loaded by the thermal history. In this model, the thermal analysis is solved with an enthalpy formulation and no convection is taken into account for simplicity. Nevertheless, the strategy could be adapted if an analysis similar to Farzadi et al [START_REF] Farzadi | Phase Field simulation of weld solidification microstructure in an Al Cu alloy[END_REF] is performed by using a thermo convective model to compute thermal fields.

Heat transfer model

A two dimensional heat plane analysis is performed. The dimension of the plate is 400x100x3 mm. The material is 6061 Aluminium Alloy. Only half of the plate is discretized due to the assumption of symetry along the centre line. Elements are quadrangle isoparametric elements in the middle and outer part of the plate. A transition region is meshed with triangular elements. The arc heat input to the workpiece is modelled as a flux with a gaussian distribution with r the distance from the arc center:

q(r) = q max exp(-3( (x -x 0 ) 2 + (y -y 0 -V t) 2 r 0 ) 2 ) (1)
The q max = ηU I Table 1: Parameters used in the thermal transient analysis computed with the integration of the solid fraction in the solidification interval. Because the coupling between process and metallurgy wants to be analysed, the solid fraction is considered not dependent of the solidification path in the whole domain. This approximation is strong but ease the analysis of the coupling (see section 4). The enthalpy is tabulated with respect to the temperature. All the parameters are summarized in table 1. More details on the thermal transient analysis can be found in [START_REF] Niel | Modeling hot cracking phenomena in 6061 Aluminum alloy weld with a microstructural based criteria[END_REF]16].

Mechanical model

The geometry and the mesh are the same as the thermal analysis. To estimate the mechanical fields, the important ingredients are:

• The mechanical behaviour of the material.

• The deformation in the molten weld pool.

• The solidification shrinkage.

First, the mechanical behaviour of the 6061 alloy was deeply investigated under a wide range of temperature by Maisonnette [START_REF] Maisonnette | D "Influences mcaniques et mtallurgiques de procds hautes tempratures sur un alliage 6061[END_REF]. A mixed kinematic-hardening model was chosen and restauration is taken into account with two internal variables representing the state of the microstructure. One variable represents the initial microstructure and the other represents the microstructure with dissolved precipitates. In the mushy zone, the volumetric plastic strain rate [START_REF] Zavanglios | Modeling of the mechanical behavior of semi-solid metallic alloys at high volume fractions of solid[END_REF] is not taken into account. It is questionable but it seems that structural effects plays a more important role on the plastic strain state in the mushy zone than its behaviour.

Secondly, the weld pool can hardly sustain mechanical load which can be modelled as a stress free region. Previous plastic strain of the resolidified weld metal has to be Property Ambient 500 annealed. These two phenomena are modelled by defining a coherent temperature above which stresses and plastic strains are set to zero.

Third, the solidification shrinkage has a significant impact on stresses and strains.

Because, thermal stresses of an alloy is concerned, the effect of solidification shrinkage is modelled as thermal contraction in solidification temperature range. Some mechanical parameters are summarized in table 2.

More details on the mechanical models can be found in [16].

Thermal and mechanical fields

When the analysis are finished, the fields are extracted from the results databases.

The temperature are known at nodes. Thermal fluxes and plastic strains are extrapolated from gauss points to nodes. All fields are kept at nodes with respect to time. Then, these fields could be interpolated to any subdomain. The method to interpolate these values on the microstructural domain (CA) is a natural neighbour interpolation [19]. This interpolation has the advantage to let the procedure find the nearest points to build the shape functions.

Cellular Automaton

The cellular automaton has been already used to predict microstructure during welding process in two dimension [START_REF] Rappaz | Prediction of grain structures in various solidification processes[END_REF][START_REF] Dress | Cellular Automata Modeling of Weld Solidification Structure[END_REF] and in three dimensions [START_REF] Ichikawa | Modelling of solidification and grain growth in steel welds[END_REF]. The temperature at microstructural scale is directly governed by the thermal transfer at the process scale. In welding process, the thermal solution can be uncoupled for the microstructural prediction. To investigate the interest to study an integrated modelling with microstructure prediction and process simulation, a two dimensional Cellular Automaton was implemented. The two dimensional thermal field coming from the process modelling is used to feed modelling of physical mechanisms present during solidification. The main physical phenomena are :

• Nucleation of new grain in the liquid or on a wetting grain.

• Growth of dendritic structure from existing grains or newly nucleated ones.

• Coalescence of remaining liquid films or isolated pocket.

The CA algorithm follows the development of [START_REF] Rappaz | Prediction of grain structures in various solidification processes[END_REF][START_REF] Ch | M "A coupled finite element-cellular automaton model for the prediction of dendritic grain structures in solidification process[END_REF]. Some important features are To determine CA cell that has melt, the temperature field is used. CA cell with its center at a temperature higher than the melting temperature has its state index set to MELT. Cell in the border of the weld pool has their state index set to PARTIALLY MELT. In 2D, the growing shape of a cell is a square with the diagonals aligned with the crystallographic direction. One of the main key point of CA is the propagation of the shape. The neighbour cell is captured when the parent cell has its shape that reaches the cell centre. The neighbouring cell in MELT state has its state index assigns to PARTIALLY MELT and the crystallographic orientation is set to the parent one. A new smaller shape is associated with the capture cell. The growth is stopped if a neighbouring cell is in a PARTIALLY MELT state. The shape grows according to speed of growth v g which depend on undercooling and species concentration.

Thermal fields are known at times t proc n ,t proc n+1 where n subscript stands for the time step of the process simulation. The process and metallurgical modelling haven't got the same characteristic times and must be synchronized. The time increment at the process scale is denoted ∆t. The CA algorithm time increment δt has to be limited in order to have effective prediction [START_REF] Ch | M "A coupled finite element-cellular automaton model for the prediction of dendritic grain structures in solidification process[END_REF]. Times for CA are denoted: t CA k , t CA k+1 . Fields needs to be interpolated from the process domain to the CA domain with the different time steps (figure 1). To know temperature for each time step of the CA algorithm, a temporal interpolation is performed between times t proc n-1 , t proc n and t proc n+1 using a parabolic scheme. The temperature is kept in a First In First Out queue.

T CA (t CA ) = a(t CA ) 2 + bt CA + c (2)
where a,b,c are constant. They are determined with the condition:

T CA (t proc i ) = T proc (t proc i
) with the three temperatures in the FIFO queue.

Then, the microstructure is predicted following the CA algorithm depicted on table 3 and reimplemented from [START_REF] Rappaz | Prediction of grain structures in various solidification processes[END_REF][START_REF] Ch | M "A coupled finite element-cellular automaton model for the prediction of dendritic grain structures in solidification process[END_REF]. The microstructure morphology is driven by several mechanisms that are described in welding by Kou [START_REF] Kou | Welding Metallurgy[END_REF].

Nucleation To take into account the nucleation of grains below the liquid temperature of the alloy, an athermal model is used [START_REF] Dantzig | Soldification[END_REF]:

dn d∆T = n max ∆T σ √ 2π 1 ∆T exp(- 1 2 ( ln∆T -ln∆T 0 ∆T σ ) 2 ) ( 3
)
where n is the grain density, ∆T the total undercooling, ∆T σ is the standard deviation and ∆T 0 is the mean undercooling.

Growth The growth takes place in this direction and the growth speed followed a Kurz, Giovanola et Trivedi model [START_REF] Kurz | Theory of microstructural development during rapid solidification[END_REF]:

v g = A∆T n (4)
where A and n depends on depend on species concentrations in the alloy. This law was identified from Ivantsov solution and can be identified for binary alloys. Recently, it was shown that this law can be adapted to multicomponent alloy [START_REF] Hunziker | On the formation of a centreline grain boundary during fusion welding[END_REF] by introducing local linearized multicomponent phase diagram. Normally, microsegregation in the intergranular fluid network needs to be integrated in the modelling. Difference of concentrations would lead to difference in v g .

Coalescence Rappaz et al [START_REF] Rappaz | Last-Stage Solidification of Alloys: Theoretical Model of Dendrite-Arm and Grain Coalescence[END_REF] introduced a model to take into account the fact that bridging of dendrites occurs later at grain boundaries for two grains with different crystalographic orientations. This model was used in welding process of single crystals by Wang et al [START_REF] Wang | Solidification cracking of superalloy single and bicrystals[END_REF] to investigate hot cracking. A bridging undercooling is defined for a pure metal as:

∆T b = T m L f γ gb -2γ sl δ (5)
where γ gb is the grain boundary surface energy, γ sl the solid liquid surface energy, L f the latent heat of fusion and T m the temperature of fusion and δ is the thickness of the remaining liguid and is of the order of nm. The grain boundary surface energy is given by the Read-Shockley law [START_REF] Wang | Solidification cracking of superalloy single and bicrystals[END_REF][START_REF] Rappaz | Last-Stage Solidification of Alloys: Theoretical Model of Dendrite-Arm and Grain Coalescence[END_REF].

RDG Finite Element Model

The grain substructure is mainly dendritic due to fast cooling rate. A typical weld microstructure is depicted in figure 2. The RDG model was adapted [START_REF] Niel | Modeling hot cracking phenomena in 6061 Aluminum alloy weld with a microstructural based criteria[END_REF] to take into account heterogenous fields. This model is then considered to model fluid flow along grain boundaries .The weak formulation is reminded here.

L B T KBdx{p} = L N T (β df s dT Ṫ + (1 + β)f s εp⊥ )dx (6)
N is a shape function and B is the derivative of the shape function and {p} is the vector of degree of freedom. Where L is the length of a grain boundary element, K its permeability, f s is the solid fraction, εp⊥ is the plastic strain rate perpendicular to the grain boundary. As given in [START_REF] Rappaz | Study and modeling of hot tearing formation[END_REF], the plastic strain rate is given by equation 7.

εp⊥ = K localization εp⊥ ( 7
)
where εp⊥ is the average plastic strain rate deduced from the continuous mechanical model and K localization is a localization factor. This factor can depend on several microscopic factors. This value was approximated to φ/δ where φ is the grain size and δ is proportional to the width of the grain boundary.

In the modelling εp⊥ is taken into account if it is positive and negative. Negative values has a benefit effect if solidification shrinkage is important. εp 33 is taken into account because if tensile strain rate is observed it has to be compensated by liquid feeding. θ is the orientation of the grain boundary in the macroscopic frame.

The model is able to predict pressure drop due to solidification shrinkage as well as strain localisation with the plastic strain perpendicular to the grain boundary. The right hand term of equation 6 correspond to source terms. The integration is performed by a gauss quadrature along line element representing the grain boundaries. The internal variables in the source terms (right hand side of equation 6) are interpolated from the thermal and mechanical process simulation (see section 2.1.3).

One of the important point of the model is the solid fraction evolution vs temperature.

In the model, the 6061 solid fraction evolution is taken from [START_REF] Easton | Effect of alloy composition on the dendrite arm spacing of multicomponent aluminium alloys[END_REF] where it was computed with the help thermo-calc database. The solid fraction of 6061 Alloy is then tabulated in the solidification range defined from 540 to 645 • C . The solid fraction was chosen only dependent on the temperature. In reality the solid fraction evolution will depend on segration occuring at micro scale. Microsegregation can be different from point to point.

Before refining the model by taking into account this phenomenon, this evolution of the solid fraction is chosen to be the same for the whole domain, to appreciate the reliability of the approach that couples microscale modelling and process simulation.

The other important ingredient is on the left hand side in the computation of the stiffness matrix with the approximation of the permeability. In such computation, the Carman-Kozeny is normally used to predict flow (equation 9).

K = λ 2 2 180 (1 -f s ) 3 f 2 s (9)
where λ 2 is the spacing in the direction perpendicular to the fluid flow. Here, it corresponds to the final width of primary arm-spacing and it was estimated to 10 µm with micrographies. Equation 9 is the permeability for two dimensionnal flow. In the case of solidifying microstructure, the grains are not two dimensionnal and the fluid can come from the third direction and ease the fluid flow. The permeability is often higher for a three dimensional microstructure.

Conversion from CA domain to RDG Finite elements

The RDG based on Darcy and mass-conservation equations must normally be solved in the mushy zone. Grain boundaries are detected by scanning all CA cell and by looking edges along which I C changes from cell to cell. The edge is kept in a list. This list will contains all segments that will be converted in one dimensional RDG finite element model [START_REF] Niel | Modeling hot cracking phenomena in 6061 Aluminum alloy weld with a microstructural based criteria[END_REF]. Not all edges are converted to a RDG Finite element model. To insert edges in the model, the solid fraction at the midpoint is computed. If the solid fraction is less than the coherent solid fraction, the edges is inserted in the model. At the end a continuous network is obtained where the different stiffness and right hand side can be assembled to obtain the system of equations on the whole domain:

K{p} = F ( 10
)
where {p} is the pressure vectors on the whole microstructure .

The system of equations is augmented by a set of boundary conditions. In this work, it is assumed that when the solid fraction is lower than 0.5 the pressure is imposed to the atmospheric pressure. Once the pressure is known along grain boundaries, the location with the highest pressure drop in the solidification interval can be determined.

Each segment of the grain boundaries assumed to have a continuous liquid films (under f coh ) has a 'activated' tag. Grain boundaries along grains with different orientations have still liquid films at solid fraction higher than f coh . To take into account this fact, coherent solid fraction is modified to take into account the coalescence phenomena. A new coherent solid fraction is defined based on the coalescence undercooling. The coherent temperature is defined as the temperature corresponding to the coherent solid fraction.

Then, the coalescence undercooling is substracted to finally obtain the new coherent solid fraction. Grain boundaries satisfying the criteria to have a solid fraction lower than the coherent solid fraction are not integrated in the model but has their tag 'activated' set to True [START_REF] Dantzig | Soldification[END_REF]. This can lead to isolated liquid films.

Results

The model is applied to an hot cracking test developed at the lab. The test consists in a simple autogeneous bead on 6061 Aluminium Alloy [START_REF] Niel | Modeling hot cracking phenomena in 6061 Aluminum alloy weld with a microstructural based criteria[END_REF]. The welding speed is around 500mm/min, the current is 200A and the arc voltage is around 10V. Thermal and mechanical results are investigated and details in [START_REF] Niel | Modeling hot cracking phenomena in 6061 Aluminum alloy weld with a microstructural based criteria[END_REF]. In [16], the thermal model was compared with thermocouples and weld bead size and gives good correlations. The overall configuration of the coupon is shown on figure 3. The first domain corresponds to the plate, it is the process domain. The second domain is the CA domain where all microstructure evolutions will be performed.

Microstructural prediction

In autogeneous TIG welding, the macroscopic grain structure is controlled by a combination of thermal conditions that prevails at the solid-liquid interface and crystal growth A (m/s/K n ) 0.01 n 2 ∆T σ (K)

.5 ∆T 0 (K) 5.

n max grain/mm 3 1. k 0.05 T m ( • C) 645 T s ( • C) 540 µ (P a.s -1)
1.e-3 f coh 0.92 and 0.96 rate which is directly related to welding speed. The conditions vary considerably depending on the position at the solid-liquid interface. At the trailing edge, solidification rate is high and favours formation of dendritic equiaxed grain that can be seen in the center of the weld pool. On either side of this central axial grains, columnar dendritic grains were found and grew towards the central axial region. The solidification of these grains is driven by the maximum thermal gradient or undercooling. The width of the equiaxed zone is around 4mm and columnar zone is around 2mm.

Metallurgical parameters used are given in table 4.

In figure 4, the predicted microstructure for the welding test is shown. the model reproduces the competition growth from the base material and the nucleation mechanism at the trailing edge. The size of the equiaxed grains is more or less reproduced as well the orientation of the columnar grain. The model is able to reproduce the columnar to equiaxed transition observed during the test [START_REF] Niel | Modeling hot cracking phenomena in 6061 Aluminum alloy weld with a microstructural based criteria[END_REF].

In figure 5, the black arrow is for the primary arm-dendrite direction in the columnar grain and the blue line is for the secondary arm-dendrite direction for the columnar grain.

For the equiaxed grain, the two directions are equivalent. Time analysis of microstructure generation reveals that grains are nucleated in front of columnar grains due to a high undercooling and it activates the nucleation phenomenon.

These grains have no parent grain and grows in the thermal field. The crystallographic orientation is assigned randomly and the grain type is set to equiaxed. There is then a competition with the growth of this grain and the growth of the columnar grain, explaining the columnar to equiaxed transition.

Remaining liquid drops

Remaining liquid films in tensile stresses region can be critical for hot tearing. A grain boundary element has its tag 'activated' when its solid fraction is lower that the coherency solid fraction. If grain boundaries have solid fraction higher than the coherency solid fraction it means that bridges was created between each grains leading to isolated liquid drops. With grain misorientation on each side of the grain boundary, a higher undercooling (bridging undercooling) is necessary to join both grain. This means that the In figure 6, the red large edges represented grain boundaries which have still remaining liquid films. Some grain boundaries are isolated from the intergranular fluid network.

These films can not be feed by liquid coming from the weld pool. These isolated liquid films are critical sites for crack onset.

Pressure drop

In figure 7, the result of pressure is obtained for the system of equations 10 with a coherent solid fraction of 0.92.

In figure 7, The highest pressure drop can be located and gives the critical site for crack onset. In this simulation, the highest pressure drop is in a grain in the last stage of solidification. This means that if pressure drop is sufficiently high to create a pore, a crack will certainly appears inside the columnar grain boundaries and it is located far from the weld pool. This indicates the importance of the microstructure and in particular the grain boundaries network to feed parts in the last stage of solidification with liquid. The order of magnitude is 8kPa of depression. It is a little bit larger than those found in the litterature but it is of the same order of magnitude. This prediction is in accordance with some observations during GTAW of 6061 aluminum alloy [START_REF] Niel | Modeling hot cracking phenomena in 6061 Aluminum alloy weld with a microstructural based criteria[END_REF].

Coherency solid fraction influence

The coherency solid fraction is not well defined. It can vary between 0.92 and 0.96. If the pressure computation is performed with the second solid fraction, a higher depression will be found. It means that more liquid films coexist with solid, more the succeptibility of crack onset is important. The modification of the coherent solid fraction also modifies the grain boundary incorporation in the model. If the the coherent solid fraction goes from 0.92 to 0.96, grain boundaries with lower permeability are introduced in the modelling and lead to higher pressure drop.

Comparison of figures 7 and 8 indicates a maximum pressure drop with the coherent solid fraction fixed at 0.96. The maximum pressure drop is not located at the same location due to grain boundaries activated or not in the model. By increasing the coherent solid fraction, the permeability is also decreased in equation 6 and plastic strain rate become higher and has an important effect on the strain localization phenomena.

Strain rate localization and solidification shrinkage pressure

To investigate strain rate localization influence on depression, the right hand side of equation 6 can be decomposed between a solidification shrinkage pressure ∆p sh and a strain rate (localization) pressure drop ∆p ε . Because the modeling leads to a linear system, the total pressure drop can be decomposed as:

∆p = ∆p sh + ∆p ε (11)
The pressure drop due to strain localization is mainly due to plastic strain rate. In equation 6, the right hand source term associated to ∆p ε is proportional to solid fraction and perpendicular plastic strain rate. In order to see the distribution of plastic strain rate perpendicular to the grain boundary is shown in figure 10.

The highest plastic strain rate is seen in the equiaxed grains region. It means that at constant solid fraction, the critical site due to strain localization will be in the equiaxed zone. The problem is a little bit more complicated due to the non homogeneous fields and due to structural effects. The pressure drop is the result of high strain rate localization and to the incapacity of the grain boundaries to let the metal flow and feed the grain boundary opening due to the strain localization. In other words, if there is a strong strain localization but that the fluid network can feed the opening between two grains.

The fluid flow compensate teh opening and no pore nucleates.

To appreciate the influence of strain localization to pressure drop the same simulation as the one of figure 7 is performed but with forcing the pressure drop due to shrinkage equal to zero. The result is shown in figure 10.

The highest pressure drop is located in columnar zone and it is not located where there is the highest strain rate.

Discussion

The pressure drop is influenced by the solidification shrinkage and by the plastic strain rate as well as the capacity of the fluid to flow along grain boundaries. The model is not intend to be predictive. The modeling wants to see if the coupling between process simulation, microstructure generation and fluid flow along the grain boundaries can give information on the critical location for crack onset. In this term, the model reproduced the critical site along columnar grain boundaries. The highest depression is around 10kPa and is higher that depression measured in the liquid metal probably due to rapid solidification effects that are not taken into account in the model.

It has to be outlined, that our model is nevertheless a poor representation of the morphology and of the real behaviour of the solidifying microstructure. Different points can be improved. First, The morphology of the intergranular network follows the CA grid and doesn't take into account any curvature of the grain boundaries. Secondly, the behaviour of the liquid network can be discontinuous whereas in reality the fluid network is often continuous [START_REF]Shear behavior of AA6061 aluminium in the semisolid state under isothermal and anisothermal conditions[END_REF]. If the media is three dimensional, the critical site can easily be fed in a continuous network. Third, the model is non representative due to the non-dependence of the solid fraction in function of species concentration due to micro or macro segregation and to the non dependence on the solidification path.

The model developed considers no flow between the grain boundaries and the adjacent grains. In some metal it can be a poor assumptions and some improvements in the model must be done by considering fluid flow inside the grain and a mass exchange between the grain boundary and the grain.

The important features of the model are the couplings between the thermal and mechanical modelling at the process scale, the microstructure prediction based on a cellular automaton and the pressure prediction inside the intergranular network. The article seems to demonstrate that this kind of analysis is necessary to improve the understanding of hot cracking in welding. The model can detect the location of the critical site in the columnar grains boundaries as observed in [START_REF] Niel | Modeling hot cracking phenomena in 6061 Aluminum alloy weld with a microstructural based criteria[END_REF]. If it is now not predictive it can be used to optimize the process parameters to decrease the temperature rate in the critical zone by adding a heat source for example.

Conclusion

A new approach was developed in this paper that combines process simulations, microstructural prediction and multiphasic model to determine pressure at grain boundaries. The model can be enriched by refining physical solidification models. The approach demonstrates the interest to combine different physics and scales in a model to predict hot cracking. Nevertheless, it has to keep in mind that the material behaviour under solidification is three dimensional. Some improvements can be done for the thermal transfer model by developing some thermo-convective analysis in the weld pool. In aluminium alloy, hydrogen content is known to have a critical importance on the appearance of pore

3 √ πr 2 0 9

 329 is a surface flux (W/m 2 ) related to process parameters arc voltage U and current I, V the welding speed and where η is the efficiency of the process and r 0 characterized the radius of the heat distribution.In the heat transfer, heat is lost from the surface in the form of convection and radiation on all surfaces. The thermal model is conducted considering an enthalpy formulation considering latent heat release in the solidification range of the alloy. The enthalpy is Liquid Density (kg/m 3 ) 2200 Solid Density (kg/m 3 ) 2700 Capacity (J/kg/K) 896 Conductivity (W/K/m) 200 Volumetric Latent Heat J/m3 10 Exchange coeffcient W/m 2 K 10.

  reminded. The CA domain is defined as a structured grid of cells. Each cell ν has a state index I ν and a crystallographic index I C . The state index allows to know if the cell is melt, partially melt or solid. The CA domain is chosen in a part of the macroscopic domain where the crack is assumed to occur. A reduced domain allows to accelerate computational to appreciate the interest of these kinds of modelling. The CA domain is initialized by defining a voronoi diagram on the domain. The orientation (an angle in 2D) of each voronoi cell is assigned randomly. The voronoi cell size is chosen in order to match the grain size of the base material. If CA cell center belongs to one voronoi cell, the crystallographic orientation of the voronoi cell is assigned to the CA cell.

Figure 1 :

 1 Figure 1: Synchronization of the different time steps used in the coupled model.

Figure 2 :

 2 Figure 2: Dendritic structures configuration after welding of 6061 Aluminium alloy in a equiaxed zone

Figure 3 :

 3 Figure 3: Domains induced in the modelling of the cracking test.

Figure 4 :

 4 Figure 4: Evolution of the microstructure during moving arc heat input.

Figure 5 :

 5 Figure 5: Two dimensionnal microstructural prediction. A grain is represented by a domain of same color. The microstructural observation reproduces Columnar to equiaxed transition (CET)

Figure 6 :

 6 Figure 6: Activated tag along grain boundaries. Liquid grain boundaries are depicted in red color and with large segments.

Figure 7 :

 7 Figure 7: Pressure drop in Pa in the intergranular network (coherent solid fraction set to 0.92). Other values in the model are given in table 4

Figure 8 :

 8 Figure 8: Pressure drop in Pa in the intergranular network for a coherent solid fraction set to 0.96. Others values of the model are given in table 4.

Figure 9 :

 9 Figure 9: Perpendicular plastic strain rate distribution in the grain boundaries. Only positive strain rate are reported.

Figure 10 :

 10 Figure 10: Pressure drop in the intergranular network due to plastic strain rate at grain boundaries for 0.92 coherent solid fraction with K localization =1000

Table 2 :

 2 Parameters used in the mechanical analysis

  For each cell in PARTIALLY MELT state compute the dendritic growth with the undercooling ∆T ii. For each cell in MELT state compute the nucleation with undercooling ∆T (e) if t CA = t proc jump to next process step.

1. For each process time step t proc n (a) Assign new temperature in the FIFO queue. (b) Interpolate temperature for current time step (c) Find CA cells up to liquidus temperature and assign I ν to MELT. (d) For each CA time step i.

Table 3 :

 3 CA algorithm used to predict the microstructure

Table 4 :

 4 Solidification parameters used in the microstructural prediction

and then cracks and must be integrated in the modelling.


Soma Datta

Nabendu Chaki
email: nabendu@ieee.org

Person Identification Technique using RGB Based Dental Images

Keywords: HSI color format, Wiener filtering, Opening, watershed, Snake

Dental signature captures information about teeth, including tooth contours, relative positions of neighboring teeth, and shapes of the dental work. However, this is complicated as dental features change with time. In this paper, we proposed a new, safe and low cost dental biometric technique based on RGB images. It uses three phases: image acquisition with noise removal, segmentation and feature extraction. The key issue that makes our approach distinct is that the features are extracted mainly from incisor teeth only. Thus the proposed solution is low cost besides being safe for human.

Introduction

Person identification based on human teeth is very popular in the world. Generally, dental records have been used to identify the victims of disasters like 9/11 terrorist attack, Asian tsunami [START_REF] Choorat | A single tooth segmentation using structural orientations and statistical textures[END_REF] etc. According to these facts there are many advantages of using dental biometrics. It is very possible to identify an unidentified (mutilated) body by comparing the post-mortem (PM) records against ante-mortem records. This will produce the closest match of multiple identities [START_REF] Chen | Dental biometrics: Alignment and matching of dental radiographs[END_REF]. Dental biometric system based on radiograph images are also able to identify the correct person from a large set of database where manual method fails. Dental biometrics is worked upon the radiograph images. Frequently taken x-ray is very injurious to our health. The radiation of used X-Ray changes our DNA structure that leads to cancer [START_REF] Siltanen | Statistical inversion for medical x-ray tomography with few radiographs: I. General theory[END_REF]. One question naturally comes to our mind that can we use dental biometric on living people? The answer is 'yes' off-course we can by taking the RGB based dental image. The other reasons for choosing dental biometrics are as follows 1.

It is very precise and correct. 2.

We don't need an expert to treat and compare the result.

3.

It is not easy to replicate the dental signature.

The primary contribution of this research work is to propose a three phase novel technique for RGB dental image based authentication and verification using biometrics. Our proposed method takes the RGB images of human teeth, extract the features [START_REF] Choorat | A single tooth segmentation using structural orientations and statistical textures[END_REF] like shape; contour; length etc. Depending on the features system will take the decision whether record is matched or not. This will be the extended version of this re-search work. There are so many good biometrics methods other than dental are exist. Among them retina scan is invasive in nature and it is injurious to our eyes. For finger print, iris detection [4, [START_REF] Ito | A fingerprint recognition algorithm using phase-based image matching for low-quality fingerprints[END_REF][START_REF] Jea | A minutia-based partial fingerprint recognition system[END_REF][START_REF] Javad | Biometric identification through hand geometry[END_REF] special type of device set-up is needed where as dental image can be taken by any kind of camera. Hence to identify a person at any place is possible by only carrying the camera and this application. However this system does not require any extra high cost device like retina scan requires Topcon retinal camera and its cost is $15,995. Hence it is very less cost device comparative to other biometrics devices. This method is safe to our health. For this reason dental signature based biometrics method is helpful. Here we extracted the features of five incisor teeth. Hence we can say that RGB based dental biometrics is advantageous with respect to speed, storage requirement and human health.

Related Review Work on Dental Biometrics

Before starting of this work, we have gone through many good dental biometrics methods. Most of these methods work upon the dental radiographs. Here some of them are explained shortly. Anil K. Jain and Hong Chen dedicated a concept of semiautomatic contour extraction method for shape extraction and pattern matching [START_REF] Jain | Matching of dental X-ray images for human identification[END_REF]. The main problems in their approach are if the image is too blurred, their algorithm will not work and slight angle deviation in the database image and incoming images may not be handled with this approach. Said et al. [START_REF] Shah | Automatic tooth segmentation using active contour without edges[END_REF] offered a mathematical morphology approach, which used a series of morphology filtering operations to improve the segmentation. Morphological filters like top hat and bottom filters are used for tooth segmentation. Nomir and Adbel-Mottaleb introduced a fully automated segmentation technique [START_REF] Nomir | Human identification from dental X-ray images based on theshape and appearance of the teeth[END_REF]. It starts by applying iterative thresholding followed by adaptive thresholding to segment the teeth from both the background and the bone areas. After adaptive thresholding, horizontal integral projection followed by vertical integral projection is applied to separate each individual tooth. And this method can achieve the position of each tooth precisely. All of the above method works on the radiograph dental images where as our proposed method works on dental RGB image.

Proposed Method

This research work describes up to the dental biometrics features extraction procedure. Fig 1 shows the entire system. Image acquisition and filtering is the first step. For de-noising Winer filter is used. In the first phase of this method conventional methods are used that are described in the next sub sections in short. In the segmentation phase only the tooth region are segmented based on morphological operations. Segmentation of individual tooth is done by using the 'watershed' and 'snake' algo-rithm. In the last step the features like volume, size, area etc are extracted from the incisors teeth.

Image Acquisition and Filtering

We worked upon freely available database named 'Labial Teeth and Gingiva Image Database' [START_REF]Labial Teeth and Gingiva Image Database[END_REF]. Here dental images are different for resolution, file format and image scene type. Due to easiness of our research work we transform each image into jpeg format. All the rest of the works are based on these jpeg images. The qualities of the images that are obtained from the sensor camera are not so good, some noises exist. To remove these noises we used 'Winer' filter.

Segmentation

Segmentation phase contains two sub phases. The phases are 'Segment the tooth region based on morphological operations' and 'Segment individual teeth based on Water Shed and Snake algorithm'. The next subsections describe how the teeth are extracted finally.

Information Extraction

It is not so easy to extract only tooth region from the entire input image. The input image contains white space surrounded by the green color, lips, some part of face etc. These regions should be eliminated from the ROI. It is very challenging to separate the teeth and non-teeth region based on Red, Green and Blue component. Sometimes it fails to find the exact ROI due to the variation of white color that present in the teeth. Hence we have used HSI color model [START_REF] Gonzalez | Woods, digital image processing[END_REF]. The meaning of HSI color model is Hue (H), Saturation (S) and Intensity (I). The HSI color model described more exact color than the RGB color. Though HIS color model is non-uniform in perception, still it is one of the most popular color models for color image processing. The another advantage of using HSI color model is that it has human-intuitional advantages such as color image enhancement, fusion, skin area detection, segmentation, color based

      - <= = otherwise G B if H θ θ 0 360 (1)
Where

          - - + - - - + - = - ) )( ( ) )( ( )] ( ) [( 2 1 cos 1 B G B R B R G R B R G R θ � = 1 - � (�����) * [min (R, G, B)] (2)
I = (R+G+B)/3 (3)
Where, H stands for Hue i.e. pure color, S for saturation, i.e. the degree by which the pure color is diluted using white light and I for intensity i.e. Gray level.

Region of Interest (ROI) Extraction

After information extraction, the initial segmentation is done based on the hue and intensity component of the image. In this phase the hue and intensity component are multiplied pixel-wise. Now this is converted into binary image using self updating threshold method [START_REF] Liu | Document image binarization based on texture features[END_REF]. It stores the intensities of the pixels in an array. The threshold is calculated by using total mean and variance. Based on this threshold value each pixel is set to either 0 or 1. i.e. background or foreground. Thus here the change of image takes place only once. Suppose f(i,j) is the gray-scale value at pixel (i, j), and Ti is the segmentation threshold value at step i. To obtain a new threshold value, we have to threshold the original image using Ti to separate the image into teeth areas and non-teeth areas, where µ i 0 and µ i 1 are the mean gray values for the two areas.

(

) ∑ ∈ = pixel no j i f j i i _ / )) , ( ) , (( 0 µ (4) ( ) ∑ - ∈ = pixel no j i f background j i i _ / )) , ( ) , ( ( 1 µ (5)
The threshold value for step i+1 can be obtained as In Fig3.a the name plate portion of the image is mostly visible. This is not required at all, hence to remove this unwanted region we used the following procedure re-move_nameplate. Noises have been removed using erosion followed by dilation [START_REF] Gonzalez | Woods, digital image processing[END_REF] technique. Erosion erodes contours so small contours are removed but it also shrinks the desired contours. Dilation recovers the looses areas. Now the contours contains only tooth region. Algorithm3 is used to find the final ROI which is a subset of input image. Algorithm 3 accepts the saturation components (Sat_Mat) and binary mask (B_Mask) of tooth region which is cleaned by erosion and dilation techniques. The dilation of A by B, denoted, is defined as the set operation.

Procedure: Fainting_Nameplate Assumption Presence of the Nameplate as specified in [11] Input: hue and intensity image of the given teeth image, threshold_value obtain in Binarization algorithm

      ≠ = ⊕ φ A B z B A z  | ( 7
)
Where Ø is the empty set and B is the structuring element. In words, the dilation of A by B is the set consisting of all the structuring element origin locations where the reflected and translated B overlaps at least one element of A. The erosion of A by B, denoted, is defined as

} ) ( { φ = = Θ c z A B z B A  (8)
Step 1 temp= image(hue)*image(intensity) // Pixel wise multiplication of hue and saturation and store it into temp variable.

Step 2 Take the size of temp and store them in two different variable, [r, c] = size(temp).

Step

Start loop i=1 to r Start loop j=1 to c if(g(i,j)>=1) temp2(i,j)= threshold_value / Ti+1; End Loop End Loop Step4 End

Here, erosion of A by B is the set of all structuring element origin locations where no part of B overlaps the background of A.

Procedure: Find_ROI

Assumption: Nameplate should be removed from the binary image. Input: saturation matrix, binary matrix after applying Remove_Nameplate Output: exact teeth image In this step the individual teeth are segmented using mainly watershed [START_REF] Kass | Snakes:Active Contour Models[END_REF] and snake active contour model [START_REF] Roerdink | The Watershed Transform -Definition Algorithm and parallelization Strategies[END_REF] followed by some morphological preprocessing. In dental biometry it is not required to find the features of all the teeth. We have taken five incisors tooth for which features have been extracted. These teeth are shown in Fig 8 . After opening operation we have pointed out the incisors tooth.

Result

After applying procedure Feature_Extraction we got the following data for five teeth that mentioned in the table given below. This method can be applied on all of the remaining tooth. If we apply this algorithm to all 32 teeth then processing time will be much more higher than the proposed method and storage space requirement will be again coonsiderably high. However that case accuracy will be much higher. Similarly we can find the features of the tooth. Here for simplicity we have mentioned the features of incisors_central_1, incisors_central_2, incisors_central_3, inci-sors_central_1 and incisors_lateral (discussed in introduction).

Discussion and conclusion

In this research work we proposed a system by which the features of the incisors tooth have been extracted. We have worked upon Labial Teeth and Gingiva Image Database [20]. We have applied the algorithm upon 270 data from set 1 to set 6 in Labial Teeth and Gingiva Image Database. Our proposed method can correctly find out the features from these sample dental images. In the next phase of this research work these features will be matched with the database. To extract the tooth contour from the entire tooth region we have used watershed and snake algorithm. Sometimes it happens that these methods are unable to detect exact tooth contour due to the less 'gap-valley' area. Even if the proposed dental biometrics approach is a low cost solution, dental signatures do change with time. A dental signature change due to accidents, dental work and as human grows up. Hence biometrics method based on dental signature is not 100% accurate. Our future plan is to develop low cost, multimodal biometrics method based on dental signature and ear pattern which will be more accurate towards identifying a person.

Fig. 1 .

 1 Fig. 1. Block diagram of the system.

  Image Acquisition, Filtering Extracted features are shape, size, volume etc of incisors teeth Segment the tooth region based on morphological operations Segment individual teeth based on Water Shed and Snake algorithm Segmentation Phase object detection etc. The RGB to HSI color model conversation is done by using the following equations:-

  Fig 2.a shows the corresponding saturation image of the input image in fig 3. Saturation component contains a lot of information for both tooth and background region. So it is very difficult to extract only teeth region from background from saturation component. As a consequence hue and intensity components are used to segment ROI and after that saturation component is used for feature extraction.

Fig. 2 .

 2 Fig. 2. Input Image Fig 2.a Corresponding saturation Image

  image Output: Binary image of the given teeth image

Fig. 3 .

 3 Fig. 3. a. Output after using Binarization. Fig 3.b: After applying Fainting_Nameplate

  Output: Teeth image where Nameplate area is filled with threshold_value Step 1. Select an initial estimate of the threshold T. Step 2. Calculate the mean grey values �0 ��� �1 of the partition, R1, R2. Step 3. Partition the image into two groups, R1, R2, using the threshold T. Step 4. Select a new threshold as T= (�0 � �1)/2 \ Step 5.Repeat step 2-3 until the mean value �0, �1 values in successive iteration are equal. Step 6. End Again we used procedure Binarization upon Fig 3.b to get the image where teeth region and name plate both are present. Next task is to remove the name plate region. It is done by subtracting fig 3.a from fig 4.a, shown in fig 4.b.

Fig. 4 .

 4 Fig. 4. a After applying 2 times Binarization Fig 4.b Only teeth region with noise In fig 5.b some noise are present.Noises have been removed using erosion followed by dilation[START_REF] Gonzalez | Woods, digital image processing[END_REF] technique. Erosion erodes contours so small contours are removed but it also shrinks the desired contours. Dilation recovers the looses areas. Now the contours contains only tooth region. Algorithm3 is used to find the final ROI which is a subset of input image. Algorithm 3 accepts the saturation components (Sat_Mat) and binary mask (B_Mask) of tooth region which is cleaned by erosion and dilation techniques. The dilation of A by B, denoted, is defined as the set operation.

Fig. 5 .

 5 Fig. 5. The actual ROI after applying Find_ROI on fig 5.b

  Fig 5 is converted into its equivalent binary image using algorithm1 shown in Fig 6. Now we have segmented each contour using opening operation [12] shown in Fig 7.

Fig. 6 .

 6 Fig. 6.a. Binary image of Fig 5 Fig 6.b. After applying Opening on Fig 5

Fig. 7 .

 7 Fig. 7. After applying Watershed incisor are Fig. 8. Segmented contour after label matrix are segmented Technique. Now each contour has been identified by applying the label matrix technique [12] shown in Fig 8. The boundary of the contour is detected by using boundary detection algorithm as follows.

Fig. 9 .

 9 Fig. 9. incisors_central_1 Contour Fig. 10. incisors_central_2 ContourIn the last phase the features are extracted teeth wise. We have extracted features of the incisor teeth from right maxilla, left maxilla, right mandible and left mandible.The binary contour of each tooth (teeth_contour) is given as the input of the system. The algorithm returns the features of the teeth.

Table 1 .

 1 The calculated featured of individual teethRelative positions of CG points vary with the width of the teeth. CG points are helpful to determine the relative position.

	Teeth Name	Width	Height Volume	CG-Row CG-Col
	incisors_central_1	208	256	38100	240	427

Acknowledgement

Authors are thankful to Dr. Satyabrata Biswas for providing relevant knowledge and the "TEQIP Funded CoE Project", at Computer Science & Engineering Department, Calcutta University, for providing infrastructural facilities during the progress of the work.


Sagarika Chowdhury

Rajat Kumar Pal

Goutam Saha
email: dr_goutamsaha@yahoo.com

A Novel Double Fault Diagnosis and Detection Technique in Digital Microfluidic Biochips

Keywords: biochip, droplet, LOC, micro-fluidic technology, fluidic constraints

This paper presents a rigorous offline double fault diagnosis as well as a detection technique for Digital Microfluidic Biochips (DMFBs). Due to the underlying mixed technology biochips exhibit unique failure mechanisms and defects. Thus, offline and online test mechanisms are required to certify the dependability of the system. In this paper, the proposed algorithm detects double faults anywhere in the chip satisfying the dynamic fluidic constraints and improves the fault diagnosis time to an extent.

Introduction

An integrated microfluidic device incorporates many of the necessary components and functionality of a typical room-sized laboratory onto a small chip [START_REF] Erickson | Integrated Microfluidic Devices[END_REF]. These composite micro-systems, also known as lab-on-a-chip (LOC) (or bio-MEMS), offer a number of advantages over the conventional laboratory procedures and enable the handling of small amounts, e.g., micro-and nanolitres of fluids [START_REF] Chakrabarty | Design Automation Challenges for Microfluidics based Biochips[END_REF]. Droplet routing on the surface of the microfluidic biochip has been attracting much attention in recent years as it is one of the key issues to make use of the digital microfluidic device efficiently [START_REF] Zhang | A Droplet Routing for Fault-Tolerant Digital Microfluidic Devices[END_REF].

Microfluidic biochips have been characterized for the detection of faults [START_REF] Su | Testing of Droplet based Microfluidic Systems[END_REF], [START_REF] Su | Test Planning and Test Resource Optimization for Droplet based Microfluidic Systems[END_REF], [START_REF] Su | Testing and Diagnosis of Realistic Defects in Digital Microfluidic Biochip[END_REF], [START_REF] Srinivasan | A Digital Microfluidic Biosensor for Multianalyte Detection[END_REF]. The test planning problem was formulated in terms of Euler circuit in [START_REF] Su | Testing and Diagnosis of Realistic Defects in Digital Microfluidic Biochip[END_REF], [START_REF] Srinivasan | A Digital Microfluidic Biosensor for Multianalyte Detection[END_REF]. In [START_REF] Su | Defect Tolerance based on Graceful Degradation and Dynamic Reconfiguration for Digital Microfluidics based Biochips[END_REF], a functional testing has been proposed, referring to [START_REF] Su | Testing of Droplet based Microfluidic Systems[END_REF], [START_REF] Su | Ensuring the Operational Health of Droplet based Microelectrofluidic Biosensor Systems[END_REF], [START_REF] Su | Concurrent Testing of Droplet based Microfluidic Systems for Multiplexed Biomedical Assays[END_REF] that targets the functional operations of the microfluidic modules. Su et al. have proposed a defect tolerance based on graceful degradation and dynamic reconfiguration [START_REF] Xu | Functional Testing of Digital Microfluidic Biochip[END_REF]. A network flow based routing algorithm has been proposed in [START_REF] Kerkhoff | Testing of microelectronic-biofluidic systems[END_REF] for the droplet routing problem on biochips [START_REF] Yuh | BioRoute: A Network-Flow-Based Routing Algorithm for the Synthesis of Digital Microfluidic Biochips[END_REF]. An efficient diagnosis technique has been enhanced in [START_REF] Xu | Parallel Scan-Like Test and Multiple-Defect Diagnosis for Digital Microfluidic Biochips[END_REF] by Xu et al. such that multiple defect sites can be efficiently located using parallel scan-like testing. A more advanced multiple fault detection technique has been proposed in [START_REF] Chowdhury | Multiple Fault Detection Technique for Digital Microfluidic based Biochip[END_REF] by Chowdhury et al. in much less time compared to some previous techniques. However, it is not supported by the concept of Fluidic constraints, which has been explained in [START_REF] Roy | Modelling, Detection and Diagnosis of Multiple Faults in Cross Referencing DMFBs[END_REF].

This paper addresses the issue of double fault diagnosis and detection technique in digital microfluidic based biochips through a graph-theoretic formulation.

The rest of the paper is organized as follows: Section 2 presents the preliminaries of microfluidic arrays, their defect characterization, and the graph-theoretic formulation. Section 3 explains the proposed technique. Experimental results are reported in Section 4. Conclusions are drawn in Section 5.

Preliminaries

Structure of a Microfluidic Array

In digital microfluidic biochips, each droplet can be independently controlled by the electrodynamic forces generated by an electric field [START_REF] Yuh | BioRoute: A Network-Flow-Based Routing Algorithm for the Synthesis of Digital Microfluidic Biochips[END_REF]. Compared to the first generation biochips (analog), droplets can move anywhere in a 2D array to perform the desired chemical reactions, and electrodes can be re-planned for the different bioassays.

There are three key components in a biochip: 2D microfluidic array, dispensing ports/reservoirs, and optical detectors [START_REF] Yuh | BioRoute: A Network-Flow-Based Routing Algorithm for the Synthesis of Digital Microfluidic Biochips[END_REF]. The 2D microfluidic array contains a set of basic cells, which handle droplet movement. The dispensing ports/reservoirs handle droplet generation, and the optical detectors are used for reaction detection.

Defect Characterization

As has been described in [START_REF] Su | Testing of Droplet based Microfluidic Systems[END_REF], faults in digital microfluidic systems can be classified as being either catastrophic or parametric. Catastrophic faults cause a complete breakdown of the system while parametric faults degrade the system performance. To detect a fault, we need to pass a droplet across the cells so that it can traverse the whole path and reach towards the sink. If there is any defect within the microarray, the droplet gets stuck there. Otherwise, it reaches the sink at a predefined time.

With the use of electrowetting phenomenon, droplets can be moved to any location in the given 2D array. However, it should be taken care that each of the droplets is satisfying the observable fact of fluidic constraints described in [START_REF] Roy | Modelling, Detection and Diagnosis of Multiple Faults in Cross Referencing DMFBs[END_REF]. Only then the incidence of droplet collision is reduced.

Proposed Technique

Detection of a single fault in the biochip has reached at its end. However, diagnosis and detection of multiple faults are complex. Thus, let us start with more than one fault. Assume that, there are at most two faults in the biochip. The proposed technique takes a greedy approach to solving double fault detection technique. It tries to visit all the boundary nodes of Gmn during the first and second pass (that are P1 and P2) starting from source to sink. Remaining edges and nodes are traversed in subsequent passes P3, P4, and so on through certain movement patterns, which are explained below.

Movement Patterns

Column traversal from the source to the sink is based on one kind of movement expressed in Fig. 1 as Down-Left-Down-Right (DLDR). Let vi,j be the current node during a traversal.

See Fig. 1(a) to understand the movement Down-Left-Down-Right (DLDR) as follows, where the value of i (j) increases from top to bottom (left to right):

D L D R vi,j vi+1,j vi+1,j1 vi+2, j-1 vi+2, j
The sink can be reached from the source by the following appropriate sequence of movements mentioned above. Example 1 presents such a journey from source to sink. The source is at position (2, 1), and the sink is at [START_REF] Srinivasan | A Digital Microfluidic Biosensor for Multianalyte Detection[END_REF][START_REF] Srinivasan | A Digital Microfluidic Biosensor for Multianalyte Detection[END_REF]. The graph has 36 nodes and 60 edges. The journey from source to sink with the proposed pathway is indicated with solid arrows.

Fig. 2.

Four base nodes at a graph dual of G4×4.

Strategy

We locate four fixed nodes as the 'base' nodes from where the journey with the prescribed movement patterns begins. It is proposed that, irrespective of the source position, the traversal of the cells starts from any one of the nearest base point following the complete procedure.

3.3

The Complete Procedure

The complete process is as follows for a test array of size M×N.

Procedure Begin /* Assume source is at Base 1 and sink at (M, N) */

Step I: Boundary Test 1: A test droplet is dispensed from the source, and it traverses the boundary region clockwise and moves to the sink, as shown in Fig. 3(a).

Step II: Boundary Test 2: The second test droplet is dispensed with a delay of two units' time slice from the source and it traverses the boundary region anti-clockwise and moves to the sink, as shown in Error! Reference source not found.(b). If BT1 and BT2 fail, it ensures at most two faults are at the boundary region. Hence, go for Detect_Fault_UB algorithm and Detect_Fault_SB algorithm. Skip steps three and four. Otherwise, proceed to step three.

Step III: Row Test: Two iterations of parallel scan-like test with one row shift are carried out, having the time delay of two units during dispense of each droplet. The test movement pattern of row test is shown in Fig. 4. Thus, in iteration 1, all the even rows are traversed and in iteration 2 all the odd rows are covered. After the first iteration, there should be a delay of two units' time. After row test, the defected rows can be identified easily, following the arrival report of the droplets at the sink. Now to further reduce the testing area of the chip, let us identify the region of the affected rows following algorithm Reduce_Area (Fig. 5). Algorithm Reduce_Area:

1. If the affected rows are in the upper region, i.e. between row #1 to floor (M/2) of the chip, then divide it into two halves in such a way, so that, the divider passes just next to the one row, which is affected last (Fig. 5(a)). 2. Else if the affected rows are in the lower region, i.e. after floor (M/2) row number of the chip, then divide the chip into two halves in such a way, so that the divider passes just before the row, which is affected first (Fig. 5(b)). 3. Else divide the chip into two halves keeping faulty rows in both the halves (Fig. 5(c)), so that column test can be carried out in those two halves in parallel.

If the row test is satisfactory but BT1 and/or BT2 are not adequate, then consider Case 1.

Else if, row test is not agreeable, but BT1 and BT2 are satisfactory, then consider Case 2 and go to Step IV.

Else if row test is not acceptable and so for BT1 and/or BT2, then consider Case 3 and go to Step IV.

Else it ensures that the row test, BT1 and BT2 all are good enough. Thus, the chip is free from any fault.

Step IV: Column Test: Delay for two units. Repeat parallel scan-like test (for two iterations) for the columns, following DLDR movement patterns. Three units of time delay against dispersion of every droplet and four units of time delay during second iteration are required. These movements are explicitly shown in Fig. 6. End Procedure

Analysis and Detection of Faults

The total time required for the entire process is the sum of the required time for Step I,

Step II, Step III, and Step IV if they are carried out sequentially.

Step III and Step IV are to be performed if only they are required. As there can be double faults anywhere in the chip, including boundary, let us discuss their possible positions on the chip sequentially.

Case 1: Assuming two faults are anywhere at the boundary, i.e. the rest of the chip is fault free. Thus, the faults may be in the following locations:

Algorithm Detect_Fault_UB: /* Detect Fault(s) in Upper Boundary region*/

Assume the source at (1, 1) location and the sink at (M, N) location.

1. Disperse a droplet from source to the second column. 2. Go downwards by one row.

3. Go to the last column of that row and then to the sink. 4. After dispensing a droplet, a next droplet is dispersed after two units' delay. 5. Repeat Steps 1 through 4 incrementing the column number by one until a droplet fails to reach the sink in proper time. 6. When a droplet does not reach to the sink, then the interleaving column of the first row is detected as the defective cell.

Once an affected cell is found, do the following: 7. Send the next droplet just before to the affected cell. 8. Go downwards by one row. 9. Follow the right way and go upwards, so that the very next cell to the affected one at the boundary get touched. 10. Go to the last column of that row and then to sink. 11. After dispensing a droplet, a next droplet is dispersed after two units' delay. 12. Repeat Steps 7 through 11 as shown in Fig. 7, until a droplet succeeds to reach the sink in proper time.

Following the last droplet that may fail to reach the sink in time, we can straightforwardly identify the defective cell. It is justified that the intersecting points of row and side boundary are treated as the faulty cell, as we have assumed that there can be at most two faults in the chip. E.g., if the fault is at (M3, N) location, then more than one droplet set for row test must fail to reach the sink. In that case, for the last droplet, which fails to reach the sink, the intersecting point of row path and side boundary can be detected as faulty (Fig. 8).

Case 2: Boundary is passable. Thus, two faults are anywhere in the chip, other than the boundary.

Go for the Complete Procedure as discussed in Section 3.3. After performing column test, we can have several scenarios to be discussed further.

1. If two adjacent column tests fail to pass droplets against one affected row, then there must be one common intersecting point, which is affected definitely (Fig. 9(a)).

Sink Iteration 2 Iteration 1

Sink Source Source 2. Else, if two adjacent column tests fail to pass droplets, against two of the affected rows, then we are getting six considerable points among which two may be faulty (Fig. 9(b)). Thus, detection can be done afterward. Fig. 10 shows it clearly. 3. Else if one column test fails to pass droplet against one affected row, then two considerable points are there, of which one has not gone for column test. Thus, that nontraversed point is the faulty one (Fig. 9(c)). 4. Else if one column test fails, against two of the affected rows, then there are exactly two infected cells to be identified for sure (Fig. Case 3: One fault is at the boundary and the other is anywhere in the chip. Consider an example 2, where the graph dual of G7×7 depicts that it has an error at cell #2, i.e. on the upper boundary and at cell #25. Now, according the proposed method, BT1 will fail while BT2 will succeed to reach towards the sink in specified time. Therefore, we go for Row test and definitely, we will get row #4 as faulty. Hence, after performing algorithm Reduce_Area, it looks like the picture as shown in Fig. 11(a The detection time is compared with some existing technique. It shows that the proposed one is much superior compared to the other. Moreover, the proposed technique satisfies the dynamic fluidic constraints as well.

Experimental Results

Extensive offline testing has been done with a large number of arrays varying from 4×4 to 60×60 electrodes. Table I reports the details of the microarray. Table II reports the performance of the proposed technique, viz. the existing technique [START_REF] Xu | Parallel Scan-Like Test and Multiple-Defect Diagnosis for Digital Microfluidic Biochips[END_REF]. The performance of the proposed technique is divided into two columns: Proposed (min) and Proposed (max).

The proposed method not only diagnoses double faults in the chip, but it also detects the location. To test N×N target array, Boundary Test 1 and Boundary Test 2 are carried out first, and the methods take 2N+2 units' time. Next, row test is performed; in each iteration, this takes N units of time. Thus, up to this step, the proposed technique takes 4N units' time slice. After that, the N×N target array be partitioned into two halves, and the column test is performed in parallel. Here, if one has to go for column test, then it takes 2N/2, i.e. N units time slice in each iteration. Hence, as a whole it takes 6N units time slice. Therefore, the total fault diagnosis procedure includes 6N steps, i.e. O(N), and compared to parallel scan-like test and multiple defect diagnosis [START_REF] Xu | Parallel Scan-Like Test and Multiple-Defect Diagnosis for Digital Microfluidic Biochips[END_REF], which has 8N steps, the time needed for this diagnosis is reduced. Columns 6 and 7 of Table II show the percentage (%) improvement in each case. This is defined as: % improvement = ((Existing Time -Proposed Time) × 100) / Existing Time It can be seen that the proposed method has an improvement over time at about 12.5% minimum and 37.5% maximum. Table II

Conclusion

Efficient fault detection in a microfluidic biochip is an indispensable activity, as it is often used to operate in critical circumstances. In this paper, an advanced double fault detection technique has been proposed, which yields better results compared to some existing methods. Also, it satisfies the phenomenon of dynamic fluidic constraints. Thus, the chances of collision between two droplets are reduced. The proposed algorithm can detect any fault in the boundary region too. Therefore, further case studies can be done so that, not only double fault at the boundary, but more than two faults can be identified.

Fig. 1 .

 1 Fig. 1. (a) Movement patterns. (b) A path from source to sink following the proposed movement pattern of DLDR.Example 1: Fig.1(b) shows a graph dual G77 digital microfluidic array. The source is at position (2, 1), and the sink is at[START_REF] Srinivasan | A Digital Microfluidic Biosensor for Multianalyte Detection[END_REF][START_REF] Srinivasan | A Digital Microfluidic Biosensor for Multianalyte Detection[END_REF]. The graph has 36 nodes and 60 edges. The journey from source to sink with the proposed pathway is indicated with solid arrows.

1 ,

 1 For the graph dual Gm×n these base nodes are B1 B2 = v1,n, B3 = vm,1, and B4 = vm,n. These four base nodes are shown in Fig.2for Gm×n.

Fig. 3 .

 3 Fig. 3. (a) Boundary Test I, where droplet #1 and time slice T1. (b) Boundary Test II where droplet #2 and time slice T4.

 2 .Fig. 4 .

 24 Fig. 4. Row_Test. (a) Iteration 1: droplet #3 at time slice T7, droplet #4 at time slice T10, and droplet #5 at time slice T13. (b) Three droplets follow the last column to the sink. (c) Iteration 2: droplet #6 at time slice T18 and droplet #7 at time slice T21. (d) Two droplets follow the last column to the sink. 3. Follow algorithm Detect_Fault_UB. 4. Next, delay for two units of time and send another droplet, from second row towards the end of that row and then downwards to sink, following side boundary. If it reaches properly, no fault in the right side region. Otherwise, there is a fault. Then follow algorithm Detect_Fault_SB .

Fig. 5 .

 5 Fig. 5. Chip divider algorithm. (a) Affected rows are in the upper region. (b) Affected rows are in the lower region. (c) Affected rows are in any of the regions.

Fig. 6 .

 6 Fig. 6. Column_Test following DLDR movement patterns.

Fig. 7 .

 7 Fig. 7. Detection of double faults at the upper boundary of the chip. (a) First droplet fails to reach the sink, so the fault is at cell #2. (b) Second droplet fails to reach the sink, but cannot say whether the fault is on cell #3. (c) Third droplet fails to reach the sink, but cannot say whether the fault is on cell #4. (d) Fourth droplet succeeds to reach the sink. Thus, cell #4 is certainly faulty.

Fig. 8 .

 8 Fig. 8. Traversing through (a) to (d) while detecting the fault at the side boundary.

  ). Now, perform algorithm Detect_Fault_UB and it will detect a fault at cell #2. At this instant, droplet #D1 of Column Test at Iteration 2 fails to reach the sink. Thus, the intersecting point of row #4 and the path traversed by D1 during column test at iteration 2 is cell #25, as discussed in Case 2(b). Fig.11(b) shows this clearly.

Fig. 9 .

 9 Fig. 9. Fault analysis for Case 2, where indicates affected row, indicates column traversed in the first iteration, and indicates column traversed in the second iteration.

Fig. 10 .

 10 Fig. 10. Detection of double fault for Case 2(b).

Fig. 11 .

 11 Fig. 11. Double Fault Detection through Example 2 for Case 3.

Table 1 .

 1 shows the case-wise result, as discussed earlier. It shows that the proposed technique achieves improvement to a great extent. Test case details

	Sl. No.	Size	Source	Sink	Total Edge
	1	44	1, 1	4, 4	24
	2	55	1, 1	5, 5	40
	3	66	1, 1	6, 6	60
	4	1010	1, 1	10, 10	180
	5	2020	1, 1	20, 20	760
	6	2025	1, 1	20, 25	980
	7	3030	1, 1	30, 30	1740
	8	4040	1, 1	40, 40	3120
	9	5050	1, 1	50, 50	4900
	10	6060	1, 1	60, 60	7080

Table 2 .

 2 Performance of the proposed technique

				Time Slice Required			% improvement
	Sl. No.	Existing [14]	Case 1 Diag-Detec-	Proposed Case 2 Diag-Detection	Case 3 Diag-Detec-	in Diagnosis (Min) (Max)
			nosis	tion	nosis	(max)	nosis	tion
	1	32	20	24	24	32	28	28	12.5	37.5
	2	40	25	30	30	40	35	35	12.5	37.5
	3	48	30	36	36	48	42	42	12.5	37.5
	4	80	50	60	60	80	70	70	12.5	37.5
	5	160	100	120	120	160	140	140	12.5	37.5
	6	200	125	150	150	200	175	175	12.5	37.5
	7	240	150	180	180	240	210	210	12.5	37.5
	8	320	200	240	240	320	280	280	12.5	37.5
	9	400	250	300	300	400	350	350	12.5	37.5
	10	480	300	360	360	480	420	420	12.5	37.5