Aurore Back

Anaïs Crestetto

Ahmed Ratnani

Eric Sonnendrücker

An axisymmetric PIC code based on isogeometric analysis

Introduction

The Vlasov equation describes the evolution of charged particles in an electromagnetic field, which can consist of an applied and a self-consistent field. The latter being computed using Maxwell's equations. So we consider the system of the Vlasov-Maxwell equations. Physical problems described by this system are varied and it is necessary to develop methods adapted for each one of them.

We are interested here in the emission of electrons in a diode with hemispherical cathode. The problem is a priori three-dimensional and we can use cylindrical coordinates (z, r, θ). Moreover this problem is such that the unknowns do not depend on θ so that we can assume 2D axisymmetric geometry.

The IsoPIC code that we describe in this paper is a (z, r) axisymmetric Vlasov-Maxwell solver. Thanks to symmetry in θ direction, it is enough to give the description of the domain for a section (for example θ = 0). The isogeometric analysis approach will be applied, and we will construct a mapping F transforming a square into the section.

For the Maxwell equations, we consider in a first stage only the transverse electric mode (denoted by TE). The electric and magnetic field components considered are E r , E z and B θ . To solve them, we use spline finite elements. This article begins with their definition, their properties and their application to solve our system.

Since their introduction, B-splines have had a large success because their implementation is based on fast and stable algorithms. They are used in industry and in academic research for interpolation, computer aided design and data fitting. Recent work by Hughes and his co-authors [START_REF] Hughes | Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement[END_REF][START_REF] Bazilevs | Isogeometric Analysis, toward Integration of CAD and FEA[END_REF][START_REF] Hughes | Efficient quadrature for nurbs-based isogeometric analysis[END_REF] and the introduction of isogeometric analysis added another dimension to their use, creating an interface between simulation and modeling.

Before the recent work of Hughes [START_REF] Hughes | Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement[END_REF], the use of splines as basis functions in the finite element method was rare and mostly limited to uniform B-splines using periodic conditions, although the web-splines developed by Hoellig and others proposed a strategy for dealing with boundary conditions [START_REF] Höllig | Weighted extended b-spline approximation of dirichlet problems[END_REF][START_REF] Höllig | Finite Element Methods with B-Splines[END_REF]. The idea of isogeometric analysis, using geometric transformations based splines and NURBS is easier for most applications. Compared with usual finite elements, the main change due to the isogeometric analysis is undoubtedly the emergence of the k-refinement, a strategy that can increase the regularity of functions with the interfaces of the mesh to reduce the number of degrees of freedom.

Modern techniques of finite elements for Maxwell's equations are based on the ideas of differential geometry and more precisely the existence of discrete spaces that define an exact De Rham sequence. Following the pioneering ideas of Bossavit [START_REF] Bossavit | Mixed finite elements and the complex of Whitney forms[END_REF][START_REF] Bossavit | Computational electromagnetism and geometry : building a finite-dimensional "Maxwell's house" (part 5)[END_REF] a complete theory has been developed [START_REF] Hiptmair | Finite elements in computational electromagnetism[END_REF][START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF]. Buffa and al. [8,[START_REF] Buffa | Isogeometric analysis in electromagnetics: theory and testing[END_REF] have extended the isogeometric analysis to the steady state Maxwell's equations, providing a discrete exact De Rham sequence, in discrete spaces based on B-splines. The isogeometric analysis approach for studying unsteady Maxwell's equations has been studied by two of the authors [18].

The present article describes an axisymmetric Vlasov-Maxwell Particle In Cell (PIC) method based on isogeometric analysis. It is organized as follows: After recalling briefly the principles of isogeometric analysis, we introduce an axisymmetric Maxwell solver based on this approach. Then, we present the emission and the motion of particles in generalized coordinates on a patch. We study the conditions for extraction of particles in a diode with hemispherical cathode. And finally numerical results obtained with the code are proposed.

Isogeometric Analysis

In this section, we give a brief description of using B-splines in a finite element method. More details can be found in [START_REF] Hughes | Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement[END_REF][START_REF] Bazilevs | Isogeometric Analysis, toward Integration of CAD and FEA[END_REF].

B-Splines

Let T = (t i ) 1 i N +k be a non-decreasing sequence of knots. Definition 2.1 (B-Spline). The i-th B-spline of order k is defined by the recurrence relation:

N k j = w k j N k-1 j + (1 -w k j+1 )N k-1 j+1 , where w k j (x) = x -t j t j+k-1 -t j , N 1 j (x) = χ [t j ,t j+1 [ (x).
We note some important properties of a B-splines basis:

• B-splines are piecewise polynomial of degree p = k -1.

• Positivity.

• Compact support; the support of N k j is contained in [t j , .., t j+k ].
• Partition of unity:

N i=1 N k i (x) = 1, ∀x ∈ R.
• Local linear independence.

• If a knot t has a multiplicity m then the B-spline is C (p-m) at t.

Let (P i ) 1 i N ∈ R d be a sequence of control points, forming a control polygon.

Definition 2.2 (B-Spline curve). The B-spline curve in R d associated to T = (t i ) 1 i N +k and (P i ) 1 i N is defined by: Remark 2.3. In figure 1, we show an example of a B-spline curve. The basis here is simply Bernstein polynomials. In figure 2, we show the impact of inserting knots starting from a given knot vector. We remark the regularity of the new basis, for the repeated knot 3 4 . This refers as the k-refinement. The use of elements with high regularity at the interface of cells, allows us to reduce the dimension of the associated discrete space. More details on the impact of the k-refinement can be found in [18,[START_REF] Hughes | Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement[END_REF][START_REF] Beirao Da Veiga | Some estimates for h/p/k refinement in isogeometric analysis[END_REF].

M(t) = N i=1 N k i (t)P i .
Remark 2.4. In order to use B-splines to treat 2D or 3D problems, we need to use multivariate B-splines, which can be constructed by a tensor product. This will generate a mapping F that maps a rectangular grid Patch, onto the physical domain Ω. In figure 4, we show the impact of the h-refinement. Remark that using refinement strategies does not affect the mapping F.

Grid generation

Starting from a coarse grid, a new grid can be constructed by refinement, which can be done in three different ways:

• using the patch parameter h, by inserting new knots. This is the h-refinement, it is the equivalent of mesh refinement of the classical finite element method.

• using the degree p, by elevating the B-spline degree. This is the p-refinement, it is the equivalent of using higher finite element order in the classical FEM.

• using the regularity of B-splines, by increasing / decreasing the multiplicity of inserted knots. This is the k-refinement. This new strategy does not have an equivalent in the classical FEM. To define the computational mesh, we use alternatively h and p-refinement, see figure 3. The minimal degree of the basis functions is imposed by the domain design. When inserting knots, we can use uniformly-spaced knots or non uniformly-spaced ones.

Spline finite elements on general domains

In order to use the spline finite elements on general domains, we shall use the ideas of isogeometric analysis and consider domains defined with NURBS that can be obtained via CAD software and use these NURBS to map a rectangular domain on the physical domain.

Let Q be a cell in the physical domain. Q is the parametric associated cell and such that Q = F( Q). Let J F be the Jacobian of the transformation F, that maps any parametric domain point (ξ, η) into physical domain point (x, y). For any function v of (x, y), we associate its representation in the parametric domain

v((ξ, η)) = v • F((ξ, η)) = v((x, y)).
The basis functions R i will not be affected by these changes, the reader can always know if we are working in the physical or parametric domain thanks to (x, y) = F(ξ, η), x = α(ξ, η) and y = β(ξ, η).

Then,

α 1 = ∂α ∂ξ , α 2 = ∂α ∂η , β 1 = ∂β ∂ξ , β 2 = ∂β ∂η .
We have for the determinant of the Jacobian det(J F ) = α 1 β 2α 2 β 1 and for J F -1

J F = α 1 α 2 β 1 β 2 , J F -1 = 1 ∆ β 2 -α 2 -β 1 α 1 .
Let u be a (scalar or vector) function defined on the physical domain. When we use the patch coordinates, we will write it u, idem for the used spaces.

Boundary conditions

The boundary conditions are very easy to handle. Thanks to the mapping, the boundary of the physical domain is mapped into the boundary of the patch. Periodic boundary conditions can be handled by using periodic splines with an adapted knot vector [START_REF] Boor | A practical guide to splines[END_REF]. Perfect conductor and Silver-Müller boundary conditions are also implemented in a straighforward manner.

3 Axisymmetric Variational formulation of the 2D Maxwell equations

The 2D Maxwell equations, in the TE mode, read

           -∂E ∂t + curl B = J, ∂B ∂t + curl E = 0, div E = ρ.
For the axisymmetric geometry, we have:

E = E z E r , B = B θ , curl E = ∂ z E r -∂ r E z , div E = 1 r ∂ r (rE r ) + ∂ z E z and curl B θ = 1 r ∂ r (rB θ ) -∂ z B θ .
Variational formulation: As usual, to derive the variational formulation of the Maxwell equations, we will take the weak form of one equation (Ampère or Faraday), and keep the other one in the strong form. Multiplying each of the Maxwell equations by a test function and integrating on the physical domain Ω, we obtain for the Ampère equation

- ∂ ∂t Ω E • Ψ dX + Ω curl B • Ψ dX = Ω J • Ψ dX, ∀Ψ ∈ H(div, Ω).
And, for the Faraday equation,

∂ ∂t Ω B φ dX + Ω curl E φ dX = 0, ∀φ ∈ H 1 (Ω).
Now by using the Green formula

Ω (curl G) • F dX = Ω G curl F dX - Γ (G × n) • F dS , ∀ F ∈ H(curl, Ω), ∀ G ∈ H 1 (Ω), (1)
and if we are dealing with perfectly conducting boundary conditions, we get

∂ ∂t Ω B φ dX + Ω E • curl φ dX = 0.
In the case of axisymmetric geometry, the measure dX considered is simply dX = rdrdz.

Operators expressions: We have,

E • Ψ = E z Ψ z + E r Ψ r , J • Ψ = J z Ψ z + J r Ψ r , E • curl φ = E z 1 r ∂ r (rφ) -E r ∂ z φ, and curl B • Ψ = Ψ z 1 r ∂ r (rB) -Ψ r ∂ z B. Therefore, - ∂ ∂t Ω (E z Ψ z + E r Ψ r ) rdrdz + Ω (Ψ z ∂ r (rB) -rΨ r ∂ z B) drdz = Ω (J z Ψ z + J r Ψ r ) rdrdz,
and,

∂ ∂t Ω B φ rdrdz + Ω (E z ∂ r (rφ) -rE r ∂ z φ) drdz = 0.

Construction of the finite element spaces

An important feature of the functional spaces we chose for the variational formulation is that they form an exact sequence. Depending on the variational formulation we choose, we need to work with different exact sequences. As we have relaxed the Faraday equation, the following function spaces are involved

curl div H 1 (Ω) -→ H(div, Ω) -→ L 2 (Ω) ∪ ∪ ∪ V -→ W div -→ X (2)
In order to keep the specific features of Maxwell's equations at the discrete level we need to construct finite dimensional subspaces endowed with the same structure. The involved discrete spaces are denoted by V ⊂ H 1 (Ω), W div ⊂ H(div, Ω) and X ⊂ L 2 (Ω). In our case, we shall look for (E h , H h ) ∈ W div × V .

Spline finite elements on patch grids

We shall now start constructing the actual subspaces V ⊂ H 1 (Ω), W div ⊂ H(div, Ω) and X ⊂ L 2 (Ω). Our discrete spaces will be constructed using B-splines. The key point of our method is the use of the recursion formula for the derivatives:

N p i ′ (t) = p N p-1 i (t) t i+p -t i - N p-1 i+1 (t) t i+p+1 -t i+1 . (3)
It will be convenient to introduce the notation

D p i = p N p-1 i (t)
t i+p -t i . Then the recursion formula for derivatives (3) simply becomes

N p i ′ (t) = D p i (t) -D p i+1 (t). ( 4
)
Let us first consider a rectangular domain Ω. We consider the following discrete functional spaces:

V = span{N p i (x)N p j (y), 1 ≤ i ≤ N x , 1 ≤ j ≤ N y }, W div = span N p i (x)D p j (y) 0 , 0 D p i (x)N p j (y) , 1 ≤ i ≤ N x , 1 ≤ j ≤ N y } , X = span{D p i (x)D p j (y), 1 ≤ i ≤ N x , 1 ≤ j ≤ N y }.
As proved in the article by Buffa et al. [START_REF] Buffa | Isogeometric analysis in electromagnetics: B-splines approximation[END_REF], these spaces verify the same exact property as the spaces they approximate.

Let us now define coordinate changes conserving either the curl or the divergence of a vector field. Let us start from a vector field Ψ(ξ, η) = (Ψ (1) (ξ, η), Ψ (2) (ξ, η)) T defined on the parametric domain Q. Using the transformation formula:

Ψ = (-Ψ (2) • G ∂G 1 ∂y + Ψ (1) • G ∂G 2 ∂y ), (Ψ (2) • G ∂G 1 ∂x -Ψ (1) • G ∂G 2 ∂x ) T
for the vector fields of W div which conserves the divergence, using the diffeomorphism

G = F -1 , we get Ψ (1) = 1 ∆ (α 2 Ψ (2) + α 1 Ψ (1)
), and

Ψ (2) = 1 ∆ (β 2 Ψ (2) + β 1 Ψ (1)
) . So let Ψ = (Ψ (1) , Ψ (2) ) T be a function in W div , and Ψ = (Ψ (1) , Ψ (2) ) T be a function in W div . To preserve the divergence, the corresponding space of W div on the physical domain is

W div = {Ψ = ( 1 ∆ (α 1 Ψ (1) + α 2 Ψ (2) ), 1 ∆ (β 1 Ψ (1) + β 2 Ψ (2) )) T , Ψ ∈ W div },
which can also be written

W div = span 1 ∆ N p i (ξ) D p j (η) α 1 β 1 , 1 ∆ D p i (ξ) N p j (η) α 2 β 2 .

Discrete equations

Let W div = span{ -→ ψ b , b ∈ Γ E } the basis for W div , the discrete space associated to the electric field. Let V = span{ϕ b ′ , b ′ ∈ Γ H } the basis of V , it discretizes the space associated to the magnetic field.

The mass matrices for the electric and magnetic field are respectively

M W b,b ′ = Ω - → ψ b • - → ψ b ′ rdrdz, and M V b,b ′ = Ω ϕ b ϕ b ′ rdrdz.
The matrix for the curl is

K b,b ′ = Ω curl ϕ b • - → ψ b ′ rdrdz = Ω (r∂ r ϕ b + ϕ b )ψ z b ′ -(r∂ z ϕ b )ψ r b ′ drdz = Ω r(∂ r ϕ b ψ z b ′ -∂ z ϕ b ψ r b ′ )drdz + Ω ϕ b ψ z b ′ drdz.
This is the contribution of two terms. The first one "r(∂

r ϕ b ψ z b ′ -∂ z ϕ b ψ r b ′
)" which differs from the cartesian formulation by a multiplication by r. The second one "ϕ b ψ z b ′ " is due to the axisymmetric geometry.

Domain parametrization using Splines/NURBS curves

The parametrization of the domain by curves defined with aid of B-splines and NURBS is a recent topic. Indeed, it is about creating the mesh of the boundary of the domain from a description in terms of curves. Few articles that discuss the topic provide strategies to build these meshes [START_REF] Aigner | Swept volume parameterization for isogeometric analysis[END_REF][START_REF] Cohen | Analysis-aware modeling: Understanding quality considerations in modeling for isogeometric analysis[END_REF][START_REF] Martin | Volumetric parameterization and trivariate b-spline[END_REF]20].

In our case, the problem is different. We must be able to create a mesh from the description of the boundary and this mapping must be C 1 .

The easiest approach, which involves using NURBS to describe the quarter of circle of the diode, leads to a mapping which is continuous at the crossing of some nodes and not C 1 , see Fig. 5. We thus opted for a version in the spirit isoparametric analysis using only B-splines. The elements used to describe the boundary improve this description by increasing the degree of splines. A quarter of the circle is constructed by interpolation, by choosing a degree, determined by the user. Remark 3.1. In our numerical tests, we have token the weights equal to 1, and changed the knot vector T = {000 Remark 3.2. In figure 8, we give the description used to construct the domain, using NURBS.

A version of the domain, obtained by keeping the control points and weights equal to 1 (in which case the NURBS become B-splines), is presented later. This domain is used for validation, before looking a more advanced version.

PIC method for Vlasov equation

All the equations are solved on the patch. Since we have to change coordinates to map the physical domain into the patch, we must compute the Maxwell equations, the equations of motion, the densities of current and charge and the emission of particles in a general coordinates system.

The PIC Method

To solve the Vlasov equation

∂f ∂t + V • ∇ X f -(E + V ∧ B) • ∇ V f = 0,
we use the PIC method [19], [START_REF] Birdsall | Plasma Physics Via Computer Simulation[END_REF]. We consider a set of N macro-particles with X k being as the position, V k as the velocity, and ω k as the weight. The particles represent the distribution function f . We approach f by a sum of Dirac functions centered in the positions and in the velocities of particles:

f (X, V, t) ≈ f N (X, V, t) = N k=1 ω k δ(X -X k (t))δ(V -V k (t)). ( 5
)
The motion of particles is described by the equations of motion in which the electric and magnetic fields are used. They are obtained by solving the Maxwell equations which involves the charge density ρ and the current density J. We present the different steps of the PIC method:

Mesh Generation

Computation of fields on the mesh

Initialization of positions and velocities / /

Interpolation of charge and current densities on the mesh The most important steps that we develop in the following are:

-the description of equations of motion, -the computation of charge density and current density, -the condition for the emission of the particles.

The equations of motion

In our case, we move the particles in the patch, because it is easier to localize the particles on it. Let us explain how to obtain the equations of motion in generalized coordinates.

The equations of motion in Lagrangian mechanics are the Lagrange equations, also known as the Euler-Lagrange equations. There are available in any coordinate system:

d dt ∂L ∂ q (q, q, t) = ∂L ∂q (q, q, t),
where L(q, q, t) is the Lagrangian.

To write and solve the equations of motion in any coordinate system, we must know the Lagrangian to obtain the Lagrange equations. We start from its expression in cylindrical coordinates, where denoting by A the vector potential and φ the scalar potential, it has the form

L(z, r, θ, ż, ṙ, θ, t) = 1 2 m( ṙ2 + r 2 θ2 + ż2 ) -e(A r ṙ + A θ θ + A z ż -φ(z, r, θ)).
Now, we define the map F which transforms the cylindrical coordinates F (ξ, η, θ) = (z, r, θ) into a new coordinate system (ξ, η, θ). It is denoted by F (ξ, η, θ) = (z, r, θ). But, since we work in 2D axisymmetric geometry, we have θ = 0, and so the map is simplified by F (ξ, η) = (z, r). In this coordinate system, the Lagrangian has the form

L(ξ, η, θ, ξ, η, , θ, t) = 1 2 m(M ξ ξ2 + M η η2 + 2 M ξη ξ η) -e A ξ ξ + A η η -φ(ξ, η, θ) ,
with (A ξ , A η , A θ ) the components of the vector potential in the coordinate system (ξ, η, θ) and with

M ξ = ( ∂r ∂ξ ) 2 + ( ∂z ∂ξ ) 2 , M η = ( ∂r ∂η ) 2 + ( ∂z ∂η ) 2 , M ξη = ∂r ∂ξ
∂r ∂η + ∂z ∂ξ ∂z ∂η .

From this new Lagrangian, we can deduce the equations of motion in the new coordinate system:

det(J) d ξ dt + ξ2 K ξ,η + η2 K η,η + 2 η ξ K ηξ,η = - e m det(J) (((E + q ∧ B) | ξ ) M η -(E + q ∧ B) | η ) M ξη ) det(J) d η dt -ξ2 K ξ,ξ -η2 K η,ξ -2 η ξK ξη,ξ = - e m det(J) (((E + q ∧ B) | η ) M ξ -(E + q ∧ B) | ξ ) M ξη ) with K ξ,ξ = H ξ V ξ , K ξ,η = H ξ V η , K η,η = H η V η , K η,ξ = H η V ξ , K ξη,ξ = H ξη V ξ , K ξη,η = H ξη V η ,
where

H ξ = ∂ 2 z ∂ 2 ξ ∂ 2 r ∂ 2 ξ , H η = ∂ 2 z ∂ 2 η ∂ 2 r ∂ 2 η , H ξη = ∂ 2 z ∂ξ∂η ∂ 2 r ∂ξ∂η , V ξ = ∂r ∂ξ -∂z ∂ξ , V η = ∂r ∂η -∂z ∂η . det(J) = ∂z ∂ξ ∂r ∂η -∂z ∂η ∂r ∂ξ
is the Jacobian of the change of coordinates F , the components of E, B are the components of the electric field and the magnetic field in (ξ, η, θ) and q = ξ η .

Numerically, we consider ξ and η as independent variables with ξ = dξ dt and η = dη dt . We then solve the resulting first order system of ordinary differential equations using the second-order Runge-Kutta method.

Numerically, we can impose a negative electric field at the entry of the diode and allow it to propagate. If it is strong enough, it extracts electrons of the cathode, they go towards the anode. Our domain is meshed, so we look at each cell to see if the field satisfies the conditions of the emission of particles.

Extraction conditions

At a given frequency, particles are created in each cell Ω respecting the following conditions:

-the cell Ω touches the cathode, -the normal electric field to the cathode on Ω is greater than a threshold value.

Weights of created particles are positive and such that the following relation is respected: div E = ρ, and in the same time, the normal electric field to the cathode is zero according to the Child-Langmuir law. We integrate this relation on a cell Ω :

Ω div E dΓ = Ω ρ dΓ.
The Stokes formula leads to: Ω div E dΓ = ∂Ω E • n dγ. Besides, the approximation (5) gives:

Ω ρ dΓ = particle k∈Ω
-ω k because particles are electrons, their charge is negative but their weight is positive. We obtain the relation

∂Ω E • ndγ = ∂Ω 0 E • n dγ + ∂Ω 1 E • n dγ = particle k∈Ω -ω k ,
where ∂Ω 0 is the boundary of Ω touching the cathode and ∂Ω 1 is the union of other boundaries of Ω. We want to make ∂Ω 0 E • n dγ vanish, the weight of injected particles is imposed such that we have

∂Ω 1 E • n dγ = particle k∈Ω -ω k .
6 Numerical results

Validation of the Maxwell solver

Verification of the order of convergence

We test the order of convergence of the Maxwell solver on a square in axisymmetric geometry. Using the classical FEM discretization, the expected theoretical order for a refinement of order p is p + 1 for the magnetic field and p for the electric field. Remark 6.1. We have noticed that we do not get this order using H-div formulation. The order obtained is more of the form:

• the convergence order for B is :

p -1 2 , if p is even p + 1 2 , if p is odd
• the convergence order for E z is p + 1,

• the convergence order for E r is the minimum of the two previous ones.

However, using H-rot formulation we recover the classical order.

We solve numerically the Maxwell equations on a rectangle z ∈ (0, L) and r ∈ (0, R). Let m be the mode, the analytic solution is

H(z, r) = ω c λ p R ǫ 0 µ 0 J 1 ( λ p R r) cos( mπz L ) cos(ωt), E z (z, r) = 1 µ 0 ( λ p R ) 2 J 0 ( λ p R r) cos( mπz L ) sin(ωt), E r (z, r) = 1 µ 0 mπ L λ p R J 1 ( λ p R r) sin( mπz L ) sin(ωt),
where J n are Bessel functions of the first kind order n. λ p is the p th zero of J 0 . We plot the error in L2 norm between the computed solution and the analytic solution based on the maximum diameter of the mesh, denoted by h. For the tests we take [0, 1]×[0, 1] for the domain and N the number of cells equal to 4, 8, 16, 32 and 64. The time step is h × 10 -5 , small enough so that the order of the time scheme, which is 2 in our case, does not matter. We perform respectively 25, 50, 100, 200 and 400 iterations. The p-refinement is between degrees 0 and 3. We add the degree used for the domain, here it is equal to 1. The curves are represented with logarithmic scale in the figure 6.

Emission of particles in the diode

We studied qualitatively the extraction of particles in the domain representing the diode. We impose at the beginning of the diode a tension that increases linearly until reaching a threshold value, and is then constant. Particles are emitted at the cathode under the conditions of extraction described previously. We take the time step dt = 1.0 × 10 -2 and perform 25 000 iterations. We represent the electric field E r and the particles at the time T=25, 50, 75, 100, 125, 150, 175 and 200 in the figure 7.

Conclusion and perspectives

We have developed the IsoPIC code which solves the Vlasov-Maxwell equations in 2D axisymmetric geometry. It is based on isogeometric analysis. It is applied to emit electrons in a diode with a hemispherical cathode, which we approximated using splines. The solution of the Vlasov equation is performed with a PIC method by moving the particles emitted at the cathode.

The Maxwell solver in axisymmetric geometry gives good results. The theoretical orders of convergence are checked on a square and the fields evolve correctly in the diode. Qualitatively, we also have good results for the particles. The extraction seems to be good but we have to add a confinement for which the TM mode, we have not implemented up to now, is needed.

We noticed that the conjugate gradient method is not good for our problem because the mass matrix of the electric field in axisymmetric geometry is not well conditioned and this algorithm does not converge fast enough. We did not notice this problem in cartesian coordinates where the matrices are better conditioned. So we used the direct solver Pastix, which has the advantage to be faster than iterative solvers for this kind of problems in 2D.

The equations of motion are written and implemented in a general coordinate system. We have tested them in particular cases and we have compared the results with polar coordinates. The use of generalized coordinates is a choice for this study, despite the bigger calculation time. We would like to test later other possibilities, like cartesian coordinates with a mapping/inverse mapping to handle the mesh.

This work needs to be continued addressing the following topics. On the one hand, we can rewrite the Maxwell equations in order to not have to invert the mass matrix of the electric field. On the other hand, we have to compute a Poisson solver, in order to consider more general cases. In fact, if there are particles at initial time (which is not the case in this study), we have to solve the Poisson equation initially. Then, the conservation of charge is obtained thanks to the exact De Rham sequence. Finally, we can parallelize the Maxwell solver and the motion of particles using GPU (Graphics Processing Unit), to reduce computation time. This kind of parallelization is indeed very efficient in our costly computations, which do not require a lot of memory.

Figure 1 :

 1 Figure 1: (left) A B-spline curve and its control points, (right) B-splines functions used to draw the curve. N = 9, p = 2, T = {000, 1 4

Figure 2 :

 2 Figure 2: Using h-refinement with p = 2, T = {000, 111}, T = {000, 1 2 , 111} and T = {000, 1 2 , 3 4

Figure 3 :

 3 Figure 3: Grid generation. (left) after h-refinement, N 1 = 17, N 2 = 5, (right) after p-refinement, p 1 = p 2 = 3.

Figure 4 :

 4 Figure 4: Mapping from the patch to the physical domain: (left) initial patch, (right) patch after h-refinement in the η direction.

1

 1

3

  111} into T = {000 0.2 0.4 0.6 0.8 111}, in order to have a mapping which is C 1 .

Figure 5 :

 5 Figure 5: Diode: (left) the Jacobian of the mapping and the control points, (right) the mesh.

Figure 6 :

 6 Figure 6: L 2 norm for the electric fields E z , E r and the magnetic field B θ , for splines with degrees 1, 2, 3, 4.

Figure 7 :

 7 Figure 7: Electric field and electrons at time T=25, T=50 (first line), T=75, T=100 (second line), T=125, T=150 (third line), T=175 and T=200 (last line).

The authors would like to acknowledge financial support by CEA Gramat and thank Mr Vézinet for helpful discussions

The Dirac mass with a change of variables

The Dirac mass at a point x k is denoted by δ(xx k ). It has the following properties:

and g(x k ) = M δ(xx k )g(x)dx, for any continuous function g.

Let us consider the change of coordinates F (ξ, η) = x, M is a grid of the physical domain such as F ( M ) = M , and Jac(ξ, η) the jacobian of F . Since the integral value does not change when we change the coordinates we can deduce that

4.4 Computing J and ρ with a change of variables

Replacing f by its sum of Dirac function ( 5), we have:

Numerically, we have to compute the integral of these functions in space. It is easier to do it on the patch. With help of previous part and the equality (6), we deduce

and

where M is a cell of a physical domain and M a cell in the patch such as F ( M ) = M .

5 Particles emission

Short description of a diode

A diode is constituted of two semiconductors: a cathode rich in electrons and an anode which lacks them. If we apply a positive tension at the anode and a negative one at the cathode, such that the created potential drop is greater than a threshold value, electrons are extracted of the cathode and move towards the anode. This phenomena allows an electric current to pass. [18] Ahmed Ratnani and Eric Sonnendrücker. Arbitrary high-order spline finite element solver for the time domain maxwell equations. J. Sci. Comput. to appear.

[19] E. Sonnendrücker. Modèles cinétiques pour la fusion. notes du cours de M2 (2008).

[20] G. Xu, B. Mourrain, R. Duvigneau, and A. Galligo. Optimal analysis-aware parameterization of computational domain in isogeometric analysis. submitted.

1 st Knot vector {0 0 0 1


Alois Würger

Leidenforst gas ratchets driven by thermal creep

Keywords: numbers: 47.61.-k, 47.15.Rq, 44.20.+b, 68.03.-g numbers:

We show that thermal creep is at the origin of the recently discovered Leidenfrost ratchet, where liquid droplets float on a vapor layer along a heated saw-tooth surface and accelerate to velocities of up to 40 cm/s. As the active element, the asymmetric temperature profile at each ratchet summit rectifies the vapor flow in the boundary layer. This mechanism works at low Reynolds number and provides a novel tool for controlling gas flow at nanostructured surfaces.

Liquid spilled on a hot surface rapidly evaporates. At the Leidenfrost temperature well above the boiling point, however, one observes long-lived droplets that levitate due to the excess pressure resulting from the permanent feed of vapor at the bottom. Their contact-free suspension makes such Leidenfrost droplets very mobile. Linke et al. observed that, when placed on a millimeter-sized brass ratchet, the droplets rapidly accelerate to a speed of about 10 cm/s [1]. Very recently, Lagubeau et al. found that the same effect occurs for a piece of solid dry ice, and thus is not related to properties of the liquid phase [2]. Even more surprisingly, Ok et al. reported that reducing the ratchet profile to 200 nm, has little effect on the droplet velocity [3].

Unlike other self-propulsion mechanisms based on chemical or thermal gradients [4][5][START_REF] De Gennes | Capillarity and Wetting Phenomena[END_REF][START_REF] Chaudhury | [END_REF][8], this motion is not directed along an applied field but rather arises from the asymmetric surface structure of the solid support. This sawtooth profile acts as a rectifier transferring momentum on the interstitial vapor; the resulting gas flow advects the floating Leidenfrost droplet. The experimental findings [1][2][3] suggest that there is a common principle that works for both liquids and solids, and independently of the height of the ratchet profile. Several ideas have been put forward, relying on non-uniform Laplace pressure and Marangoni forces in the droplet, surface vibrations, or rectification of the radial vapor flow through the non-linear term of the Navier-Stokes equation [1][2][3]9]; yet none of them explains all of the mentioned experiments. In particular, the submicron ratchets of Ok et al. [3] exclude non-linear hydrodynamics as the dominant mechanism, as illustrated in the left panel of Fig. 1: For small profile D, the gas velocity and the effective Reynolds number Re in the ratchet layer are much smaller than at midheight where Re ∼ 1 [2]; thus rectification is expected to disappear for D ≪ h 0 , whereas the data of Ok et al. rather show the opposite behavior.

In his 1879 attempt to explain Crookes' radiometer experiment and building on Reynolds' theory for thermal transpiration [10], Maxwell showed the existence of thermal creep velocity along a solid-gas interface [11], where ν is the kinematic viscosity of the vapor and ∇ T the parallel component of the temperature gradient. Kinetic theory relates this gas flow to the non-uniform density and velocity distribution: Molecules coming from the cold side and hitting the surface at a given point, are more frequent but carry lower momentum than those from the hot side, thus resulting in an off-diagonal component of the surface stress and the boundary velocity v C [11]. Thermal creep drives aerosol thermophoresis [12], repels air-suspended particles from a hot surface [13], and operates in small-scale gas flow devices such as thermally actuated microcantilevers and Knudsen pumps [14][15][16][17][18][19].

v C = 3 4 ν ∇ T T , (1)
The present Letter points out the role of thermal creep for self-propelling Leidenfrost droplets and, in particular, analyzes the flow around a ratchet summit. The essential argument is illustrated in the right panel of Fig. 1. In the cleft below the droplet, there is a strong temperature gradient of several tens of Kelvin per micron. Because of the asymmetric profile, the horizontal component of the creep velocity has a finite mean value; the resulting gas flow along the ratchet surface drags the droplet toward the right. Note that this argument does not rely on the existence of the outward gas flow illustrated in the left panel of Fig. 1.

Our detailed analysis relies on Stokes hydrodynamics. In anology to thermal transport in colloidal dispersions [20,21], the droplet velocity is derived from the overall force balance on a closed surface. In the absence of external forces in horizontal direction one has

σ xn dS = 0, (2)
where σ xn is the stress pulling in x-direction on the area element dS with normal n. The stress tensor

σ ij = σ ′ ij - P δ ij comprises a viscous part σ ′ ij = η(∂ i v j + ∂ j v i )
and the excess pressure P . A non-uniform flow velocity v in the cleft of width h creates a stress of the order ηv/h, which by far exceeds the viscous drag at the remaining part of the droplet surface ∼ ηv/R. Thus the surface integral may be limited to the part between droplet and support; it closely follows the ratchet profile beyond the boundary layer, as indicated by the dashed line in Fig. 1.

The velocity profile in the cleft comprises two contributions of different origin. The first one, due to evaporation at the bottom of the droplet, is the outside gas flow in the left panel of Fig. 1; in the framework of Stokes hydrodynamics it does not contribute to the stress integral. Thus in the following, we consider the second velocity term, which arises from the thermal creep along the ratchet profile, as indicated by the arrows in the right panel.

The rectification mechanism is most obvious when comparing the viscous stress at the two slopes of the ratchet. The normal on the vertical part points in x direction; the corresponding diagonal element σ ′ xx = 2ηdv x /dx vanishes since v x and its derivative are zero. On the opposite side of slope m = D/L, the stress σ ′ xn is finite. The hydrostatic pressure varies little along the profile and will be discarded; then Eq. ( 2) reduces to the condition

σ ′ xz = 1 L L 0 σ ′ xz dx = 0 (3)
on the viscous drag on the slope of the ratchet tooth. If the droplet is immobile, the shear stress reads as σ ′ xn ≈ -ηv C /h, where h is the width of the cleft. In order to satisfy (3) the droplet moves at a velocity u, leading to 3) one readily obtains the expression for the drift velocity,

σ ′ xz = η[u -v C (x)]/h. Inserting σ ′ xz in (
u = v C /h 1/h . ( 4
)
The temperature profile is determined by the boundary conditions at the solid-gas interface, imposing continuous temperature and heat flow through the interface. Because of the important difference in thermal conductivity of the brass support and the vapor layer, κ V /κ S ∼ 10 -4 , the α FIG. 2: Temperature profile close to a solid-vapor interface. Due to the small conductivity ratio κV /κS, the isotherms (solid lines) are strongly distorted. At a distance of one molecular mean free path from the hot surface (dashed line), the temperature gradient in the vapor has a significant component parallel to the surface, ∇T , which is largest close to the upper corner.

temperature profile is strongly distorted, and the gradient is much larger in the vapor phase.

For a first estimate we calculate ∇T far from the corners, in the middle part of a ratchet tooth. In the limit κ V /κ S → 0 the brass surface is at constant temperature T S , and the profile in the cleft is given by T (x, z) = T B -(z/h)∆T , where ∆T = T S -T B and h = h 0 + xD/L. The resulting temperature gradient is perpendicular on the solid-vapor and droplet-vapor interfaces. The velocity distribution of the molecules hitting the brass surface is given by the temperature profile evaluated at one mean-free path from the brass surface, z = ℓh. At this finite distance, the gradient has a component parallel to the surface ∇T ∼ ∆T ℓ/h 2 . Replacing h with h 0 and discarding numerical factors gives a rough estimate for the drift velocity, u ∼ νℓ/h 2 0 .

(

) 5
With the mean-free path ℓ = 130 nm and the kinematic viscosity ν = 60 mm 2 /s of vapor at 300 • C, and h 0 ∼ 10µm [22], one finds u ∼ 10 cm/s, which is in qualitative agreement with experiment [1][2][3].

For the sake of a more quantitative description we refine the vapor temperature profile close to the upper corner of the ratchet, which turns out to dominate the creep flow. In analogy to the electrostatic potential of a charged polygon, a simple conformal transformation provides the expression T (r, ϕ) = T S -∆T (r/h 0 ) π/α sin(πϕ/α) [START_REF] Driscoll | Schwarz-Christoffel mapping[END_REF], where r, ϕ are polar coordinates with respect to the corner. The angle α is related to the aspect ratio m = D/L =cot α; for the ratchets of Ref. [3] one finds the exponent π/α ≈ 0.63. where ξ = -π 2 α 2 cos α 2 . Its essential feature is the weak singularity at the ratchet summit, very similar to the electric field close to a charged cusp. The molecular mean-free path ℓ provides a physical cut-off for the divergency. Now the drift velocity is evaluated in terms of Eq. ( 4), resulting in [23]

The resulting parallel component of the gradient along the dashed line close to point A reads [23]

∇T = ξ∆T ℓh -π/α 0 r 2-π/α , (6)
u = ξ ∆T T S ν h 0 ℓ h 0 π/α D/L ln(1 + D/h 0 ) (7)
with the numerical prefactor ξ ≈ 0.6 [23]. This expression confirms the estimate (5) yet shows additional dependencies on the ratchet parameters. Fig. 3 reveals a striking variation of u with D: the smaller the ratchet profile, the larger the droplet velocity. This at first sight counterintuitive result is confirmed by the experiment of Ok et al. [3]: Their data at intermediate temperatures are well fitted by a logarithmic variation, similar to [START_REF] Chaudhury | [END_REF]. Though this comparison does not account for the implicit dependence of h 0 on D, one may safely conclude on a qualitative agreement of [START_REF] Chaudhury | [END_REF] with the data. Note that for a ratchet driven by non-linear hydrodynamics, one expects the opposite behavior, i.e., a smaller velocity for small D. Indeed, from the left panel of Fig. 1 it is clear that for D ≪ h 0 , the gas velocity in the ratchet layer is small and the non-linear term (v • ∇)v of the Navier-Stokes equation insignificant. So far we have considered liquid droplets on a ratchet. The same mechanism holds for a piece of dry ice (solid CO 2 ) above its sublimation temperature floating above a hot metal surface; when graving a saw tooth profile at its lower face, Lagubeau et al. observed motion similar to the droplets discussed so far. Since the creep velocity occurs at the bottom of the dry ice, the mean velocity of hot u FIG. 4: Thermal creep below a piece of dry ice (solid CO2) floating above a hot metal surface due to sublimation [2]. Due to the ratchet profile printed at its lower side, there is a parallel temperature gradient ∇T as indicated by dashed arrows; the rectified thermal creep flow propels the disk to the left. The mean velocity of the vapor in the cleft is zero. the vapor in the cleft is zero, as illustrated in Fig. 4.

Lagubeau et al. measured the force F required to immobilize a droplet floating on a ratchet [2]. In the range R = 1...7 mm, they found values from 3 to 30 microNewton, and a power law F ∝ R β , with an exponent β ≈ 1.5. In the present work, this force is given by the integral of the shear stress over the contact area, F = πR 2 η v C /h . With the above expression for the thermal creep velocity one finds

F = πξ ηR 2 D L 2 ν ∆T T ℓ π/α h 1+π/α 0 . (8)
With the relation h 0 ∝ √ R [2], the force varies with the droplet size as β = 3 2 -π 2α ≈ 1.2; within the experimental uncertainities, this compares favorably with the measured value.

The thermal-creep mechanism described here is not limited to the motion of Leidenfrost droplets. As a straightforward application we discuss the gas pump shown in Fig. 5, which consists of two nanostructured plates at temperatures differing by ∆T . Both solid interfaces show thermal creep and thus impose a unifom gas flow across the cleft. For a sufficiently small ratchet profile, D < h 0 , Eq. ( 7) simplifies to

u 0 = ξ ν L ∆T T ℓ h 0 π/α . ( 9
)
This velocity may attain several meters per second. It turns out instructive to compare this ratchet with a Knudsen pump; in the present case, the thermal gradient is perpendicular to the gas flow, whereas both are parallel in the latter. Moreover, a Knudsen pump requires the system size to be comparable to or smaller than the mean-free path, and thus is restricted to very dilute gases. Thermal creep gives rise to a uniform gas flow at velocity u 0 , as given in Eq. ( 4). Note the opposite orientation of the saw teeth on the cold and hot side.

Although the ratchet mechanism does depend on the ratio ℓ/h 0 , it works for films that are hundred times thicker than the mean free path.

FIG. 1 :

 1 FIG. 1: Leidenfrost droplet on a ratchet. The left panel shows the gas flow due to the evaporation at the bottom of a droplet of radius R; the temperature of the solid TS is significantly above the droplet's boiling temperature TB. The arrows in the right panel indicate the thermal creep flow along the ratchet of period L and height D. The dashed line indicates the boundary layer; its thickness ℓ corresponds to the molecular mean free path. In recent experiments the ratchet height D varies from 200 nm to about 1 mm; for small D one has h0 = 10...100µm.

FIG. 3 :

 3 FIG.3: Drift velocity u as a function of the ratchet parameter D for h 0 = 10, 30, 100µm. The curves are calculated from Eq. (7) with ν = 60mm 2 /s, ℓ = 130 nm, π/α = 0.63, L/D = 4, and ∆T /T = 1/2.

cold hot u 0 FIG. 5 :

 05 FIG.5: Gas pump driven by thermal creep. The thermal gradient across the channel is given by their temperature difference ∆T and spacing h0. At the ratchet corners, there is a parallel component ∇T , as indicated by dashed arrows. Thermal creep gives rise to a uniform gas flow at velocity u 0 , as given in Eq. (4). Note the opposite orientation of the saw teeth on the cold and hot side.


Anoek Hj Maes

Paul Van

Maes Ahj

Van Diest Pj

Prof J Paul

MD Van Diest

The mucinous variant of columnar cell lesions

Keywords: columnar cell lesions, flat epithelial atypia, mucinous, core needle biopsy, breast columnar cell lesions * flat epithelial atypia * mucinous * core needle biopsy * breast

published or not. The documents may come

The mucinous variant of columnar cell lesions

Introduction

Mucinous lesions of the breast are rare. Although mucinous breast carcinomas have been recognized since 1816 [START_REF] Geschickter | Gelatinous mammary cancer[END_REF] , benign and premalignant mucinous lesions of the breast have only been described over the last few decades. A mucocele-like tumor was first described by Rosen in 1986 as multiple cysts with mucinous secretions lined by flattened or cuboid epithelium with cytologically bland cells, often showing extrusion of mucinous contents in the intervening stroma ("mucinous dissection") [START_REF] Rosen | Mucocele-like tumors of the breast[END_REF] . Other authors published variants of mucocelelike lesions in a morphological continuum with mucinous usual ductal hyperplasia (UDH), mucinous atypical ductal hyperplasia (ADH), mucinous papillary lesions, mucinous ductal carcinoma in situ (DCIS) and mucinous or colloid carcinoma [START_REF] Chinyama | Mammary mucinous lesions: congeners, prevalence and important pathological associations[END_REF][START_REF] Farshid | Mucocele-like lesions of the breast: a benign cause for indeterminate or suspicious mammographic microcalcifications[END_REF][START_REF] Ro | Mucocelelike tumor of the breast associated with atypical ductal hyperplasia or mucinous carcinoma. A clinicopathologic study of seven cases[END_REF][START_REF] Tan | Mucinous breast lesions: diagnostic challenges[END_REF][START_REF] Weaver | Mucinous lesions of the breast. A pathological continuum[END_REF] .

Columnar cell lesions (CCLs) have recently drawn much attention as precursor lesions within the low nuclear grade breast cancer spectrum [START_REF] Aulmann | Invasive Tubular Carcinoma of the Breast Frequently is Clonally Related to Flat Epithelial Atypia and Low-grade Ductal Carcinoma In Situ[END_REF][START_REF] Dabbs | Molecular alterations in columnar cell lesions of the breast[END_REF][START_REF] Dessauvagie | Characterization of columnar cell lesions of the breast: immunophenotypic analysis of columnar alteration of lobules with prominent apical snouts and secretions[END_REF][START_REF] Feeley | Columnar cell lesions of the breast[END_REF][START_REF] Schnitt | Columnar cell lesions of the breast[END_REF][START_REF] Simpson | Columnar cell lesions of the breast: the missing link in breast cancer progression? A morphological and molecular analysis[END_REF] that are frequently diagnosed in breast screening programs as they are often associated with microcalcifications on mammography [START_REF] Chinyama | Mammary mucinous lesions: congeners, prevalence and important pathological associations[END_REF][START_REF] Farshid | Mucocele-like lesions of the breast: a benign cause for indeterminate or suspicious mammographic microcalcifications[END_REF][START_REF] Kim | Benign and malignant mucocele-like tumors of the breast: mammographic and sonographic appearances[END_REF][START_REF] Leibman | Mucocelelike lesions of the breast: mammographic findings with pathologic correlation[END_REF][START_REF] Tanaka | Indeterminate calcification and clustered cystic lesions are strongly predictive of the presence of mucocele-like tumor of the breast: a report of six cases[END_REF] . Only two case reports have described a mucocele like tumor and columnar cell hyperplasia occurring in a morphologic continuum [START_REF] Coyne | Columnar cell hyperplasia with intraluminal crystalloids and features of a mucocoele-like lesion[END_REF][START_REF] Fadare | Mucocele-like tumor and columnar cell hyperplasia of the breast occurring in a morphologic continuum[END_REF] and one paper has reported on mucin extravasation and mucocele-like lesions accompanying CCLs [START_REF] Begum | Mucin extravasation in breast core biopsies--clinical significance and outcome correlation[END_REF] , suggesting that the mucinous variant of CCL may be very rare. However, to the best of our knowledge, we are the first to systematically describe the incidence of mucin production in CCLs. Moreover, to study the possible precursor role of mucinous CCLs in the mucinous pathway towards mucinous carcinoma, we reviewed mucinous carcinomas for the presence of mucinous CCLs in the vicinity of the invasive cancer and compared them to ductal carcinomas.

The breast always produces mucins, but most of them are membrane bound. In mucinous lesions, there is a switch from membrane bound mucin towards the production of secreted gel forming mucins [START_REF] Kim | Diversity of mucin genes, structure, function, and expression[END_REF][START_REF] Matsukita | Expression of mucins (MUC1, MUC2, MUC5AC and MUC6) in mucinous carcinoma of the breast: comparison with invasive ductal carcinoma[END_REF] . Mucins are high molecular-weight glycoproteins expressed by different types of epithelium and can be categorized into two groups: membrane bound (MUC1, MUC3A, MUC3B, MUC4, MUC12, MUC13, MUC17) and secreted or gel forming (MUC2, MUC5AC, MUC5B, MUC6) [START_REF] Diaz | Expression of epithelial mucins Muc1, Muc2, and Muc3 in ductal carcinoma in situ of the breast[END_REF][START_REF] Rakha | Expression of mucins (MUC1, MUC2, MUC3, MUC4, MUC5AC and MUC6) and their prognostic significance in human breast cancer[END_REF] . Of all MUC genes, MUC2 has been described to be most discriminating between mucinous and ductal invasive carcinomas and the most frequently expressed in mucinous carcinomas, playing a major role in the production of the extracellular mucin [START_REF] Matsukita | Expression of mucins (MUC1, MUC2, MUC5AC and MUC6) in mucinous carcinoma of the breast: comparison with invasive ductal carcinoma[END_REF][START_REF] Rakha | Expression of mucins (MUC1, MUC2, MUC3, MUC4, MUC5AC and MUC6) and their prognostic significance in human breast cancer[END_REF][START_REF] Adsay | Pathogenesis of colloid (pure mucinous) carcinoma of exocrine organs: Coupling of gel-forming mucin (MUC2) production with altered cell polarity and abnormal cell-stroma interaction may be the key factor in the morphogenesis and indolent behavior of colloid carcinoma in the breast and pancreas[END_REF][START_REF] O'connell | Altered mucin expression is a field change that accompanies mucinous (colloid) breast carcinoma histogenesis[END_REF][START_REF] Chu | Mucin expression in mucinous carcinoma and other invasive carcinomas of the breast[END_REF] altered mucin expression in these lesions precedes mucinous carcinoma [START_REF] O'connell | Altered mucin expression is a field change that accompanies mucinous (colloid) breast carcinoma histogenesis[END_REF] . Therefore, we also investigated MUC2 expression for our mucinous CCLs. The CNBs were performed under ultrasound or stereotactic guidance, using 14-gauge needles in the UMCU, RHT, and AHN. In the MHU, 16 to 18-gauge and sporadically 11-gauge (with Mammotome) were used. We reviewed the mammography records to ascertain the reason for taking the CNB: density, microcalcifications or both. In densities group, also palpable masses, architectural distortions and asymmetries were included.

Histology

All original hematoxylin-and-eosin-stained (H&E) slides from the CNBs were reviewed by two experienced observers (AM and PJvD) and scored for presence of CCLs and atypical ductal hyperplasia originating in a CCL (ADH-CCL). CCLs were classified into four categories according to Schnitt and Vincent-Salomon: columnar cell change (CCC), columnar cell hyperplasia (CCH), CCC with atypia and CCH with atypia [START_REF] Schnitt | Columnar cell lesions of the breast[END_REF] . CCLs are variably enlarged terminal duct lobular units (TDLU), lined by enlarged columnar cells with only one or two layers of columnar type epithelium in CCC. The cells contain elongated nuclei, with an orientation along the basement membrane and containing inconspicuous or absent nucleoli.

Apical snouts are often present, frequently showing intraluminal secretions and microcalcifications. CCH has a similar appearance as CCC, but the TDLUs are lined by more than two stratified cell layers. In CCC with atypia (CCC-A) and CCH with atypia (CCH-A), cytonuclear atypia is superimposed showing relatively round or ovoid (instead of elongated) not regularly orientated nuclei. The nuclei mostly have an irregular shape, often with prominent nucleoli. According to the WHO criteria [START_REF] Tavassoli | WHO Classification Tumours of the Breast and Female Genital Organs[END_REF][START_REF] Tavassoli | WHO Classification Tumours of the Breast and Female Genital Organs[END_REF] , lesions with enlarged ducts with complex architectural patterns with arcades, bridging or micropapillae were considered as ADH or low-grade DCIS, depending on the size of the lesion and the extent of the architectural complexity and regularity. ADH lesions arising in the context of a CCL (ADH- All CCLs were reviewed again to ascertain the presence of mucin in the lesion. For the definition mucinous CCL, we used the following criteria. At first, a columnar cell lesion has to be present as described above, comprising the largest part of the lesion. Secondly, mucin must be seen in the centre of the TDLU, often combined with extravasation of the mucin in the surrounding stroma (mucinous dissection) (Figure 1). A Periodic Acid-Schiff or Alcian blue staining can be used to confirm the presence of mucin. See Figure 1, 2 and 3 for examples of mucinous CCC, mucinous CCH and mucinous ADH-CCL. The lesion has to be differentiated from a mucocele-like tumour that has (one layer of) flattened or cuboid lining epithelium. Other lesions in the differential diagnosis are mucinous atypical ductal hyperplasia and mucinous DCIS, exhibiting both the characteristic architectural and cytological atypia according to the usual criteria [START_REF] Tavassoli | WHO Classification Tumours of the Breast and Female Genital Organs[END_REF] , in combination with mucin in the cysts or ducts [START_REF] Chinyama | Mammary mucinous lesions: congeners, prevalence and important pathological associations[END_REF][START_REF] Carder | Surgical excision is warranted following a core biopsy diagnosis of mucocoele-like lesion of the breast[END_REF][START_REF] Hamele-Bena | Mammary mucocele-like lesions. Benign and malignant[END_REF] .

Mucinous CCLs in mucinous versus ductal carcinomas cases

From the UMCU archive, we retrieved surgical resection specimens of 46 pure mucinous carcinoma cases diagnosed between 1990 and 2009 and 46 consecutive invasive ductal carcinoma cases from January 2009 onwards. The slides were reviewed (by AM and PJvD) and scored for the presence of CCL with or without intraluminal mucin.

Immunohistochemistry

Mucinous CCLs with sufficient material left in the CNB blocks were stained with the Bond-Max Autostainer, using the mouse monoclonal antibody MUC2 (Santa Cruz, H-300, 1:400).

Appropriate positive and negative controls were used throughout. MUC2 was scored for the (stained) intraluminal mucin and presence of cells with intracytoplasmic staining.

Electron microscopy

One mucinous columnar cell lesion was chosen for evaluation by electron microscopy according to usual procedures for processing paraffin material.

Results

A total number of 4164 breast CNBs were reviewed. Table 1 shows the numbers of CCLs and ADH-CCL and the abnormalities at mammography leading to these CNBs. As summarized in Table 2, 17 mucinous CCLs of which 3 with atypia were found in 17 women as well as 3 mucinous ADH-CCL in 2 women (overall frequency 20/4164= 0.5%). Mucinous dissection was present in 18 out of 20 (90%) of these lesions. All of them were associated with microcalcifications at mammography (one combined with a density), and these microcalcifications were mostly recognized in the CCL at histology as well. The mean age of these women was 53.9 years (range 36 to 74 years).

Luminal mucin was confirmed by Alcian blue staining (Figure 4), and PAS-D staining showed the presence of intracytoplasmic mucin as well (Figure 5). Electron microscopy of this part of the specimen confirmed the presence of many intracytoplasmic vacuoles with a sharp boundary (Figure 6). Assuming these are filled with mucin, this suggests that intracellular mucin in the columnar cells may be the origin of the intraluminal mucin.

Moreover, MUC2 was expressed in the intraluminal mucin in 12 of 15 (80%) mucinous CCLs (Table 3, Figure 7) and showed cytoplasmic expression in some cells in 5 of 15 (33%) mucinous CCLs seen in CNBs.

Mucinous CCL was significantly more frequent (13/46, 28%) in mucinous carcinomas cases than in ductal carcinoma cases (4/46, 9%) (p = 0.004) (Figure 8). Atypical CCLs (CCL-A or ADH-CCL, both mucinous and non-mucinous CCLs) were more frequent in mucinous invasive carcinoma cases (14/19, 74%) compared to ductal carcinomas (8/19, 42%) (p = 0.049).

Three women (two with mucinous CCH-A and one with mucinous CCH in the CNB) underwent a surgical excision biopsy between 1 and 7 months. In two of the three follow up surgical excisions, CCH-A without mucin was found, and in one a mucinous CCH. During the follow-up period (mean 3.1 years, range 0.5 to 7 years), no invasive carcinoma or DCIS occurred.

Next to the subsequent excision biopsies of the mucinous CCLs in the CNBs, 68 other excision biopsies were performed of the remaining 292 CCLs and ADH-CCLs. Of these 68, five showed a mucinous CCL in the excision (1 CCC, 2 CCH, 2 CCL-A) with in one also atypical lobular neoplasia. One subsequent invasive ductal carcinoma of intermediate grade with intermediate grade DCIS showed a mucinous CCL with atypia as well. In one case, invasive lobular carcinoma had preceded the mucinous CCL four years before.

Discussion

This is the first study that systematically describes the incidence of mucinous CCL in breast CNB. 17 mucinous CCLs of 291 CCLs (5.8%) and 3 mucinous ADH-CCL of 21 ADH-CCL (14.4%) were found, together concerning 0.5% of all breast CNB. Coyne et al and Fadare et al 17, 18 were the first ones who both described such lesions in single case reports before, suggesting a potential relationship giving their coexistence and the morphologic transitions between the lesions. Begum et al [START_REF] Begum | Mucin extravasation in breast core biopsies--clinical significance and outcome correlation[END_REF] has reported about the coexistence of CCLs with mucin extravasation and mucocele-like lesions as well. Besides, mucocele-like lesions have been the subject of several other studies, often presented in coexistence with other entities like mucinous usual ductal hyperplasia, ADH, DCIS through to mucinous carcinoma [START_REF] Chinyama | Mammary mucinous lesions: congeners, prevalence and important pathological associations[END_REF][START_REF] Ro | Mucocelelike tumor of the breast associated with atypical ductal hyperplasia or mucinous carcinoma. A clinicopathologic study of seven cases[END_REF][START_REF] Weaver | Mucinous lesions of the breast. A pathological continuum[END_REF][START_REF] Hamele-Bena | Mammary mucocele-like lesions. Benign and malignant[END_REF] . These lesions are postulated to represent a morphologic and biologic spectrum [START_REF] Ro | Mucocelelike tumor of the breast associated with atypical ductal hyperplasia or mucinous carcinoma. A clinicopathologic study of seven cases[END_REF][START_REF] Weaver | Mucinous lesions of the breast. A pathological continuum[END_REF][START_REF] Leibman | Mucocelelike lesions of the breast: mammographic findings with pathologic correlation[END_REF] . So possibly an evolutionary pathway could be present from benign mucinous lesions to invasive mucinous carcinomas, just like described by Abdel-Fatah for the tubular and lobular carcinomas on the basis of histological evidence and reported genetic changes [START_REF] Powe | Morphologic and molecular evolutionary pathways of low nuclear grade invasive breast cancers and their putative precursor lesions: further evidence to support the concept of low nuclear grade breast neoplasia family[END_REF] . In these pathways, CCLs (especially those with atypia) seem to play an important role, which can be hypothesized to be true for mucinous CCLs in the mucinous pathway as well. Indeed, in the study described by Begum et al [START_REF] Begum | Mucin extravasation in breast core biopsies--clinical significance and outcome correlation[END_REF] , DCIS was found in a few excision biopsies with CCL and mucin extravasation in the initial CNBs. We have not been able to demonstrate development of mucinous DCIS or invasive mucinous cancer (or other types of DCIS or invasive cancer) from mucinous CCL in our series. This may be related to the low number of mucinous CCL patients, the relatively short follow-up interval and the fact that only very few mucinous CCLs were treated by surgical excision biopsy.

We found significantly more mucinous CCLs in mucinous carcinoma cases than in ductal carcinoma cases, confirming that mucinous CCLs indeed might be important in the developmental pathway towards mucinous carcinoma. This is underlined by the fact that CCLs (both mucinous and non-mucinous) were more frequently atypical in mucinous than in ductal carcinoma cases.

MUC2 is typically expressed in the cytoplasm of tumour cells in mucinous carcinoma [START_REF] Matsukita | Expression of mucins (MUC1, MUC2, MUC5AC and MUC6) in mucinous carcinoma of the breast: comparison with invasive ductal carcinoma[END_REF][START_REF] Rakha | Expression of mucins (MUC1, MUC2, MUC3, MUC4, MUC5AC and MUC6) and their prognostic significance in human breast cancer[END_REF][START_REF] Adsay | Pathogenesis of colloid (pure mucinous) carcinoma of exocrine organs: Coupling of gel-forming mucin (MUC2) production with altered cell polarity and abnormal cell-stroma interaction may be the key factor in the morphogenesis and indolent behavior of colloid carcinoma in the breast and pancreas[END_REF][START_REF] O'connell | Altered mucin expression is a field change that accompanies mucinous (colloid) breast carcinoma histogenesis[END_REF][START_REF] Chu | Mucin expression in mucinous carcinoma and other invasive carcinomas of the breast[END_REF][START_REF] Yonezawa | Expression of mucin antigens in human cancers and its relationship with malignancy potential[END_REF] and we found MUC2 expressing cells in 5 of 15 (33%) mucinous CCL cases as well. In concordance with our finding, O'Connell et al. [START_REF] O'connell | Altered mucin expression is a field change that accompanies mucinous (colloid) breast carcinoma histogenesis[END_REF] described positive cells in mucinous DCIS, ADH and UDH surrounding mucinous carcinoma. Also "normal epithelium" was described to show expression, but on review of the figures of "normal epithelium" in their 1999 paper we would suggest that this in fact concerns a columnar cell lesion. We have never seen MUC2 luminal or intracytroplasmic staining in normal breast tissue. These data suggest that changes like MUC2 expression already occur in precursor lesions. We also found intraluminal mucin staining in the present study, whereas other studies using other antibodies did not show this staining pattern. This indicates that different MUC2 antibody clones can provide different staining patterns [START_REF] Matsukita | Expression of mucins (MUC1, MUC2, MUC5AC and MUC6) in mucinous carcinoma of the breast: comparison with invasive ductal carcinoma[END_REF][START_REF] Rakha | Expression of mucins (MUC1, MUC2, MUC3, MUC4, MUC5AC and MUC6) and their prognostic significance in human breast cancer[END_REF][START_REF] Adsay | Pathogenesis of colloid (pure mucinous) carcinoma of exocrine organs: Coupling of gel-forming mucin (MUC2) production with altered cell polarity and abnormal cell-stroma interaction may be the key factor in the morphogenesis and indolent behavior of colloid carcinoma in the breast and pancreas[END_REF][START_REF] O'connell | Altered mucin expression is a field change that accompanies mucinous (colloid) breast carcinoma histogenesis[END_REF][START_REF] Chu | Mucin expression in mucinous carcinoma and other invasive carcinomas of the breast[END_REF][START_REF] Yonezawa | Expression of mucin antigens in human cancers and its relationship with malignancy potential[END_REF] .

We showed that the mucinous CCLs are typically biopsied for microcalcifications as abnormality on mammography, which has been described for mucocele-like lesions by others as well 3, 4, 14-16, 29, 32 . Columnar cell lesions in general also present by microcalcifications, which may further suggest a relationship in the origin of the lesions with both unfolded TDLU.

It is yet not fully clear what the implications are of finding a mucinous CCL in a CNB. For mucocele-like lesions found on CNB it has been suggested that excision should follow to exclude the presence of invasive carcinoma [START_REF] Jacobs | Nonmalignant lesions in breast core needle biopsies: to excise or not to excise?[END_REF] , although others suggest that CNBs represent such lesions reliably [START_REF] Renshaw | Atypical ductal hyperplasia in breast core needle biopsies: Correlation of size of the lesion, complete removal of the lesion, and the incidence of carcinoma in follow-up biopsies[END_REF][START_REF] Wang | Can core biopsy reliably diagnose mucinous lesions of the breast?[END_REF] . For columnar cell lesions with atypia, excision biopsies are usually recommended, and for those without atypia only follow-up is advised [START_REF] Feeley | Columnar cell lesions of the breast[END_REF] . Since we did not find any events during follow up, we do not suggest any other approach for mucinous CCL until further studies better address the natural course of mucinous CCLs.

In conclusion, mucinous CCL and mucinous ADH-CCL are rare lesions occurring in about 0.5% of breast CNB, usually presenting during the usual age range with microcalcifications.

These lesions might play a role in the mucinous progression spectrum, which is underlined by the high incidence of mucinous CCLs in mucinous carcinoma cases. Although further studies are needed to more extensively explore the natural course of these lesions, the short term progression risk seems to be low.

  pathology archiving system (PALGA), all breast core needle biopsies (CNB) from women were identified from the archives of the University Medical Centre of Utrecht (UMCU) from January 2001 until May 2008 (except for 2006 because these slides were misplaced), from the Rivierenland Hospital Tiel (RHT) between 2001 until 2005; from the Antonius Hospital Nieuwegein (AHN), and from the Mesos Hospital Utrecht (MHU) from 2002 until 2006, in The Netherlands.

  grade DCIS. Moreover, the presence of microcalcifications in the CCLs was noted in each specimen.

  We obtained follow-up information from all cases until November 2009 by a search in the local and national PALGA Pathology databases and hospital information systems for the last performed mammography and last moment of follow-up. The national PALGA searches reveal pathology reports from all Dutch pathology labs, making sure we would not miss a more advanced lesion in the follow up diagnosed in other centers. The H&E slides of the subsequent surgical excision biopsies from the mucinous CCLs were reviewed by the same observers (AM and PJvD). We noted the presence of a mucinous columnar cell lesion or other lesions in these excisions. Besides, in the subsequent surgical excision biopsies of all other CCL cases in the initial CNBs (apart from those with the mucinous variant), we looked for mucinous CCLs. It was noted for which CCLs invasive carcinoma or DCIS had occurred in the past.

Fig. 1 Fig. 2 Fig. 3 Fig. 4 Fig. 5 Fig. 6 Figure 7 Figure 8

 12345678 Fig. 1 Mucinous columnar cell change of the breast with mucinous dissection (20x)

  . Besides, O'Connell et al. described MUC2 expression of cells in mucinous DCIS, ADH and UDH in the vicinity of mucinous carcinoma as well, implying that the

	F o r
	P
	e
	e r
	R
	e
	v i e
	w

Table 1

 1 Number of core needle biopsies (CNBs) and mammographic abnormalities of columnar cell lesions (CCLs) and atypical hyperplasia originating in CCL (ADH-CCLs)

	Number of cases (%)

Table 2

 2 Nature of mucinous columnar cell lesions (CCL) found in a consecutive series of 4164 breast core needle biopsies and accompanying mammographic abnormalities (CCC=columnar cell change, CCH-columnar cell hyperplasia; CCL-A=CCL with atypia; ADH-CCL, atypical ductal hyperplasia originating in a CCL)

	Number of
	cases (%)

Table 3

 3 Frequency of mucinous and non-mucinous CCLs in 46 mucinous versus 46 invasive

	ductal carcinoma cases of the breast
		Ductal	Mucinous	Chi-square
		carcinoma	carcinoma	test
	Non-mucinous CCL	15	6	p = 0.004
	Mucinous CCL	4	13
	Total	19	19

Table 4

 4 Expression of MUC2 in 15 CCLs of different grade seen in breast core needle

	biopsies
			Intraluminal mucin		Cytoplasmic expression
		-	+	-	+
	Mucinous CCC	2	5	5	2
	Mucinous CCH	0	3	3	0
	Mucinous CCL-A	1	1	2	0
	Mucinous ADH-CCL 0	3	0	3
	Total 3	12 (80%)	10	5 (33%)

Published on behalf of the British Division of the International Academy of Pathology Histopathology


A Cunningham

W Wang

E Zio

D Allanson

A Wall

J Wang

Application of Delay-Time Analysis via Monte Carlo Simulation

Keywords: Delay-time, inspection, maintenance, ship operations, simulation

This paper presents a methodology for the application of delay time analysis via Monte Carlo Simulation. The aim of the paper is to demonstrate the efficacy and worth of delay-time analysis and how the application can provide engineers with more information when making maintenance decisions. The methodology has been developed and applied to two case studies.

Introduction

Maintenance is a huge area of interest and research for engineers. A number of papers based on maintenance strategy and decision have been published [START_REF] Barbera | A Condition Based Maintenance Model with Exponential Failures and Fixed Inspection Intervals[END_REF] [START_REF] Qi | Scheduling Maintenance on a Single Machine[END_REF] (Wang et al. 2000) [START_REF] Majid | Reliability Data Analysis and Modelling of Offshore Oil Platform Plant[END_REF] [START_REF] El-Haram | Practical Application of RCM to Local Authority Housing: A Pilot Study[END_REF] (Emblemsvåg & Tonnig 2003) [START_REF] Beebe | Condition Monitoring of Steam Turbines Performance by Analysis[END_REF]) [START_REF] Backlund | RCM Introduction: Process and Requirements Management Aspects[END_REF] [START_REF] Wang | A Stochastic Maintenance Management Model with Recovery Factor[END_REF]. Maintenance costs form a significant part of the overall operating costs in ship operations [START_REF] Mokashi | A Study of Reliability-Centred Maintenance in Marine Operations[END_REF]. [START_REF] Pillay | Technology and Safety of Marine Systems[END_REF] defined maintenance as the combination of all technical and administrative actions, including supervision actions, intended to retain an entity in, or restore it to a state, in which it can perform a required function. The International Safety Management (ISM) Code states that all ship operators 'should establish and implement procedures to identify equipment and technical systems the sudden operational failure of which would result in hazardous situations' (ISM 2002). In meeting these requirements the company should ensure that:  Inspections are held at appropriate intervals.

 Any non conformity is reported with its possible cause, if known.

 Appropriate corrective actions are taken.

 Records of these activities are maintained. [START_REF] Soncini | PMS on ships: Dream or Reality?[END_REF] suggested that most ship owners understand the need of having good control over accounting and purchasing and are found to be at the same level as their land based counterparts; however the same cannot be said when it comes to maintenance and stock control. [START_REF] Pinelton | Case Study: RCM -Yes, No or Maybe?[END_REF] introduced the 'maintenance concept', defined as the set of various maintenance interventions (corrective, preventive, condition based, etc.) and the general structure in which these interventions are brought together. The total cost of maintenance is difficult to calculate due to the number of factors involved. shows the 'optimum' maintenance concept presented in Pintelon (1999). The efficiency of a concept is dependent upon the input. Ultimately all maintenance concepts are dependent upon appropriate information being available concerning equipment. To enable marine engineers to make educated informed decisions concerning maintenance decisions methods must be developed which provide the marine engineer with information about unreliability, availability and downtime.

It has been shown in a previous paper [START_REF] Cunningham | Application of Monte Carlo Method in the Marine Environment[END_REF]) that Monte Carlo Simulation can be applied in the marine environment to give information about system unreliability based on system failure rates. Input variables and maintenance decision can be 'tested' within the simulation and the effects on system unreliability assessed. In this paper Delay-Time Analysis (DTA) methods will be implemented using Monte Carlo Methods to automate the process and produce results. DTA can be easily achieved through simulation methods but limited work exists that outlines a methodology to demonstrate this method. The simulation method can be used as a way of validating studies carried out using analytical DTA. The paper briefly presents the simplest delay-time model and a methodology for a simulation based approach is developed. Extensions have been made to the analytical delay-time method in order to relax the simplifying assumptions at a cost of increased mathematical complexity. It will be shown in the following paper that DTA via simulation can relax modelling assumptions in the same way, with very little added complexity.

Effectiveness Efficiency

Maintenance concept output input

Background

The Delay-Time Concept

The majority of current reliability and maintenance practice is based on time to first failure, or time between failures. [START_REF] Christer | Developments in Delay-Time Analysis for Modelling Plant Maintenance[END_REF] This example demonstrates that assuming we can model the way in which defects arrive, referred to as the arrival rate of defects kf, and their delay-time h, the delay-time analysis concept can be applied to understand the relationship between inspection frequency and system failures [START_REF] Christer | A Delay-Time Based Maintenance Model of Multi-Component System[END_REF].

Here the simplest delay time model used in the literature is briefly presented. It is asssumed that there is a complex plant, or multi-component plant which has a large number of components with many failure modes, and the correction of one defect or failure has nominal impact in the steady state upon the overall plant failure characteristics. The following assumptions are given for the basic complex plant maintenance modelling scenario:

1. An inspection takes place every T time units and requires 5. The delay time, H , of a random defect is described by a pdf.

) (h f , cdf. ) (h F
, and is independent of the initial point U . 6. Failure will be repaired immediately at an average f d .

7. The plant has operated sufficiently long since new to be considered effectively in a steady state.

8. Defects and failures only arise whilst plant is operating.

These assumptions characterise the simplest non-trivial inspection maintenance problem [START_REF] Christer | A Delay-Time Based Maintenance Model of Multi-Component System[END_REF]. It is now possible to proceed to construct the mathematical model of the relationship between T and an objective function of interest.

From assumptions 1-4, it is obvious that the number of system failures is identical and independent over each inspection interval, and it is possible to simply study the behaviour of such a failure process over one interval, e.g. the first interval various steps contained within the analysis and shows the information required before an analysis can be conducted. In the following, some important steps in the methodology will be explained and expanded upon. Where appropriate, examples will be given to aid the description.

  T f dt t F T N E 0 ) ( )] ( [ 

Development of Methodology

Delay-time simulation involves the consideration of a number of defects and associated delay-times within a given time line. It is assumed that in order for a breakdown to occur, there exists a defect, u, which is a pre-curser to failure. Each u value has an associated delay-time, h, that represents a time window, in which, if normal inspection activity occurs the defect could be recognised and the systems transition into a failed state prevented. Simulation of the delay-time involves consideration of the system over a mission time, T m . T m should be sufficiently long such that downtime due to breakdown and inspection can be considered negligible. The process involves the estimation of a suitable distribution of defects, f(u) and a suitable distribution of delay-times, f(h). The program can be described in the following steps:

1. Generate a value, U 1 , which represents a time of defect, where f(u) is the probability density function of the defect time.

2. Generate an associated delay-time, h 1 , which represents the opportunity window in which inspection could arrest a developing failure, where f(h) is the probability density function of the delay-time.

3. Perform a test to see if the defect is found at the time of inspection.

4. Generate the next defect time, u 2 , from the point u 1 and an associated delay-time h 2 .

Repeat step 3.

The process outlined above is repeated until the cumulative value, CU, is greater than the mission time T m where, DELAY ALGORITM then performs the test which decides if this particular combination of u and h leads to a failure or a breakdown. AVG_LOOP is a second counting loop which repeats the process a set number of times, N, for a given value of T. At the end of each iteration of AVG_LOOP, the array which contains the average values is updated. This averaging process allows more accuracy in the final results. The Monte Carlo Simulation returns the total expected failures for the whole mission at each value of T considered. In the following the Delay Algorithm will be explained in more detail. It takes one of three forms depending upon the analysis, section 3.1.1 presents the Delay Algorithm for perfect inspections, section 3.1.2 presents the Delay Algorithm when imperfect inspections are considered and section 3.1.3 presents the Delay Algorithm when imperfect inspections and imperfect repairs are considered.

𝐶𝑈 = 𝑢 1 + 𝑢 2 +. . . +𝑢 𝑛

Delay Algorithm -Perfect Inspection

The The algorithm performs the calculation, 𝑅𝐸𝐿_𝐼𝑁𝑇 -𝐻 = 14 -0.7 = 13.3

The test is performed to see if 𝐶𝑈 ≥ (𝑅𝐸𝐿_𝐼𝑁𝑇 -𝐻) or 𝐶𝑈 < (𝑅𝐸𝐿_𝐼𝑁𝑇 -𝐻), if the first condition is found to be true an inspection occurs, if the latter is found to be true a breakdown occurs. In the case of breakdown the counter DTE, is increased by one. In both cases of inspection and failure the present u value is added to the cumulative value CU.

Delay Algorithm -Imperfect Inspection

The flowchart form of the algorithm for imperfect inspections is shown in fig 13. The presupposition of perfect inspections has been relaxed. The algorithm now takes into account the probability of an inspection being perfect or imperfect which is preset by the analyst.

All the assumptions previously outlined for analysis still hold true apart from the assumption of perfect inspection. In the case of imperfect inspection it is assumed that at the point of inspection there is a probability, 𝑟, that a defect present will be identified. Conversely there is a probability, 1 -𝑟, that a defect will go unnoticed at inspection and will continue to develop into a full breakdown. [START_REF] Christer | Developments in Delay-Time Analysis for Modelling Plant Maintenance[END_REF] demonstrated how the analytic model can be extended to include imperfect inspections. It should be noted that imperfect inspection when using analytical methods is achieved at the cost of a significant increase in mathematical complexity. The simulation does not suffer from the same increase in complexity. It can be seen from the flowchart shown in fig 13 that the flowchart for imperfect inspection is very similar to the flowchart for perfect inspections.

The inclusion of imperfect inspections into the simulation model is achieved through the introduction of a discrete distribution which represents the probability of perfect and imperfect inspections. The distribution is made of two distinct intervals, 0,1 -𝑟 and 1 -𝑟, 1 . The random number, RI, is called where 𝑅𝐼 = 𝑈~ 0,1 and a test is performed to examine in which interval RI falls. This test decides whether a defect is recognised and repaired at inspection or unnoticed and left to develop into a breakdown failure. For the analysis in the following case studies inspections are considered imperfect 10% of the time. Woods (1984) suggest that in emergency situations this incorrect inspection rate could be as high as 60 %. The value of 10% in light of this can be considered appropriate as the inspections do not take place under emergency conditions.

Delay Algorithm -Imperfect Inspections, Imperfect Repair

Imperfect repair involves the consideration of delay-time analysis with a non-homogeneous defect arrival rate, 𝑘 𝑓 . The assumption that 𝑘 𝑓 is constant is a reasonable assumption for most systems that have been running for a sufficiently long period to be considered mature. Imperfect repair, first considered by Brown & Proschan (1983), can be closely linked to models considering 'minimal repair The model developed in this work considers the effect of minimal repair after an inspection action. It is still assumed that in the case of a breakdown failure the repair of components is perfect and the system is put back online in a 'good as new state'. After a breakdown repair the system is put back online with the original steady-state arrival rate of defects, 𝑘 𝑓 . To examine the effect of non steadystate conditions it is assumed that when a defect is identified at inspection and the defect subsequently repaired, this repair action is non-perfect. This non-perfect repair action has the effect of increasing the arrival rate of defects by 20%. The flowchart form of the algorithm for imperfect inspections with imperfect repair is shown in fig 14.

Calculate Expected Downtime over T

The Monte Carlo Simulation outlined provides the total expected number of failures over a given mission time. The equation for downtime per unit time requires the expected value of failures over T.

In order to achieve this, results given by the simulation have to be divided by N. N is equal to the total number of inspections, T, possible within the given mission time, T m , i.e. 𝑁 = 𝑇 𝑚 𝑇 .

Case Studies

With the intent of demonstrating the method for DTA via simulation two case studies are presented.

In the first the data for the case study was taken from an existing journal paper [START_REF] Pillay | Maintenance Study of Fishing Vessel Equipment Using Delay-Time Analysis[END_REF]). In the second case study a new model is presented based on a centrifugal pump, where repair data is based on OREDA data (OREDA 2002) and expert judgement.

Fishing Vessel Case Study

The delay-time model is based on the operation of a main hydraulic winch operating system on board a fishing vessel. The vessel has length overall of 60m and gross tonnage of 1266. The data for the analysis is taken directly from [START_REF] Pillay | Maintenance Study of Fishing Vessel Equipment Using Delay-Time Analysis[END_REF] and is shown in table 1. The downtime for breakdown repair takes into account any delays caused while waiting for spares to be sent to the vessel.

Initial Modelling Assumptions

When performing the analysis for the case study the following modelling assumptions were made.

 Inspections take place at regular intervals of T hours and each inspection is identical.

 The arrival rate of defects is constant and distributed according to an exponential probability density function.

 Failures are repaired instantaneously and the system is returned to a 'good as new' state.

 The mission time is set to 10 years and is sufficiently large that downtime due to breakdown and inspection during the analysis can be considered negligible.

 Inspections are perfect in that any defect present will be identified and the failure arrested within the inspection period.

Case study Results

Fishing Vessel -Perfect Inspections

Fishing Vessel -Imperfect Inspections

The analysis was conducted using a FORTRAN programme in the way outlined in the methodology previously. Fig 19 shows the results of the analysis. The programme was run using shape and scale parameters, k=10 and λ=5, which produced the most definitive result for perfect inspection. It can be seen from the graph that when imperfect inspections are considered the value of minimum DT per unit time is increased. The recommendation for the optimum inspection interval remains appropriate at T=9 hours giving a downtime per unit time of 0.041 hours.

Cooling System -Centrifugal Pump -Imperfect Inspections

The programme was run using shape and scale parameters, k=10 and λ=5, which produced the most definitive result for perfect inspection.

Fishing Vessel -Imperfect Repair

The programme was run using shape and scale parameters, k=10 and λ=5, which produced the most definitive result for perfect inspection.

Cooling System -Centrifugal Pump -Imperfect Repair

The programme was run using shape and scale parameters, k=10 and λ=5, which produced the most definitive result for perfect inspection.

Sensitivity Analysis

A sensitivity analysis provides a way of partially validating a model. For this model three axioms are detailed and must be satisfied before the sensitivity analysis can be considered complete.

 An increase in the arrival rate of defects should result in a proportional increase in the DT per unit time.

 Further increase in the arrival rate of defects should reflect a consistent increase in the DT per unit time.

 An increase in more than one input parameter should result in a larger increase in DT per unit time than that caused by an increase in a single input parameter.

The sensitivity analysis was conducted on a single case presented previously, with perfect inspections and perfect repairs. The more complex cases involving imperfect inspection and repair are extensions of this model; therefore, partial validation of this model will also provide partial validation of the more complex cases. The case study involving the input parameters for the cooling system centrifugal pump was used. The models for both the fishing vessel and the centrifugal pump both

Discussion

The analysis programme can be easily altered to consider a different set of equipment, with different input parameters. The only limitation to the simulation method is the ability of the programmer to generate random numbers distributed to different distributions. The method gains accuracy when the mission time is set at larger values. This is often at the expense of time to compute simulation results.

As computers increase in both speed and processing power this will become less of a problem, however the analyst should always give careful consideration to the suitability of the mission time.

Short mission times will produce more results in a shorter period of time but this may be at the expense of accuracy. Conversely exceptionally long mission times will produce very accurate results but may prove unrealistic in terms of an average component lifetime and may also prove impractical in terms of processing time.

When considering the results of any analysis reflection on the propriety of the modelling assumptions must be made. The assumption that all inspections that take place are perfect and that all defects are recognised and corrected is improbable. However the simulation programme can easily be amended to take into account the case of imperfect inspection. In order to examine the impact of imperfect inspections the analysis was repeated with the premise that inspections were only perfect 90% of the time. In the remaining 10% the defects went unnoticed at inspection and developed into full breakdown failures. It can be seen from the results shown that imperfect inspection intervals result in an increase of the DT per unit time. The DT per unit time increases as the amount of imperfect inspections increase, reducing the amount of imperfect inspections reduces the DT per unit time. The optimal inspection interval remains unchanged. Further increase or reduction in the amount of imperfect inspections has a similar affect of 'shifting' the curve vertically away from or towards the perfect case.

The assumption that the system is returned to 'good as new' after inspection and repair is also one that seems unrealistic. This may not prove to be the case in real life, systems may be put back into service in a degraded state after inspection or repair. This is ultimately dependant upon the experience and skill of the maintenance personnel and the quality of the replacement parts. To examine the effect of imperfect repair, the analysis was repeated with the assumption that after an inspection and subsequent corrective action the system is put back online with an arrival rate of defects increased by 20%. It can be seen from the results that this increases the level of downtime per unit time achieved.

The optimal inspection periods remain unchanged.

All of the simulations implemented in the paper use a random seed function whenever a random number is generated. This provides suibtably random results. Due to the nature of the pseudo-random number generator function used in the simulation, non random seed values would produce repeated results and the generation of the random would become some what deterministic in nature.

The results of both models which concerned imperfect inspection and imperfect repair are logical. If inspections are imperfect then there is an increased chance for system breakdown, this is reflected in the increase of downtime per unit time. In the case of imperfect repair the arrival rate of defects increases, this leads to more defects and results in an increase of the downtime per unit time. The strength of the MCM of DTA is the method's ability to deal with different situations in a logical and straightforward way. The inclusion of imperfect inspection and repair comes at the cost of a few additional lines of code. To consider the same problems using traditional analytical methods would result in a significant increase in mathematical rigour. Marine engineers, having often achieved their qualifications in a vocational system, often lack the mathematical skills necessary to perform such an analysis via analytical methods. The simulation method presented circumnavigates this knowledge gap and provides a useful tool for marine engineers in an accessible way.

The results were as expected in that a minimum downtime was found. Before running the analysis the analyst had only very general ideas of expected reults. The model is only partially validated by the sensitivity study and further work is needed to produce a theoretical result.

The method could prove to be a very useful tool in defining inspection regimes for particular pieces of equipment. For the method to be fully effective an inspection regime would have to be implemented to provide the simulation program with accurate historical failure data. The more data gathered the more accurate and effective the analysis would become. Any decisions made concerning the maintenance regime onboard will ultimately be decided by the owner/operator of the vessel. It may be the case in certain situations that it is not possible to carry out inspection at the recommended interval. It is not always convenient or even safe to take certain systems offline during passage, this would obviously be system specific and engineering judgement would play a large part in how the inspection regime could be altered or adjusted. The decision to implement DTA will depend upon existing operating and maintenance culture onboard.

Conclusion

The purpose of this work was to demonstrate an alternative method for DTA other than traditional analytical methods. Previous research work reported on DTA is often arduous in terms of the mathematical models presented. It has demonstrated the benefit of the method but the esoteric nature of the mathematical models, has often prevented engineers in industry from implementing the method.

The intention of the researchers of this work was to present a methodology of achieving the same results in a more accessible way to a wider range of engineers. Based on the evidence of the results presented the methodology outlined for performing the analysis will provide optimal inspection periods for a given set of data. This work also demonstrates the power and flexibility contained within the MCM to consider a number of different models and methods. A need is also identified for ship owners/operators to invest more time into the collation of failure data specific to their vessels.

Different vessels operating in different areas and conditions will display different failure characteristics. The collection of failure data and its use in the analysis of systems with respect to reliability and appropriate maintenance scheduling could only prove beneficial to ship operators.

Fig 1 .

 1 Fig 1. Optimum maintenance concept.

  published a review considering the developments in DTA, stating that 'maintenance concepts based on Reliability Centered Maintenance (RCM) or Total Productive Maintenance (TPM) are prescriptive and often lack scientific concept, testing, verification or validation'. Delay-time modelling is a concept which has been developed to be relevant in the operating culture of today's industry (Christer 1999). DTA provides engineers with a tool which can help to minimise downtime, D(T) of a machine or plant item, based on an inspection period, T. The delay-time concept bifurcates the failure process as shown in fig 2.

Fig 2 .

 2 Fig 2. Diagram showing the delay-time concept.

Fig 3 .

 3 Fig 3. Diagram showing the failure points.

Fig 4 .

 4 Fig 4. Diagram showing the failure points and the initial points.

  are perfect in that all (and only) defects present are identified.

3.

  Defects identified are repaired during the inspection period.

4.

  Defects arise according to a Homogeneous Poisson Process (HPP) with the rate of occurrence of defects,  , per unit time.

Fig 5

 5 Fig 5 Proposed methodology for performing a delay-time analysis via Monte Carlo Simulation

Fig 6 Fig 7

 67 Fig 6 shows the generation of a number of u values within the mission time.

Fig 9

 9 Fig 9 Flowchart representing the Delay Algorithm for perfect inspections

  at failure'(Barlow & Proschan 1965),(Blumenthal et al. 1976). Further study and extension of the Brown & Proschan model was conducted byWhitaker & Samaniego (1989). Baker &Wang (1993) considered delay-time analysis where the assumption of constant 𝑘 𝑓 is relaxed. The model considers the effect of component age on the arrival rate of defects and the consequence of inspection activity and its possible hazardous or beneficial effect on the lifetime of a component.

Fig

  Fig 13 Flowchart representing the Delay Algorithm for imperfect inspections

  Fig 15 shows a schematic of the main hydraulic piping system.

Fig

  Fig 15 Hydraulic winch operating system

Fig 17 Fig 18

 1718 Fig 17 Graph showing DT per unit time against T

Fig 19

 19 Fig 19 Graph showing DT per unit time against T considering imperfect inspections

Fig 20

 20 Fig 20 Graph showing DT per unit time against T considering imperfect inspections

Fig 21

 21 Fig 21 Graph showing DT per unit time against T considering imperfect inspections with perfect repair, imperfect inspections with imperfect repair and perfect inspections with perfect repair

Fig 22

 22 Fig 22 Graph showing DT per unit time against T considering imperfect inspections with perfect repair, imperfect inspections with imperfect repair and perfect inspections with perfect repair

  Fig 23 Graph showing the results of the sensitivity analysis

Table 1

 1 Table showing the input parameters for the analysis

	Inspection Downtime ds 0.25 hrs
	Breakdown Downtime db 108 hrs
	Arrival rate of defects kf 0.0223 hrs -1

Table 3

 3 Table showing the shape and scale parameters of the Weibull distributions used in the analysis

	k	Λ
	10	5
	8	6
	3	10
	2	20

The analysis is performed on one centrifugal pump which services the main cooling system. To carry out the analysis a number of input variables were required. The downtime due to breakdown was taken from OREDA 2002 and set equal to 168 hours or 7 days which allows for any logistical delay in spare part procurement. For the downtime due to inspection the expert opinion of an experienced marine chief engineer was used. A detailed description of the chief engineer's industrial experience and academic qualifications is listed in the appendix. Daily inspection of the centrifugal pump involves visual inspection of suction and discharge pressure, audio inspection for any abnormal noise and electrical inspection of the current being drawn by the electric motor. The chief engineer suggested that this daily inspection on average would take 10-15 minutes. In light of this the downtime due to inspection was taken as 12.5 minutes or 0.2083 hours. When considering the arrival rate of defects it is argued that the failure rate of a system and the arrival rate of defects are intrinsically linked. In order for this to be true the component or system would have to be operated under a breakdown maintenance policy. OREDA data is not presented for systems operating under a breakdown maintenance regime. However for the purpose of the analysis it is assumed that the OREDA failure data for a centrifugal pump and the arrival rate of defects are equivalent. OREDA gives the failure rate per 10 6 hours for a centrifugal pump, in all modes of failure, as 1277.00. This is based on a population of 350 pumps over 59 installations.

Estimation of Delay-time Probability Density Function

In a case study based on a specific system the probability density function of the delay-time would be estimated using historical failure data and operator questionnaires. This process in itself takes a great deal of time and logistical work. The purpose of this work was to demonstrate the simulation method of DTA, therefore the analysis was performed using a number of different Weibull distributions for the delay-time and may not represent accurately the true distributions of the delay-times for the real life systems. Table 3 shows the shape and scale parameters used for the different analyses. A number of shape and scale parameters are used to give an idea of their effect on the analysis.

There is certainly huge scope for further work especially when the simplifying assumptions are considered. In the models presented two of these assumptions were relaxed. The more interesting of the two is the assumption that defect arrival rate is constant. In the analysis the arrival rate was changed as a result of different inspection and breakdown actions, however the arrival rate always obeyed the same distribution. Further work could be done to examine the effect of changing the distribution of the arrival rate of defects throughout the analysis. There is also scope for work considering the age of components and the effect of component age on defect arrival rate.