Frédéric Blanqui

Type safety of rewrite rules in dependent types

Keywords: 2012 ACM Subject Classification Theory of computation → Type theory, Theory of computation → Equational logic and rewriting, Theory Keywords and phrases subject-reduction, rewriting, dependent types Digital Object Identifier 10.4230/LIPIcs.FSCD.2020.11

The expressiveness of dependent type theory can be extended by identifying types modulo some additional computation rules. But, for preserving the decidability of type-checking or the logical consistency of the system, one must make sure that those user-defined rewriting rules preserve typing. In this paper, we give a new method to check that property using Knuth-Bendix completion.

Introduction

The λΠ-calculus, or LF [START_REF] Harper | A framework for defining logics[END_REF], is an extension of the simply-typed λ-calculus with dependent types, that is, types that can depend on values like, for instance, the type V n of vectors of dimension n. And two dependent types like V n and V p are identified as soon as n and p are two expressions having the same value (modulo the evaluation rule of λ-calculus, β-reduction).

In the λΠ-calculus modulo rewriting, function and type symbols can be defined not only by using β-reduction but also by using rewriting rules [START_REF] Wrightson | Term rewriting systems[END_REF]. Hence, types are identified modulo β-reduction and some user-defined rewriting rules. This calculus has been implemented in a tool called Dedukti [9].

Adding rewriting rules adds a lot of expressivity for encoding logical or type systems. For instance, although the λΠ-calculus has no native polymorphism, one can easily encode higher-order logic or the calculus of constructions by using just a few symbols and rules [START_REF] Cousineau | Embedding pure type systems in the lambda-Pi-calculus modulo[END_REF]. As a consequence, various tools have been developed for translating actual terms and proofs from various systems (Coq, OpenTheory, Matita, Focalize, . . . ) to Dedukti, and back, opening the way to some interoperability between those systems [START_REF] Assaf | Dedukti: a logical framework based on the λΠ-calculus modulo theory[END_REF]. The Agda system recently started to experiment with rewriting too [START_REF] Cockx | Sprinkles of Extensionality for Your Vanilla Type Theory[END_REF].

To preserve the decidability of type-checking and the logical consistency, it is however essential that the rules added by the user preserve typing, that is, if an expression e has some type T and a rewriting rule transforms e into a new expression e , then e should have type T too. This property is also very important in programming languages, to avoid some errors (a program declared to return a string should not return an integer).

When working in the simply-typed λ-calculus, it is not too difficult to ensure this property: it suffices to check that, for every rewriting rule l → r, the right-hand side (RHS) r has the same type as the left-hand side (LHS) l, which is decidable.

The situation is however much more complicated when working with dependent types modulo user-defined rewriting rules. As type-checking requires one to decide the equivalence 11:2 Type safety of rewrite rules in dependent types of two expressions, it is undecidable in general to say whether a rewriting rule preserves typing, even β-reduction alone [START_REF] Saillard | Type checking in the Lambda-Pi-calculus modulo: theory and practice[END_REF].

Note also that, in the λΠ-calculus modulo rewriting, the set of well-typed terms is not fixed but depends on the rewriting rules themselves (it grows when one adds rewriting rules).

Finally, the technique used in the simply-typed case (checking that both the LHS and the RHS have the same type) is not satisfactory in the case of dependent types, as it often forces rule left-hand sides to be non-linear [START_REF] Blanqui | Definitions by rewriting in the calculus of constructions[END_REF], making other important properties (namely confluence) more difficult to establish and the implementation of rewriting less efficient (if it does not use sharing).

Example 1.

As already mentioned, Dedukti is often used to encode logical systems and proofs coming from interactive or automated theorem provers. For instance, one wants to be able to encode the simply-typed λ-calculus in Dedukti. Using the new Dedukti syntax1 , this can be done as follows (rule variables must be prefixed by $ to distinguish them from function symbols with the same name):

) → (τ a → τ b ) rule app $a $b ( lam $a ' $b ' $f ) $x → $f $x // β -reduction
Proving that the above rule preserves typing is not trivial as it is equivalent to proving that β-reduction has the subject-reduction property in the simply-typed λ-calculus. And, indeed, the previous version of Dedukti was unable to prove it.

The LHS is typable if f is of type τ (arr a b ), τ (arr a b ) τ (arr a b), and x is of type τ a. Then, in this case, the LHS is of type τ b.

Here, one could be tempted to replace a by a, and b by b, so that these conditions are satisfied but this would make the rewriting rule non left-linear and the proof of its confluence problematic [START_REF] Klop | Combinatory reduction systems[END_REF].

Fortunately, this is not necessary. Indeed, we can prove that the RHS is typable and has the same type as the LHS by using the fact that τ (arr a b ) τ (arr a b) when the LHS is typable. Indeed, in this case, and thanks to the rule defining τ , f is of type τ a → τ b. Therefore, the RHS has type τ b as well.

In this paper, we present a new method for doing this kind of reasoning automatically. By using [START_REF]Automation of reasoning. 2: classical papers on computational logic[END_REF], the equations holding when a LHS is typable are turned into a convergent (i.e. confluent and terminating) set of rewriting rules, so that the type-checking algorithm of Dedukti itself can be used to check the type of a RHS modulo these equations.

Outline. The paper is organized as follows. In Section 2, we recall the definition of the λΠ-calculus modulo rewriting. In Section 3, we recall what it means for a rewriting rule to F. Blanqui 11:3 preserve typing. In Section 4, we describe a new algorithm for checking that a rewriting rule preserves typing and provide general conditions for ensuring its termination. Finally, in Section 5, we compare this new approach with previous ones and conclude.

λΠ-calculus modulo rewriting

Following Barendregt's book on typed λ-calculus [START_REF] Barendregt | Lambda calculi with types[END_REF], the λΠ-calculus is a Pure Type System (PTS) on the set of sorts S = { , }: 2 Definition 2 (λΠ-term algebra). A λΠ-term algebra is defined by: a set F of function symbols, an infinite set V of variables, such that V, F and S are pairwise disjoint.

The set T (F, V) of λΠ-terms is then inductively defined as follows: A substitution is a finite map from V to T (F, V). It is written as a finite set of pairs. For instance, {(x, a)} is the substitution mapping x to a.

t, u := s ∈ S | x ∈ V | f ∈ F | λx : t,
Given a substitution σ and a term t, we denote by tσ the capture-avoiding replacement of every free occurrence of x in t by its image in σ.

Definition 3 (λΠ-calculus).

A λΠ-calculus on T (F, V) is given by: a function Θ : F → T (F, V) mapping every function symbol f to a term Θ f called its type (we will often write f : A instead of Θ f = A), a function Σ : F → S mapping every function symbol f to a sort Σ f , a set R of rewriting rules (l, r) ∈ T 2 , written l → r, such that FV(r) ⊆ FV(l).

We then denote by the smallest equivalence relation containing → = → R ∪ → β where → R is the smallest relation stable by context and substitution containing R, and → β is the usual β-reduction relation.

Example 4.

For representing natural numbers, we can use the function symbols N : of sort , and the function symbols 0 : N and s : N → N of sort . Addition can be represented by + : N → N → N of sort together with the following set of rules:

0 + y → y x + 0 → x x + (sy) → s(x + y) (sx) + y → s(x + y) (x + y) + z → x + (y + z)
Note that Dedukti allows overlapping LHS and matching on defined symbols like in this example. (It also allows higher-order pattern matching like in Combinatory Reduction Systems (CRS) [START_REF] Klop | Combinatory reduction systems: introduction and survey[END_REF] but we do not consider this feature in the current paper.)

2 PTS sorts should not be confused with the notion of sort used in first-order logic. The meaning of these sorts will be explained after the definition of typing (Definition 5). Roughly speaking, is the type of objects and proofs, and is the type of set families and predicates. (fun)

Θ f : Σ f f : Θ f (var) Γ A : s Γ, x : A x : A (x / ∈ Γ) (weak) Γ t : T Γ A : s Γ, x : A t : T (x / ∈ Γ) (prod) Γ A : Γ, x : A B : s Γ Πx : A, B : s (app) Γ t : Πx : A, B Γ a : A Γ ta : B{(x, a)} (abs) Γ, x : A b : B Γ Πx : A, B : s Γ λx : A, b : Πx : A, B (conv) Γ t : T Γ U : s Γ t : U (T U )
Throughout the paper, we assume a given λΠ-calculus Λ = (F, V, Θ, Σ, R).

Definition 5 (Well-typed terms).

A typing environment is a possibly empty ordered sequence of pairs (x 1 , A 1 ), . . . , (x n , A n ), written x 1 : A 1 , . . . , x n : A n , where the x i 's are distinct variables and the A i 's are terms.

A term t has type A in a typing environment Γ if the judgment Γ t : A is derivable from the rules of Figure 1. An environment Γ is valid if some term is typable in it.

A substitution σ is a well-typed substitution from an environment Γ to an environment Γ , written Γ σ : Γ, if, for all x : A ∈ Γ, we have Γ xσ : Aσ.

Note that well-typed substitutions preserve typing: if Γ t : T and Γ σ : Γ, then Γ tσ : T σ [START_REF] Blanqui | Théorie des types et récriture[END_REF].

A type-checking algorithm for the λΠ-calculus modulo (user-defined) rewriting rules is implemented in the Dedukti tool [9].

We first recall a number of basic properties that hold whatever R is and can be easily proved by induction on [START_REF] Blanqui | Théorie des types et récriture[END_REF]:

Lemma 6. (a) If t is typable, then every subterm of t is typable. (b) is not typable. (c) If Γ t : T then either T = or Γ T : s for some sort s. (d) If Γ t : then t is a kind, that is, of the form Πx 1 : T 1 , . . . , Πx n : T n , . (e) If Γ t : T , Γ ⊆ Γ and Γ is valid, then Γ t : T .
Throughout the paper, we assume that, for all f , Θ f : Σ f . Indeed, if Θ f : Σ f does not hold, then no well-typed term can contain f . (This assumption is implicit in the presentations of LF using signatures [START_REF] Harper | A framework for defining logics[END_REF].)

More importantly, we will assume that → is confluent on the set T (F, V) of untyped terms, that is, for all terms t, u, v ∈ T (F, V), if t → * u and t → * v, then there exists a term w ∈ T (F, V) such that u → * w and v → * w, where → * is the reflexive and transitive closure of →.

F. Blanqui 11:5

This condition is required for ensuring that conversion behaves well with respect to products (if Πx : A, B Πx : A , B then A A and B B ), which in particular implies subject-reduction for → β .

This last assumption may look strong, all the more so since confluence is undecidable. However this property is satisfied by many systems in practice. For instance, → is confluent if the left-hand sides of R are algebraic (Definition 8), linear and do not overlap with each other [START_REF] Van Oostrom | Confluence for abstract and higher-order rewriting[END_REF]. This is in particular the case of the rewriting systems corresponding to the function definitions allowed in functional programming languages such as Haskell, Agda, OCaml or Coq. But confluence can be relaxed in some cases: when there are no type-level rewriting rules [START_REF] Barbanera | Modularity of strong normalization in the algebraic-λ-cube[END_REF] or when the right-hand sides of type-level rewriting rules are not products [START_REF] Blanqui | Definitions by rewriting in the calculus of constructions[END_REF].

When → is confluent, the typing relation satisfies additional properties. For instance, the set of typable terms can be divided into three disjoint classes: the terms of type , called kinds, of the form Πx 1 : A 1 , . . . , Πx n : A n , ; the terms whose type is a kind, called predicates; the terms whose type is a predicate, called objects.

Subject-reduction

A relation preserves typing (subject-reduction property) if, for all environments Γ and all terms t, u and A, if Γ t : A and t u, then Γ u : A.

One can easily check that → β preserves typing when → is confluent [START_REF] Blanqui | Théorie des types et récriture[END_REF]. Our aim is therefore to check that → R preserves typing too. To this end, it is enough to check that every rule l → r ∈ R preserves typing, that is, for all environments Γ, substitutions σ and terms A, if Γ lσ : A, then Γ rσ : A.

A first idea is to require that: The condition (*) is sufficient if the rule left-hand side is a non-variable simply-typed first-order term [START_REF] Barbanera | Modularity of strong normalization in the algebraic-λ-cube[END_REF], a notion that we slightly generalize as follows: Definition 8 (Pattern). We assume that the set of variables is split in two disjoint sets, the algebraic variables and the non-algebraic ones, and that there is an injection from algebraic variables to non-algebraic variables.

A term is algebraic if it is an algebraic variable or of the form f t 1 . .

. t n with each t i algebraic and f a function symbol whose type is of the form Πx

1 : A 1 , . . . , x n : A n , B.
A term is an object-level algebraic term if it is algebraic and all its function symbols are of sort .

A pattern is an algebraic term of the form f t 1 . . . t n where each t i is an object-level algebraic term.

The distinction between algebraic and non-algebraic variables is purely technical: for generating equations (Definition 10), we need to associate a type x to every variable x, and we need those variables x to be distinct from one another and distinct from the variables F S C D 2 0 2 0 11:6 Type safety of rewrite rules in dependent types used in rules. To do so, we split the set of variables into two disjoint sets. The ones used in rules are called algebraic, and the others are called non-algebraic. Finally, we ask the function to be an injection from the set of algebraic variables to the set of non-algebraic variables.

In the rest of the paper, we also assume that rule left-hand sides are patterns. Hence, every rule is of the form f l 1 . . . l n → r, and we say that a symbol f ∈ F is defined if there is in R a rule of the form f l 1 . . . l n → r.

However, the condition (*) is not satisfactory in the context of dependent types. Indeed, when function symbols have dependent types, it often happens that a term is typable only if it is non-linear. And, with non-left-linear rewriting rules, → is generally not confluent on untyped terms [START_REF] Klop | Combinatory reduction systems[END_REF], while there exist many confluence criteria for left-linear rewriting systems [START_REF] Van Oostrom | Confluence for abstract and higher-order rewriting[END_REF].

Throughout the paper, we will use the following simple but paradigmatic example to illustrate how our new algorithm works: Example 9. Consider the following rule to define the tail function on vectors:

tail n (cons x p v) → v
where tail : Πn : N, V (sn) → V n, V : N → , nil : V 0, cons : R → Πn : N, V n → V (sn) and R : .

For the left-hand side to be typable, we need to take p = n, because tail n expects an argument of type V (sn), but cons x p v is of type V (sp).

Yet, the rule with p = n preserves typing. Indeed, assume that there is an environment Γ, a substitution σ and a term A such that Γ tail nσ (cons xσ pσ vσ) : A. By inversion of typing rules, we get V (nσ)

A, Γ A : s for some sort s, V (spσ) V (snσ) and Γ vσ : V pσ. Assume now that V and s are undefined, that is, there is no rule of R of the form V t → u or st → u. Then, by confluence, pσ nσ. Therefore, V pσ A and Γ vσ : A.

Hence, that a rewriting rule l → r preserves typing does not mean that its left-hand side l must be typable [START_REF] Blanqui | Definitions by rewriting in the calculus of constructions[END_REF]. Actually, if no instance of l is typable, then l → r trivially preserves typing (since it can never be applied)! The point is therefore to check that any typable instance of l → r preserves typing.

A new subject-reduction criterion

The new criterion that we propose for checking that l → r preserves typing proceeds in two steps. First, we generate conversion constraints that are satisfied by every typable instance of l (Figure 2). Then, we try to check that r has the same type as l in the type system where the conversion relation is extended with the equational theory generated by the conversion constraints inferred in the first step. For type-checking in this extended type theory to be decidable and implementable using Dedukti itself, we use to replace the set of conversion constraints by an equivalent but convergent (i.e. terminating and confluent) set of rewriting rules.

Inference of typability constraints

We first define an algorithm for inferring typability constraints and then prove its correctness and completeness.

F. Blanqui 11:7

Definition 10 (Typability constraints). For every algebraic term t, we assume given a valid environment ∆ t = y 1 : , y 1 : y 1 , . . . , y k : , y k : y k where y 1 , . . . , y k are the free variables of t.

Let ↑ be the partial function defined in Figure 2.

It takes as input a term t and returns a pair (A, E), written A[E],

where A is a term and E is a set of equations, an equation being a pair of terms (l, r) usually written l = r.

A substitution σ satisfies a set E of equations, written σ |= E, if for all equations a = b ∈ E, aσ bσ.

y ↑ y[∅] f : Πx 1 : T 1 , . . . , Πx n : T n , U t 1 ↑ A 1 [E 1 ] t n ↑ A n [E n ] f t 1 . . . t n ↑ U σ[E 1 ∪ . . . ∪ E n ∪ {A 1 = T 1 σ, . . . , A n = T n σ}]
where σ = {(x 1 , t 1 ), . . . , (x n , t n )} Example 11. In our running example tail n (cons

x p v) → v, we have cons x p v ↑ V (sp)[E 1 ] with E 1 = { x = T , p = N , v = V p}, and tail n (cons x p v) ↑ V n[E 2 ] with E 2 = E 1 ∪ { n = N , V (sp) = V (sn)}. Lemma 12. If Γ Πx 1 : T 1 , . . . , Πx n : T n , U : s then, for all i, Γ i-1 T i : and Γ n U : s, where Γ i = Γ, x 1 : T 1 , . . . , x i : T i .
Proof. Since → is confluent and left-hand sides are patterns, s s iff s = s . The result follows then by inversion of typing rules and weakening.

In particular, because Θ f : Σ f for all f , we have: Corollary 13. For all function symbols f : Πx 1 : T 1 , . . . , Πx n : T n , U and integer i, we have

Γ i-1 f T i : and Γ n f U : Σ f where Γ i f = x 1 : T 1 , . . . , x i : T i .
Lemma 14. For all environments Γ, terms t, x 1 , T 1 , . . . , x n , T n , U, T and substitutions σ for x 1 , . . . , x n , if Γ t : Πx 1 : T 1 , . . . , Πx n : T n , U and Γ tx 1 σ . . . x n σ : T , then U σ T and Γ σ : ∆ n where ∆ n = x 1 : T 1 , . . . , x n : T n .

Proof. Let σ i = {(x 1 , t 1 ), . . . , (x i-1 , t i-1 )}. We proceed by induction on n.

Case n = 0. By equivalence of types.

Case n > 0. By inversion of typing rules and weakening, Γ, ∆ n-1

T n : , Γ tx 1 σ n-1 . . . x n-1 σ n-1 : Πx n : A, B, Γ x n σ : A and B{(x n , x n σ)} T . By induction hypothesis, Γ σ n-1 : ∆ n-1 and (x n : T n σ n-1 )U σ n-1 Πx n : A, B. By substitu- tion, Γ
T n σ n-1 : . By confluence, T n σ n-1 A and U σ n-1 B. Therefore, by conversion, Γ x n σ : T n σ and Γ σ : ∆ n . Now, x n can always be chosen so that

U σ = U σ n-1 {(x n , x n σ)}. Therefore, U σ B{(x n , x n σ)} T .

Lemma 15.

(Correctness) For all algebraic terms t, terms T and sets of equations E, if t ↑ T [E] then, for all valid environments Γ, substitutions θ such that Γ θ : ∆ t and θ |= E, we have Γ t θ : T θ. Type safety of rewrite rules in dependent types (Completeness) For all environments Γ, patterns t, substitutions θ and terms A, if Γ tθ : A, then there are a term T , a set of equations E and a substitution θ extending θ such that t ↑ T [E], θ |= E, Γ θ : ∆ t and A T θ.

Proof.

(Correctness) By induction on t. Case t = y. Then, T = y and E = ∅. By assumption, we have Γ yθ : yθ. Therefore, Γ tθ : T θ.

Case t = f t 1 . . . t n with f : Πx 1 : T 1 , . . . , x n : T n , U , t 1 ↑ A 1 [E 1 ], . . . , t n ↑ A n [E n ]. Then, T = U σ and E = E 1 ∪ . . . ∪ E n ∪ {A 1 = T 1 σ, . . . , A n = T n σ} where σ = {(x 1 , t 1 ), . . . , (x n , t n )}. By Lemma 12, we have Γ i-1 f T i : .
By induction hypothesis, for all i, we have Γ t i θ : A i θ and A i θ T i σ θ. We now prove that, for all i, Γ T i σ θ : and Γ x i σ θ : T i σ θ, hence that Γ σ : Γ i f , by induction on i.

* Case i = 1. Since T 1 : , T 1 is closed and T 1 σ θ = T 1 . Therefore, by weakening, Γ T 1 σ θ : and, by conversion, Γ x 1 σ θ :

T 1 σ θ. * Case i > 1. By induction hypothesis, Γ σ θ : Γ i-1 f . Since Γ i-1 f T i : , by
substitution, we get Γ T i σ θ : . Therefore, by conversion, Γ x i σ θ : T i σ θ. Hence, Γ σ θ : Γ n f . Now, since Γ n f f x 1 . . . x n : U , by substitution, we get Γ t : U σ θ. (Completeness) We first prove completeness for object-level algebraic terms t such that Γ A : , by induction on t.

Case t = y. We take T = y, E = ∅ and θ = θ ∪ {( y, A)}. We have t ↑ T [E], θ |= E and A T θ. Now, Γ y θ : y θ and Γ y θ : . Therefore, Γ θ : ∆ t . Case t = f t 1 . . . t n with f : Πx 1 : T 1 , . . . , x n : T n , U . By Lemma 12, for all i, we have Γ f x i : T i and Γ f T i : , where Γ f = x 1 : T 1 , . . . , x n : T n . By Lemma 14 because Θ f : Σ f for all f , we have A U σθ and Γ σθ : Γ f . Hence, by substitution, for all i, we have Γ t i θ : T i σθ and Γ T i σθ : . Therefore, by induction hypothesis, there are A i , E i and θ i extending θ such that

t i ↑ A i [E i ], θ i |= E i , Γ θ i : ∆ ti and T i σθ A i θ i . Then, let T = U σ, E = E 1 ∪ . . . ∪ E n ∪ {(A 1 , T 1 σ), . . . , (A n , T n σ)}, y θ = yθ if y ∈ FV(t)
, and y θ = y θ i where i is the smallest integer such that y ∈ FV(t i ). Then, we have t ↑ T [E] and A U σθ = T θ. If y ∈ FV(t i ) ∩ FV(t j ), then y θ i = yθ = y θ j since θ i and θ j are both extensions of θ. Now, if Γ y θ i : y θ i and Γ y θ j : y θ j then, by equivalence of types, y θ i y θ j . Therefore, θ |= E and Γ θ : ∆ t . Let now t be a pattern. By definition, t is of the form f t 1 . . . t n with f : Πx 1 : T 1 , . . . , x n : T n , U and each t i an object-level algebraic term. As we have seen above, for all i, we have Γ t i θ : T i σθ and Γ T i σθ : . Therefore, by completeness for object level algebraic terms, there are A i , E i and θ i extending θ such that

t i ↑ A i [E i ], θ i |= E i , Γ θ i : ∆ ti and T i σθ A i θ i .
We can now conclude like in the previous case.

Example 16. In our running example tail n (cons x p v) → v, we have seen that

cons x p v ↑ V (sp)[E 1 ] with E 1 = { x = T , p = N , v = V p}, and tail n (cons x p v) ↑ V n[E 2 ] with E 2 = E 1 ∪ { n = N , V (sp) = V (sn)}.
This means that, if σ is a substitution and (tail n (cons x p v))σ is typable, then σ |= E 2 . In particular, V (spσ) V (snσ).

Type-checking modulo typability constraints

For checking that the right-hand side of a rewriting rule l → r has the same type as the left-hand side modulo the typability constraints E of the left hand-side, we introduce a new

F = F ∪ {x | x ∈ FV(l)} ∪ { x | x ∈ FV(l)} V = V -({x | x ∈ FV(l)} ∪ { x | x ∈ FV(l)}) Θ = Θ ∪ {(x, x) | x ∈ FV(l)} ∪ {( x, ) | x ∈ FV(l)} Σ = Σ ∪ {(x, ) | x ∈ FV(l)} ∪ {( x, ) | x ∈ FV(l)} R = R ∪ E ∪ E -1
, where l = r ∈ E -1 iff r = l ∈ E. We denote by l,E the conversion relation of Λ l,E , and by l,E its typing relation.

Λ l,E is similar to Λ except that the symbols of {x | x ∈ FV(l)} ∪ { x | x ∈ FV(l)} are not variables but function symbols, and that the set of rewriting rules is extended by E ∪ E -1 which, in Λ l,E , is a set of closed rewriting rules (rules and equations are synonyms: they both are pairs of terms).

Lemma 18. For all patterns l, sets of equations E and substitutions σ in Λ, and for all terms t, u in Λ l,E , if σ |= E and t l,E u, then tσ uσ. 4Proof. Immediate as each application of an equation (a, b) ∈ E ∪ E -1 can be replaced by a conversion a b.

Theorem 19. For all patterns l, sets of equations E, and terms T, r in Λ, if l ↑ T [E] and l,E r : T , then l → r preserves typing in Λ.

Proof. Let ∆ be an environment, σ be a substitution and A be a term of Λ such that ∆ lσ : A. By Lemma 15 (completeness), there are a term T , a set of equations E and a substitution σ extending σ such that t ↑ T [E ], σ |= E , ∆ σ : ∆ t and A T σ. Since ↑ is a function, we have T = T and E = E.

We now prove that, if Γ l,E t : T , then ∆, Γ σ t σ : T σ, by induction on l,E (note that σ replaces function symbols by terms).

(fun)

l,E Θ f : Σ f l,E f : Θ f . By induction hypothesis, we have ∆ Θ f σ :

Σ f σ = Σ f . Case f ∈ F. Then, f σ = f , Θ f σ = Θ f = Θ f and Σ f = Σ f .
By inverting typing rules, we get Θ f : Σ f . Therefore, by (fun) and (weak), ∆ f Hence, we have ∆ r σ : T σ. Since FV(r) ⊆ FV(l), we have r σ = rσ. Since ∆ lσ : A, by Lemma 6, either ∆ A : s for some sort s, or A = and lσ is of the form Πx 1 : A 1 , . . . , Πx k : A k , . Since l is a pattern, l is of the form f l 1 . . . l n . Therefore, ∆ A : s and, by (conv), ∆ rσ : A.

: Θ f , that is, ∆ f σ : Θ f σ. Case f = x ∈ FV(l). Then, f σ = xσ and Θ f σ = x σ. Therefore, ∆ f σ : Θ f σ since ∆ x σ : x σ. Case f = x with x ∈ FV(l). Then, f σ = x σ and Θ f σ = . Therefore, ∆ f σ : Θ f σ since ∆ x σ : . (conv) Γ l,E t : T T l,E U Γ l,E U : s Γ l,

F. Blanqui 11:13

Then, any total order on function symbols allows to prove preservation of typing. For instance, by taking f > → > a > a > b > b, we get the rewriting rules f → τ a → τ b, τ a → τ a, τ b → τ b and x → τ a, so that one can easily check that, modulo these rewriting rules, f x has type τ b. Therefore, app a b (lam a b f ) x → f x preserves typing.

Note that the result does not depend on the total order taken on function symbols. For instance, if one takes f > → > a > a > b > b (flipping the order of a and a ), we get the rewriting rules f → τ a → τ b, τ a → τ a , τ b → τ b and x → τ a . In this case, f x has type τ b as well. Flipping the order of b and b would work as well.

Related works and conclusion

The problem of type safety of rewriting rules in dependent type theory modulo rewriting has been first studied for simply-typed function symbols by Barbanera, Fernández and Geuvers in [START_REF] Barbanera | Modularity of strong normalization in the algebraic-λ-cube[END_REF]. In [START_REF] Blanqui | Definitions by rewriting in the calculus of constructions[END_REF], the author extended these results to polymorphically and dependently typed function symbols, and showed that rule left-hand sides do not need to be typable for rewriting to preserve typing. This was later studied in more details and implemented in Dedukti by Saillard [START_REF] Saillard | Type checking in the Lambda-Pi-calculus modulo: theory and practice[END_REF]. In this approach, one first extracts a substitution ρ (called a pre-solution in Saillard's work) from the typability constraints of the left-hand side l and check that, if l is of type A, then the right-hand side r is of type Aρ (in the same system). For instance, from the simplified set of constraints E = { x = T , p = N , v = V p, n = N , p = n} of our running example, one can extract the substitution ρ = {(n, p)} and check that v has type (V n)ρ = V p. However, it is not said how to compute useful pre-solutions (note that we can always take the identity as pre-solution). In practice, the pre-solution is often given by the user thanks to annotations in rules. A similar mechanism called inaccessible or "dot" patterns exists in Agda too [START_REF] Norell | Towards a practical programming language based on dependent type theory[END_REF].

An inconvenience of this approach is that, in some cases, no useful pre-solution can be extracted. For instance, if, in the previous example, we take the original set of constraints E = { x = T , p = N , v = V p, n = N , V (sp) = V (sn)} instead of its simplified version E , then we cannot extract any useful pre-solution.

In this paper, we proposed a more general approach where we check that the right-hand side has the same type as the left-hand side modulo the equational theory generated by the typability constraints of the left-hand side seen as closed equations (Theorem 19). A prototype implementation is available on: https://github.com/wujuihsuan2016/lambdapi/tree/sr.

To ensure the decidability of type-checking in this extended system, we propose to replace these equations by an equivalent but convergent rewriting system using Knuth-Bendix completion [15] (which always terminates on closed equations), and provide conditions for preserving the termination and confluence of the system when adding these new rules (Lemma 27). This approach has also the advantage that Dedukti itself can be used to check the type safety of user-defined Dedukti rules.

We also showed that, for the algorithm to work, the typability constraints sometimes need to be simplified first, using the fact that some function symbols are injective (Theorem 22). It would be interesting to be able to detect or check injectivity automatically (see [START_REF] Wu | Checking the type safety of rewrite rules in the λπ-calculus modulo rewriting[END_REF] for preliminary results on this topic), and also to find a simplification procedure more general than the one of Figure 3.

  constant symbol T : TYPE // Dedukti type for representing simple types constant symbol arr : T → T → T // arrow simple type constructor injective symbol τ : T → TYPE // interprets T elements as Dedukti types rule τ ( arr $x $y ) → τ $x → τ $y // ( Curry -Howard isomorphism ) // representation of simply -typed λterms symbol lam : Π a b , (τ a → τ b ) → τ ( arr a b ) symbol app : Π a b , τ ( arr a b

FFigure 1

 1 Figure 1 Typing rules of the λΠ-calculus modulo rewriting

Example 7 .

 7 (*) there exist ∆ and B such that ∆ l : B and ∆ r : B. But this condition is not sufficient in general as shown by the following example: Consider the rule f (xy) → y with f : B → B. In the environment ∆ = x : B → B, y : B, we have ∆ l : B and ∆ r : B. However, in the environment Γ = x : A → B, y : A, we have Γ l : B and Γ r : A.

Figure 2

 2 Figure 2 Typability constraints

F

Definition 17 .

 17 Given a pattern l and a set of equationsE such that R contains no variable of {x | x ∈ FV(l)} ∪ { x | x ∈ FV(l)} 3 , we define a new λΠ-calculus Λ l,E = (F , V , Θ , Σ , R )where:

https://github.com/Deducteam/lambdapi

This can always be done by renaming variables.

Note that, here, we extend the notion of substitution by taking maps on V ∪ F .

https://github.com/Deducteam/lambdapi

Example 20. We have seen that cons x p v ↑ V (sp)[E 1 ] with E 1 = { x = T , p = N , v = V p}, and tail n (cons x p v) ↑ V n [E 2 ] with E 2 = E 1 ∪ { n = N , V (sp) = V (sn)}. After the previous theorem, the rewriting rule defining tail preserves typing if we can prove that l,E2 v : V n where, in Λ l,E2 , v is a function symbol of type v and sort , and types are identified modulo and the equations of E 2 . But this is not possible since v : V p and V p l,E2 V n. Yet, if σ |= V (sp) = V (sn) and V and s are undefined then, by confluence, σ |= p = n and thus σ |= V p = V n. We therefore need to simplify the set of equations before type-checking the right-hand side.

Simplification of typability constraints

In this section, we show that Theorem 19 can be generalized by using any valid simplification relation, and give an example of such a relation. We have seen in the previous example that, thanks to confluence, σ |= p = n whenever σ |= sp = sn and s is undefined. But this last condition is a particular case of a more general property:

Definition 21 (Valid simplification relation). A relation on sets of equations is valid if, for all sets of equations

Definition 23 (I-injectivity). Given f : Πx 1 : T 1 , . . . , Πx n : T n , U and a set I ⊆ {1, . . . , n}, we say that f is I-injective when, for all t 1 , u 1 , . . . , t n , u n , if f t 1 . . . t n f u 1 . . . u n and, for all i / ∈ I, t i u i , then, for all i ∈ I, t i u i .

For instance, f is {1, . . . , n}-injective if f is undefined. The new version of Dedukti allows users to declare if a function symbol is I-injective (like the function τ in Example 1), and a procedure for checking I-injectivity of function symbols defined by rewriting rules has been developed and implemented in Dedukti [START_REF] Wu | Checking the type safety of rewrite rules in the λπ-calculus modulo rewriting[END_REF]. For instance, the function symbol τ of Example 1, which is defined by the rule τ (arr x y) → τ x → τ y, can be proved to be {1}-injective.

Clearly, I-injectivity can be used to define a valid simplification relation. In fact, one can easily check that the following simplification rules are valid too: Lemma 24. The relation defined in Figure 3 is a valid simplification relation.

Proof. We only detail the first rule which says that, if some substitution σ validates some equation t = u, that is, if tσ uσ, then σ validates any equation t = u where t and u are reducts of t and u respectively. Indeed, since t is a reduct of t, t t. Similarly, u u. Therefore, by stability of conversion by substitution and transitivity, t σ u σ.

11:11

Example 25. We can now handle our running example. We have

The above simplification relation works for the rewriting rule defining tail but may not be sufficient in more general situations:

Example 26. Let D be the set of equations {f ct = ga, f cu = gb, a = b} and assume that f is {2}-injective. Then the equation t = u holds as well, but D cannot be simplified by the above rules because it contains no equation of the form f ct = f cu.

We leave for future work the development of more general simplification relations.

Decidability conditions

We now discuss the decidability of type-checking in Λ l,E and of the simplification relation based on injectivity, assuming that → β ∪ → R is terminating and confluent so that typechecking is decidable in Λ. In both cases, we have to decide l,E , the reflexive, symmetric and transitive closure of

, where E is a set of closed equations.

As it is well known, an equational theory is decidable if there exists a convergent (i.e. terminating and confluent) rewriting system having the same equational theory: to decide whether two terms are equivalent, it suffices to check that their normal forms are identical.

In [15], Knuth and Bendix introduced a procedure to compute a convergent rewriting system included in some termination ordering, when equations are algebraic. Interestingly, this procedure always terminates when equations are closed, if one takes a termination ordering that is total on closed terms like the lexicographic path ordering > lpo wrt any total order > on function symbols (for more details, see for instance [START_REF] Baader | Term rewriting and all that[END_REF]).

For the sake of self-contentness, we recall in Figure 4 a rule-based definition of closed completion. These rules operate on a pair (E, D) made of a set of equations E and a set of rules D. Starting from (E, ∅), completion consists in applying these rules as long as possible. This process necessarily ends on (∅, D) where D is terminating (because D ⊆ > lpo ) and confluent (because it has no critical pairs).

We leave for future work the extension of this procedure to the case of non-algebraic, and possibly higher-order, equations.

If we apply this procedure to the set E of equations (assuming that they are algebraic), we get that l,E is the reflexive, symmetric and transitive closure of → β ∪ → R ∪ → D , where → β ∪ → R and → D are both terminating and confluent. However, termination is not modular in general, even when combining two systems having no symbols in common [START_REF] Toyama | Counterexamples to termination for the direct sum of term rewriting systems[END_REF].

There exists many results on the modularity of confluence and termination of first-order rewriting systems when these systems have no symbols in common, or just undefined symbols F S C D 2 0 2 0 11:12 Type safety of rewrite rules in dependent types (see for instance [START_REF] Gramlich | Modularity in term rewriting revisited[END_REF] for some survey). But, here, we have higher-order rewriting rules that may share defined symbols.

So, instead, we may try to apply general modularity results on abstract relations [START_REF] Doornbos | On the union of well-founded relations[END_REF]. In particular, for all terminating relations P and Q, P ∪ Q terminates if P steps can be postponed, that is, if P Q ⊆ QP * . In our case, we may try to postpone the D steps: Lemma 27. For all sets of higher-order rewriting rules R and D, we have that Assume that t → D u at position p and u → βR v at position q. If p and q are disjoint, then these reductions can be trivially permuted: t → βR → D v. The case p ≤ q (p prefix of q) is not possible since D is closed (c) and no right-hand side of D is βR-reducible (d). So, we are left with the case q < p:

Case u → β v. The case p = q1 is not possible since no right-hand side of D is headed by an abstraction (d). So, t| q is of the form (λx : A, b)a and the D step is in A, b or a. Therefore, t → β → * D v. Case u → R v, that is, when u| q = lσ where l is a left-hand side of a rule of R. The case p = qs where s is a non-variable position of l is not possible because no non-variable subterm of a left-hand side of R unifies with a right-hand side of D (e). Therefore, since l is left-linear (b), t| q is of the form lθ for some substitution θ, and the D step occurs in some xθ. Hence, t → R → * D v.

Example 28. As we have already seen, the typability conditions of l = tail n (cons x p v) is the set of equations

as total order on function symbols, the Knuth-Bendix completion procedure yields with > lpo the rewriting system


D Bresch

V Calvez

E Grenier

P Vigneaux

Editors J.-F Gerbeau

S Enault

D Lombardi

P Poncet

M Thiriet

MUCUS DYNAMICS SUBJECT TO AIR AND WALL MOTION *

This study presents a numerical investigation of basic interactions between respiratory mucus motion, air circulation and epithelium ciliated cells vibration. One focuses on identification of meaningful rheological parameters, physiological and numerical simulation dimensioning. These preliminary results are crucial before the study of more general configurations of respiratory mucus motion. The numerical study presented in this work aims at providing a first numerical tool able to simulate the effects of mucus mobility and its ability to carry out pathogens or to deliver aerosol therapy to membrane wall cells. Momentum diffusion is identified as the dominant effect, as expected in this micrometer scale configuration, and its associate momentum diffusion operator is shown to be extremely stiff. Furthermore, epithelium vibration is shown to be much more efficient than air circulation for mucus propulsion.

Résumé. Ce travail présente une étude numérique des interactions élémentaires entre le mucus de l'appareil respiratoire, la circulation de l'air et le mouvement vibratoire des cellules épithéliales ciliées tapissant la membrane pulmonaire. On s'intéresse en particulier à l'identification des paramètres les plus importants parmi les données rhéologiques, le dimensionnement physiologique et la configuration numérique. Ces résultats préliminaires sont cruciaux afin d'envisager l'étude de cas plus généraux de transport du mucus respiratoire. Un tel outil numérique aura pour but de quantifier la mobilité du mucus et sa capacité à expulser les pathogènes ou à délivrer un aérosol à travers membrane. La diffusion a été identifiée comme étant le phénomène dominant, comme attendu à l'échelle du micromètre, et son opérateur de diffusion associé a été observé comme étant extrêmement raide. D'autre part, il a été montré que la vibration de l'épithélium est bien plus efficace que la circulation de l'air pour assurer le transport du mucus.

CONTEXT

The present work aims at modeling and simulating motion of respiratory mucus film in lungs. Understanding mechanisms involved in mucus mobility is a key point for improving airway clearance, crucial for patients with respiratory disease. Human respiratory ducts are coated by epithelial cells that yield a carpet of cilia vibrating in the airway surface fluid to propulse the mucus stratum with entrapped particles outside the ventilatory system. Mucus is a visco-elastic fluid secreted by the respiratory epithelium (see figure 1) that protects tracheo-bronchial tree mucosa from dehydration and traps inhaled particles (allergens, carcinogens, dust, micro-organisms, and inflammatory debris) that come into contact with it. Ciliary motions that are associated with mucus propulsion in the human respiratory tract cleans airways, as it flows from either the tracheo-bronchial tree or upper airways toward the pharynx, where it is swallowed (or expectorated). Strong mucus mobility is necessary to provide satisfying airway clearance and healthy behavior of respiratory system. Otherwise either proliferation of pathogen agents is observed in stagnant mucus, either lack of mucus film let pathogens slip through pulmonary membrane. Long time targeted diseases are cystic fibrosis (or mucoviscidosis), cilio-toxic affections (chronic bronchitis, virus and/or infections), and aerosol-therapy (lung cancer). Such pathologies impact directly fluid features, such as viscosity main modulus, as shown on figure 2.

Mucus delivery and motion are investigated in the framework of flow and transport in the respiratory tract that also include gas transport under several alternative modes of ventilation, aerosol transport, and surfactant dynamics in the terminal airways. Most -if not almost all -studies of transport in airways deal with gas convection in large respiratory ducts and diffusion in bronchioles and alveoli [START_REF] Baffico | Multiscale modelling of the respiratory tract[END_REF][START_REF] Thiriet | The effect on gas mixing of a He-O2 mixture in chronic obstructive lung diseases[END_REF]. Experimental and numerical tests focused on unsteady 1D flow in collapsible tracheas and bronchi during forced expiration, a usual manoeuvre daily used in lung function testing [START_REF] Thiriet | Experiments on flow limitation during forced expiration in a monoalveolar model[END_REF][START_REF] Thiriet | Experimental and theoretical models of flow during forced expiration. Pressure and pressure history dependence of flow rate[END_REF]. In parallel, some studies aimed at investigating tracheal deformation both experimentally [START_REF] Thiriet | Banc d'analyse optique de la déformation de conduites souples[END_REF] and numerically [START_REF] Begis | A finite element model of tracheal collapse[END_REF]. Results obtained from measurements and simulations were validated by in vivo observations. 3D flows in simple representative models of airways were afterward studied [START_REF] Conca | A computational model of Navier-Stokes equations with imposed pressure and velocity fluxes[END_REF][START_REF] Thiriet | A pulsatile developing flow in a bend[END_REF]. Nowadays, numerical simulations are carried out in fluid domains extracted from medical image processing [START_REF] Fetita | An image-based computational model of oscillatory flow in the proximal part of tracheobronchial trees[END_REF] and compared to in vivo measurements [START_REF] De Rochefort | In vitro validation of CFD simulation in human proximal airways reconstructed from medical images with hyperpolarized helium-3 MRI phase contrast velocimetry[END_REF]. However, due to impaired mucus-associated respiratory pathologies as well as development of minimally invasive aerosol-based therapeutic procedures, there is a need in knowledge and simulations of mucus delivery and motion in patient-specific domains. Many attempts have been made to study separately epithelium, mucus features, and air flow in lungs, but coupling of movements of these interacting compartments has never been investigated. The present manuscript aims at providing a two-dimensional preliminary study on numerical simulation dimensioning, scalability, and dominant effects. Such simulations are crucial in order to predict achievable size of future threedimensional simulations, and also to consider meaningful model for future work. Section 2 describes the fluid model, and section 3 gives problem dimensioning and boundary conditions. Discretization and numerical algorithms are described in section 4, and section 6 presents the results obtained, strategies and dominant effects.

MULTI-FLUID AND VORTICITY FORMULATION OF NAVIER-STOKES EQUATION

The methodology used in order to model incompressible multi-fluid, is to consider one fluid with variable viscosity and density with possible jumps of these quantities. This is valid with both diffuse and sharp interfaces, and with possible extension to miscible fluids. The Navier-Stokes equations in their velocity-vorticity formulation will be used. This allows simplifying pressure from these equations and thus avoiding a special treatment for continuity of normal component of stress tensor at interface. Furthermore, in a first approach, mucus is assumed to be a Newtonian fluid.

Lagrangian velocity-vorticity methods, in three-dimensional context, have been proved to be robust, efficient and linearly scalable in complex geometry [START_REF] Poncet | Analysis of an immersed boundary method for three-dimensional flows in vorticity formulation[END_REF] under complex boundary conditions [START_REF] Poncet | Analysis of direct three-dimensional parabolic panel methods[END_REF] with spatial [START_REF] Poncet | Spatially distributed control for optimal drag reduction in cylinder wakes[END_REF] or temporal [START_REF] Poncet | Spatially distributed control for optimal drag reduction in cylinder wakes[END_REF] strong variations, and even in multi-scale context [START_REF] Ossmani | Efficiency of multiscale hybrid grid-particle vortex methods[END_REF]. Consequently, these methods offer a good framework for future work and have been chosen for the present study, in a simplified context. In order to do so, the method developed in [START_REF] Coquerelle | A vortex level set method for the two-way coupling of an incompressible fluid with colliding rigid bodies[END_REF] for fluid solid interaction is extended in the present work to multi-fluid.

In the present work vorticity-velocity equations will be considered, but in their Eulerian formulation, in order to provide pure grid-based schemes and thus help implementation. Moreover, one can expect that the present microfluidic context makes the momentum diffusion stability condition much stronger than transport CFL condition, and consequently the use of Lagrangian method for this kind of simulation would be questionable.

One considers a two-dimensional domain of height L and width A, as displayed on figure 3. The Newton's second law of mechanics reads ∂ρu ∂t

+ ∇ • (ρu ⊗ u) = ∇ • σ + f ( 1
)
where u is the velocity field, ρ the volumic mass, σ = -pI + µ2D(u) + λ(∇ • u)I, µ the dynamic viscosity, p the pressure and g the gravity, and D(u) = (∇u + ∇u T )/2 is the strain-rate tensor. Moreover, the mass conservation reads ∂ t ρ + ∇ • (ρu) = 0. Now, two formalisms relying on vorticity will be derived from this equation.

For an incompressible fluid of constant density and viscosity, one gets the usual Navier-Stokes equations :

∂u ∂t + u • ∇u -ν∆u = (f -∇p)/ρ (2)
where ν = µ/ρ. Introducing vorticity ω = ∇ ∧ u and taking the curl of equation [START_REF] Baffico | Multiscale modelling of the respiratory tract[END_REF] gives

∂ω ∂t + u • ∇ω -ω • ∇u -ν∆ω = 0 (3)
since ∇ ∧ ∇p ≡ 0, under the condition that f derives from a potential (such as gravity). The compatibility relations are ∇ • u = 0 (incompressibility), ∇ • ω = 0 (by definition of vorticity). Physical boundary conditions related to the Navier-Stokes equations are often no-slip boundary conditions u = 0. When wall models are involved, non-homogeneous conditions can be considered and/or with first order conditions, such as Navier conditions for solid wall or normal stress for fluid-solid interaction and deformable walls. Virtual boundary conditions can also be periodic, no-slip-through or Neumann. In order to avoid boundary conditions on velocity while the equation is on vorticity, tangential components of boundary conditions on velocity can be converted in Neumann or Fourier-Robin conditions on vorticity [START_REF] Poncet | Analysis of direct three-dimensional parabolic panel methods[END_REF].

Furthermore, for an incompressible fluid of variable density and viscosity, on which the external force is only gravity f = ρg, equation [START_REF] Adams | MUDPACK: Multigrid Fortran Software for the Efficient Solution of Linear Elliptic Partial Differential Equations[END_REF] gives

∂u ∂t + u • ∇u - 1 ρ ∇ • (2µD(u)) = g -∇p/ρ (4)
one which continuity of normal stress σ • n at interface can be applied. Taking the curl of this last equation gives

∂ω ∂t + u • ∇ω -ω • ∇u -∇ ∧ 1 ρ ∇ • (2µD(u)) = 1 ρ 2 ∇ρ ∧ ∇p (5)
since ∇ ∧ ∇p = 0 and ∇ ∧ g = 0. Since ∇p is defined through equation ( 4), equation ( 5) can be written

∂ω ∂t + u • ∇ω -ω • ∇u -∇ ∧ 1 ρ ∇ • (2µD(u)) = 1 ρ ∇ρ ∧ g - ∂u ∂t -u • ∇u + 1 ρ ∇ • (2µD(u)) (6)
In Boussinesq approximation, the right hand side of equation can be reduced to (∇ρ ∧ g)/ρ (see [START_REF] Coquerelle | A vortex level set method for the two-way coupling of an incompressible fluid with colliding rigid bodies[END_REF] for instance). This equation is valid in 3D as well as in 2D. Moreover, in the present context of two-dimensional simulations, one has no stretching, that is to say ω • ∇u ≡ 0. Equation ( 6) is then finally reduced to

∂ω ∂t + u • ∇ω -∇ ∧ 1 ρ ∇ • (2µD(u)) = 1 ρ ∇ρ ∧ g (7)
The full governing equations, involving this last one, viscosity and density transport and boundary conditions are displayed in equation ( 15) of section 4. Transport is kept in the model since it can be non neglectible under some circumstances, such as strong velocities at boundaries or strong alteration of air/mucus interface topology [START_REF] Craster | Surfactant transport on mucus films[END_REF].

Furthermore, in case of jumps of µ and ρ, calculus has to be performed carefully since multiplication of Dirac and non-continuous functions may occur. Two arguments are required to have rigorous expression of diffusion terms. On the one hand, density and viscosity can jump only together (at the fluid interface) and are both affine function of a level-set 1. Mucus and air rheology parameters (assumed to be Newtonian).

(characteristic of species, see equation ( 17) of section 4), so ρ -1 ∇µ has the same regularity as µ -1 ∇µ = ∇ ln µ which are well defined as generalized function whatever the regularity of µ. On the other hand, the Boussinesq approximation makes the µD(u) with the same regularity as the stress σ, that is to say continuous when no surface tension is involved (the present work does not involve it in the governing equations). Consequently, ∇ • µD is a smooth generalized function (which may be singular, but without jump derivation), and ρ -1 ∇ • µD is well defined, even at interface. Moreover, if future work requires on non-Boussinesq hypothesis (such as short wavelengths generated at interface), equation [START_REF] Chorin | Numerical study of slightly viscous flow[END_REF], involving (u, ω, p) models in the spirit of [START_REF] Riaz | Vorticity interaction mechanisms in variable-viscosity heterogeneous miscible displacements with and without density contrast[END_REF], will has to be used instead of equation [START_REF] Conca | A computational model of Navier-Stokes equations with imposed pressure and velocity fluxes[END_REF]. If surface tension is also involved, then stress is no longer continuous and more general vorticity formulations will be required in order to satisfy mathematical definition of diffusion, such as ζ = ∇ × (ρu).

DIMENSIONAL AND BOUNDARY CONDITIONS SETUP

We consider a small box of height L and width A, as displayed on figure 3. As the dimension of the box in negligible compared to the size of the bronchiole (in the third generation L = 8µm has to be compared to the radius R = 3mm), we can reasonably forget the curvature of the bronchiole in our model, and define differential operators in Cartesian coordinates (and not in curvilinear). In the transverse direction, we suppose periodic boundary conditions, the domain size for transverse motion being chosen as A = 2.5 10 -5 m = 25µm (a few epithelial cells).

Rheology is provided by the data on table 1. In the future, variable rheology of mucus will be considered in order to take into account mucus hydration, mucin hydrophilic feature and exogen ratio. Currently, mucus rheology will be considered as constant in space and time.

We set a boundary condition on the velocity V ref , on the top virtual boundary (see figure 3) . Moreover, for this preliminary study, we did not consider the epithelium layer with the complex interaction between cilia, goblet cells, mucin and mucus. We use a virtual boundary at the bottom in the spirit of [START_REF] Craster | Surfactant transport on mucus films[END_REF] (cf. figure 3), and a boundary condition on the velocity, u cell , that model the output of the epithelium box.

Top reference velocity

Due to small size of bronchioles, air flow is supposed to be a Poiseuille flow in a cylindrical duct of radius R, satisfying full no-slip condition at r = R -H, H being the mucus height, with mucus assumed to have a velocity neglectful toward air velocity. Indeed, this configuration leads to a small Womersley number W o = R 2π/νT = 0.93, thus the frequency of pulsations is sufficiently low, so that the flow will be very nearly in phase with the pressure gradient, and will be given to a good approximation by Poiseuille's law. We can notice that a large radius or a faster respiratory cycle would lead to a state for which Poiseuille's law would not be applicable.

Maximum velocity at the center of bronchiole, assumed to be known and denoted V max (t), leads to the velocity V ref at the top of computational box, which is the velocity at location r = R -L (see figure 3). The reference velocity V ref at the top of computational box is then obtained by means of Poiseuille formula :

V ref (t) = v max (t) 1 - R -L R 2 ≈ 2L R V max (t) ≈ 4L πR 3 Q(t) (8)
where Q(t) is the flow rate, or equivalently V max (t) is the maximum velocity at the duct center. In the present configuration (typical size of third generation), V ref (t) (or V max (t) up to a constant) is the function displayed on figure 4. The interpolation methodology used to obtain this function V max (t) is developed in the following section 3.2.

Maximum air velocity during respiratory cycle

In order to model the reference velocity at the top virtual boundary, we need to find a function f (t) that mimic the air flow into the lung during the respiratory cycle. In fact, the air flow function is linked to the air velocity in the lumen and if the air velocity profile is supposed to follow a Poiseuille profile, the maximum air velocity can be found throught the following equation:

f (t) = S v • ndx = 2πv max (t) R-H 0 1 - r R 2 rdr ≈ π 2 R 2 v max (t) (9)
The respiratory cycle can be split in two phases: an inspiration one during T i and an expiration one during T e . In the following we suppose that our bronchal tree is perfect, that is to say the breathed volume V n I in the bronchial of the n-th generation is linked to the breathed volume in the trachea V I :

V n I = 2 -n V I
We have searched the air flow function in the form

f (t) = Aθ (t) sin(θ(t)) (10)
In order to fit the biologic phenomenon of breath, we impose the following constraints:

(1) At time t = 0, t = T i and t = T e +T i the air flow is zero. So, we want θ(0) = 0, θ(T i ) = π and θ(T i +T e ) = 2π.

(2) The breathed volume is given V I . That imposes the amplitude of the function: 2A = V I .

(3) We also want to control the derivatives at the begining d and end f of the respiratory cycle with d and f close to each other.

θ (0) = d A θ (T i + T e ) = f
A We interpolate the function θ by a polynomial function of degree 4:

θ(t) = αt 4 + βt 3 + γt 2 + δt
The previous conditions lead to a simple linear system:

    0 0 0 1 T 4 i T 3 i T 2 i T i T 4 T 3 T 2 T 4T 3 3T 2 2T 1         α β γ δ     =       2 d V I π 2π 2 f V I      
The numerical input parameters are T = 4s, T e /T i = 1.4 thus T i = 1.667s and T e = 2.333s, V I = 500mL, which leads to the function displayed on figure 4.

Wall condition and cilia model

A cilium has an planar asymmetrical motion:

• propulsion phase during a time T prop .

• rest phase during a time T rest .

• recovery phase during a time T return : the cilium bends while its return. During this time, the cilium has a mean velocity u m . We assume that during the propulsion phase, the cilium has a sinusoidal movement, and after the velocity is zero during T rest + T return .

u cell (t) = Ttotπum 2Tprop sin( πt Tprop ) sur t ∈ [0, T prop ] 0 sur t ∈ [T prop , T tot ] (11)
This leads to the curve displayed on figure 5 for cilium velocity, using values displayed on table 2. In our simulation we use a box of 10 cilia. A model for the motion of a set of cilia is required, and we have decided to take a traveling wave of our cilium motion, so that this "kind of hola" motion could be optimal for mucus transport. Future work will require to investigate effect of cilia synchronization of mucus transport, and to build a more realistic model based on these investigations. In practice, cilia have a periodic motion with a frequency that lies between 4 and 25 Hz; The value 20 Hz has been selected for our simulations.

Non-dimensional system

In order to improve a priori conditioning of the discretization of problem [START_REF] Conca | A computational model of Navier-Stokes equations with imposed pressure and velocity fluxes[END_REF], the unknowns and the parameters are rescaled in the following way: where x * and u * are the standard distance and velocity of the problem.

x = x x * u = u u * t = t t * = tu * x * (12)
ρ = ρ ρ * µ = µ µ * (13)
The Froud number, the dynamical pressure and the Reynolds number are introduced:

F r = u * √ gx * , p = p ρ * u * 2 , Re = ρ * x * u * µ *
We get, as a rescaling of equation ( 4) :

∂ t u + u • ∇u = - 1 ρ ∇p + êg F r 2 + 1 Reρ ∇ • (µ2D(u)) ( 14
)
where êg is the unit vector of the gravitation. With these data, we found Re ≈ 4 • 10 -6 and F r 2 ≈ 1, 27 • 10 -2 . It seems, as expected in most of microfluidic configurations, that the diffusion term is the most important phenomenon in the model.

NUMERICAL ALGORITHM

One considers equation ( 7) with boundary conditions described in the last section :

                                     ∂ω ∂t + u • ∇ω -∇ ∧ ρ -1 ∇ • (2µD(u)) = ∇ ln ρ ∧ g in Ω×]0, T [ ∇ • u = 0, ω = ∇ ∧ u in Ω×]0, T [ ∂ρ ∂t + u • ∇ρ = 0 and ∂µ ∂t + u • ∇µ = 0 in Ω×]0, T [ u(x, 0) = u 0 (x) on Ω ρ(x,
with aftermath non-dimensionalization. Initial mucus location is a film of height H = L/3 at the bottom of the computational box. First, we introduce a vorticity-to-velocity operator in the spirit of [START_REF] Poncet | Analysis of an immersed boundary method for three-dimensional flows in vorticity formulation[END_REF] : on a domain Ω =]0, A[×]0, L[, one takes at each time a vorticity field ω : Ω → R and gets a velocity field u

: Ω → R 2 .
This linear application is denoted u = Aω, and it is defined by means of the two following embedded Poisson equations :

                               -∆Ψ = ∇ ∧ u = ω on ]0, A[×]0, L[ ∆Φ = 0 on ]0, A[×]0, L[ Ψ(0, x 2 ) = Ψ(A, x 2 ) x 2 ∈ [0, L] Φ(0, x 2 ) = Φ(A, x 2 ) x 2 ∈ [0, L] ∂ x2 Ψ + ∂ x1 Φ = u • τ = V ref on x 2 = L ∂ x2 Ψ + ∂ x1 Φ = u • τ = u cell on x 2 = 0 -∂ x1 Ψ + ∂ x2 Φ = u • n = 0 on x 2 = 0 (16)
whose solutions satisfy Ψ, Φ : Ω → R, and the velocity is finally computed by u = ∇ ∧ Ψ + ∇Φ. For the Laplacian as well as for the derivations, centered second order finite difference schemes are used, except for fluxes at boundaries for which second order one-sided schemes are used.

Using such a vorticity-to-velocity operation A, one can use time splitting algorithm in order to split transport, gravity and momentum diffusion in equation [START_REF] Poncet | Analysis of an immersed boundary method for three-dimensional flows in vorticity formulation[END_REF]. These three steps, for which initial condition is the final value of the step performed before, are as follows :

(1) Convection -Let φ be a level-set function characteristic of mucus (that is to say of value 1 in mucus and 0 otherwise). The following transport equations are solved using centered second order finite differences in space and second-order Runge-Kutta scheme in time :

     ∂ω ∂t + Aω • ∇ω = 0 ∂φ ∂t + Aω • ∇φ = 0 (17)
both in Ω×]0, T [. CFL condition δt δx √ 3/ u ∞ of this scheme is non restrictive compared to the one induced by diffusion. Density ρ and dynamic viscosity µ are built afterward by ρ(x) = ρ mucus φ(x) + ρ air (1 -φ(x)) and µ(x) = µ mucus φ(x) + µ air (1 -φ(x)). Input/output boundary conditions on φ and ω are periodic for left and right domain boundary, and Dirichlet conditions and imposed at top and bottom boundaries.

(2) Gravity -The following dynamical system can be integrated analytically in time and be discretized by centered second order finite difference scheme :

∂ω ∂t = ∇ ln ρ ∧ g in Ω×]0, T [ (18)
whose source is localized at fluid interface. (3) Momentum diffusion -Over a time step [t n , t n + δt], the following momentum diffusion equation has to be solved :

∂ω ∂t -∇ ∧ ρ -1 ∇ • (2µD(u)) = 0 in Ω×]t n , t n + δt[ (19)
with initial condition ω(•, t n ) given as the final value of gravity step. In order to avoid three successive numerical differentiations, one prefers to introduce a Stokes equation whose initial condition is u(•,

t n ) = Aω(•, t n ) : ∂u ∂t - 1 ρ ∇ • (2µD(u)) + ∇γ = 0 in Ω×]t n , t n + δt[ (20)
since ρ is stationary over this sub-step (thus ∂ρ/∂t ≡ 0), with u = V ref , u cell or periodic boundary conditions. Relation ω(•, t n + δt) = ∇ ∧ u(•, t n + δt) gives the final vorticity. In the momentum diffusion step, γ is not the pressure, and it is not used to ensure boundary conditions on u, since most of them are already satisfied by definition of A : This is only as a projector on divergence free fields, and periodic and homogeneous Neumann boundary conditions are set at boundaries. Moreover, stability condition of the momentum diffusion equation reads

Max µ mucus ρ mucus , µ air ρ air δt δx 2 1 2 (21)
With the rheological data considered herein (see table 1), one gets (as expected) a tiny time step δt of around 10 -10 when Euler or RK2 time integrators are used (Explicit Euler scheme gives a CFL at δt 5 10 -10 ). An implicit scheme is then required. In order to solve the implicit momentum diffusion step, the GMRES method is used to solve the resultant linear system (and not a fast solver for standard Poisson equation since diffusion coefficients are not constant). Simulations have shown that the time of convergence of the GMRES algorithm is non linear with respect to the time step. In our configurations, it seems optimal to take the time step δt around 10 -4 .

Furthermore, in order to check consistency and convergence, the cilium velocity is chosen independent of x. As the initial condition is also independent of x, a first assumption is to suppose that all the functions are independent of x. One gets then ∂ x2 (µ∂ x2 u x1 ) = 0 and this equation can be easily solved analytically:

u x (y) = u cell + C µmucus y pour y ∈ [0, H] u x (y) = v ref + C µair (y -L) pour y ∈ [H, L] (22)
with

C = (v ref -u cell ) H µ mucus + L -H µ air -1
5. VALIDATION

Non uniform decelerating Couette flow

In this section a test case is set up in order to validate the scheme presented above. In particular the dynamics of a plane shear flow is taken into account. The domain is The convergence properties of the algorithm were analyzed. In particular the L 2 norm of the error is investigated as the space resolution is varied. The error is defined as follows:

[0, 4] × [0, 1] ⊂ R 2
E = 1 0 e(t) dt = 1 0 Ω (U x -U (e) ) 2 dΩ 1/2 dt, (30)
where U x is the numerical solution for the x component of the velocity field, U (e) is the exact solution.

In Fig. 7.a the time evolution of the L 2 space error is plotted as function of time for different resolutions. It has a maximum at approximately one fith of the simulated evolution, and it is quite small in norm. In Fig. 7.b the logarithmic plot of the error E is plotted as function of the space resolution. The scheme is convergent and the convergence rate shown by the experiments results a bit faster than the first order one. Several numerical experiments were done varying both the space and the time resolution, using the same linear algebra as in last section. In Fig. 8.a a longitudinal section (x = 0.5) of the exact solution is shown at different times. Let us note that whether a proper exponent is chosen for the upper bound velocity deceleration, the solution is the first Legendre polynomial, scaled, i.e. it coincides with the profile of the stationary Couette flow. This is due to the fact that the momentum change induced by the deceleration is compensated by the effect of the varying viscosity.

In Fig. 8.b a zoom of the solution at time T=0.3 is shown, when the space resolution is N x = 200, N y = 200. The agreement between the solution is good, no relevant error appeared in the numerical tests. We analyzed the convergence properties, as in the previous test case. In particular the L 2 norm of the error is investigated as the space resolution is varied. In Fig. 9.a, the time evolution of the L 2 space error is plotted as function of time for different resolutions. As in the previous case it shows a maximum at one fith of the time history. In Fig. 9.b, the logarithmic plot of the error E is shown. In the case of non-uniform diffusion the convergence rate is not exactly of the first order.

RESULTS

With the algorithm described in the previous section, the parameters and boundary conditions given in section 3, one gets that gravity and convection contribution are less than 1% : the dominant term and effect is the momentum diffusion.

The figure 10 is a result obtain with the first scheme. In the upper part, the velocity has a Poiseuille profile as expected. In the lower part (under the red dashed line), the shape of the velocity field is more complex: the cilium velocity field produces a vortex. The figure 10 also shows the shape of the vorticity in our numerical box : vorticity produced by the variation of cilium velocity u cell (x), is diffused in the mucus almost instantly, while it propagates through the interface to the air.

At the chosen time step δt = 10 -4 , 500 time steps are required for a 20Hz cilia cycle, and 4 10 4 time steps are required for a full 4s respiratory cycle. The typical spatial resolution is δx 1 = 1.25µm, and δx 2 = 100nm. These setting looks like the minimum requirements for having stable and physical solutions. This will be essential for dimensioning future three-dimensional simulations.

Furthermore, a crucial point for future work is to identify the stiffness of the problem, in order to apply adequate numerical algorithms. Since it has been shown above that momentum diffusion is the dominant effect, focus is given on the stiffness of operator :

H ρ,µ (u) = - 1 ρ ∇ • 2µ ∇u + ∇u T 2 ( 39
)
whose spectrum is displayed on figure 11 with several zoom levels. It depends on functions ρ and µ, and thus depends on air/mucus interface topology. The ratio between the largest and smallest eigenvalue modules is 5 10 9 , which means that H ρ,µ is an extremely stiff momentum diffusion operator. In future, state-of-the-art "stiff" integrators will be used, such as exponential integrators and/or Krylov-based time integrators [START_REF] Tokman | Efficient Integration of Large Stiff Systems of ODEs with Exponential Propagation Iterative (EPI) Methods[END_REF][START_REF] Tokman | Three-dimensional Model of the Structure and Evolution of the Coronal Mass Ejections[END_REF]. Another noticeable point is that air/mucus interface remains visibly flat (deviation amplitude can be neglected toward the film height), despite the strong vorticity generation. This phenomenon would be even stronger if surface tension was involved in the model. In the present study, standard (but stable) Eulerian schemes have been chosen for solving convection and air/mucus interface transport. Such schemes have been chosen more for convenience than for their ability to track these quantities accurately. Nevertheless, the use of Lagrangian methods would help to track more accurately fluid interface and would improve accuracy of transport terms (for both vorticity and level-set defining species). On this kind of problems, motivation for Lagrangian methods is not due to stability improvement such as in aerodynamic flow simulation at higher Reynolds numbers [START_REF] Poncet | Topological aspects of the three-dimensional wake behind rotary oscillating circular cylinder[END_REF]. Moreover, mucus motion follows the velocity induced by epithelium, and it is not significantly affected by air velocity (the ratio between their Reynolds numbers is 50, see table 1 for instance).

CONCLUSION AND PERSPECTIVES

The present study has brought precious preliminary results on dimensioning and strategies for future three-dimensional computation of mucus motion in tracheo-bonchial tree. Essential rheological parameters and characteristic lengths have been used in order to carry out minimum requirements on simulation size.

As expected, like in most microfluidic simulations, momentum diffusion is the dominant effect. The momentum diffusion operator as been shown to be extremely stiff. Adequate time integrators, such as exponential integrators, will be used in the future in order to minimize computational time of momentum diffusion.

Vorticity formulation allows to avoid management of continuity of normal stress at interface (which has to be satisfied if pressure is computed) and to consider the mucus, the pericelial mucus and the air as one fluid of sharply variable rheology. Computing accurately transport of species is then crucial in order to localize correctly mucus. In the future, more elaborate Lagrangian methods will be used for transport computing, despite the convection is not the dominant effect in global mucus motion.

The dominant input parameter has been shown to be the epithelium cilia displacement. Accurate computation of epithelium dynamics will be required in order to get correct bottom boundary conditions for mucus. Mechanical properties of cilia will also be investigated.

Moreover, the rheological model of mucus, considered as Newtonian in the present study for convenience, has to be improved in order to take into account visco-elasticity and hydration level. Visco-plasticity feature of mucus is also an open question that will also deserve attention in further works.

FIGURE 1 .

 1 FIGURE 1. Real lung geometry and membrane neighborhood with air/mucus/epithelium configuration.

FIGURE 2 .

 2 FIGURE 2. Different lung mucus viscosity values, altered by different pathologies, and modifying airway clearance feature (data from [13]).

FIGURE 3 .

 3 FIGURE 3. Setting of the model : Computational box and related boundary conditions (to the left), and Poiseuille flow setup for reference velocity at the top of the box (to the right).

FIGURE 4 .

 4 FIGURE 4. Reference velocity V ref during the respiratory cycle.

FIGURE 5 .

 5 FIGURE 5. Cilium velocity with respect to time (left picture) and naive model for velocity of 10 cilia with respect to time and space (right picture).

  FIGURE 6. a) Exact solution at different times, b) Exact and numerical solution at T=0.3, zoom.

FIGURE 8

 8 FIGURE 8. a) Exact solution at different times, b) Exact and numerical solution at T=0.2, zoom.

FIGURE 10 .

 10 FIGURE 10. Shape of velocity with u cell (x) as bottom boundary condition (left picture), and resulting vorticity (right picture). Dashed line is the interface air/mucus.

FIGURE 11 .

 11 FIGURE 11. Spectrum of momentum diffusion operator, quantifying stiffness induced by momentum diffusion, with successive zoom level from top left picture (a) displaying whole spectrum, to (d) displaying zero neighborhood. Picture (e) shows eigenvalue modules.

TABLE 2 .

 2 Characteristic time, space and velocity values for epithelium cilia motion.

  0) and µ(x, 0) constant by pieces with respect to values of table1

u(x, t) = V ref (x, t) on top boundary u(x, t) = u cell (x, t) on bottom boundary u(x, t) is periodic on lef t and right boundaries

* This project NUMULUNG has been funded by INSA Toulouse, University Pierre and Marie Curie (Paris VI) and INRIA (Project REO).

constant density, instead, the latter is supposed to be a function of the space coordinate. The flow being incompressible, the density is purely advected so that the mass conservation reduces to:

where v = (u, v) is the velocity field. Several hypothesis will be assumed in the following:

(1) the flow is invariant in x direction, i.e. ∂ x • = 0;

(2) u(x, 0, t) = 0; u(x,

The above assumptions holding, let ρ = ρ(y, t) be the density distribution, the conservation of momentum can be simplified and a simple diffusion problem is found:

Let us observe that, since the equation for the divergence is implicitly satisfied (v = 0 and u = u(y, t)) and the boundary conditions are sufficient to provide the motion, no pressure gradient is needed and, as consequence, the pressure can be set to p = 0 (or whatever constant value). Moreover, since the motion is purely parallel to the x direction, applying the conservation equation lead to:

We look for a separate solution so that:

where the partial derivative with respect to time is denoted with a dot while the partial derivative with respect to the space direction is denoted with a prime. Let us introduce the separation constant λ. The problem can be splitted into two sub-problems:

The first problem is a classical initial value problem and its solution is the exponential, λ being the inverse of a charateristic time scale. The second problem is a second order Sturm-Liouville boundary value problem. Several analytical solutions can be found for different expressions of ρ(y). Let us recall a classical result of mathematical physics:

where

, the following solution is obtained:

(29) If the initial condition for the problem is u(y, 0) = y exp(1/2(1 -y 2 )) the solution written above is the unique solution for the problem.

Several numerical experiments were done varying both the space and the time resolution. A krilov subspace solver (GMRES) and a biconjugate gradient stabilized solver were tried for the solution of the diffusion substep, with analogous results. In Fig. 6.a) a longitudinal section (x = 0.5) of the exact solution is shown at different times. Let us note that the stationary solution for a variable density Couette flow is a linear profile, as in the constant density case: the non-stationary boundary conditions induce a space variation in the inertia term, making the profile curve. In Fig. 6.b) a zoom of the solution at time T=0.3 is shown, when the space resolution is N x = 200, N y = 200. The agreement between the solution is good, no relevant error appeared in the numerical tests.

Space varying viscosity decelerating Couette flow

In this section a test case is set up in order to validate the scheme presented above. In particular the dynamics of a space varying viscosity Couette flow is taken into account. This flow allows us to validate the non-linear diffusion operator.

The domain is

, and the time is t ∈ [0, 1]. The fluid is now considered of constant density, but the dynamic viscosity is supposed to be a function of the space coordinate. The mass conservation reduces to the incompressibility of the fluid:

where v = (u x , u y ) is the velocity field.

Several hypothesis will be assumed in the following:

(1) the flow is invariant in x direction, i.e. ∂ x • = 0;

(2) u x (x, 0, t) = 0; u x (x, 1, t) = f (t); ∀x (3) u y (x, 0, t) = 0; u y (x, 1, t) = 0. The equation for the divergence leads to ∂ y u y = 0, and with the boundary conditions, we obtain u y (y) = 0.

The above assumptions holding, let ν = ν(y) be the density distribution, the conservation of momentum can be simplified and a simple diffusion problem is found:

Let us observe that the equation for y gives ρg = ∂ y p (33) thus p = ρgy + cste. We search solution u x as a function of separate variables:

and then we obtain: θ(t)

that is to say

Now, we notice that the n-th Legendre polynomial L n is solution of :

When n is odd, L n is null at the origin. We can take ν(y) = 1 -y 2 , λ n = -(2n + 1)(2n + 2) < 0 and

We recall some Legendre polynomials of odd order:

• P 1 (x) = x,

• P 3 (x) = 1 2 (5x 3 -3x),

• P 5 (x) = 1 8 (63x 5 -70x 3 + 15x),

• P 7 (x) = 1 16 (429x 7 -693x 5 + 315x 3 -35x),

• P 9 (x) = 1 128 (12155x 9 -25740x 7 + 18018x 5 -4620x 3 + 315x).