E Ciurea

C Schiopu

The Minimum Flows in Bipartite Dynamic Networks. The Static Approach

Keywords: bipartite dynamic network flow, minimum flow, bipartite static network flow

In this paper we study minimum flow algorithms for bipartite dynamic networks. We resolve this problem by rephrasing into a problem in bipartite static network. In a bipartite static network the several minimum flow algorithms can be substantially improved. The basic idea in this improvement is a two arcs pull rule. In the final of the paper we present an example.

Introduction

The theory of flow is one of the most important parts of Combinatorial Optimization. The static network flow models arises in a number of combinatorial applications that on the surface might not appear to be optimal flow problems at all. The problem also arises directly in problems as far reaching as machine scheduling, the assignment of computer modules to computer processor, tanker scheduling etc. [1]. However, in some applications, the time is an essential ingredient [3], [START_REF] Ford | Flow in Networks[END_REF], [START_REF] Salehi | Minimum flow problem on network flows with time-varying bounds[END_REF], . In this case we need to use dynamic network flow model. On the other hand, the bipartite static network also arise in practical context such baseball elimination problem, network reliability testing etc. and hence it is of interes to find fast flow algorithms for this class of networks [START_REF] Ahuja | Improved algorithms for bipartite network flows[END_REF], [START_REF] Gusfield | Fast algorithms for bipartite network flow[END_REF].

The computation of a maximum flow in dynamic networks with lower bound function e = 0 has been investigated by many researchers [3] , [START_REF] Ford | Flow in Networks[END_REF], [START_REF] Skutella | An Introduction to Network Flows Over Time[END_REF], [START_REF] Wilkinson | An algorithm for universal maximal dynamic flows in network[END_REF] and other. Unfortunately, the minimum flow problem in dynamic networks was not treated as often as the maximum flow problem. In the paper [START_REF] Ciurea | An algorithm for minimal dynamic flow[END_REF] is presented an algorithm for minimum flow in stationary dynamic networks with lower bounds e(i, j) ≥ 0 and upper bounds q(i, j) = ∞. Salehi H.F., et al. [START_REF] Salehi | Minimum flow problem on network flows with time-varying bounds[END_REF] has considered minimum flow problem on a class of dynamic networks called dynamic generative networks, where travel and transmission time are instantaneous but flow generation and consumption are dynamical. The authors consider the lower bounds of the arcs leaving source node are increasing functions of the time, the lower bounds of the arcs entering the sink node are decreasing functions of the time and the lower bounds of the other arcs and upper bounds for all arcs are constant.

In this paper we present the minimum flow problem in bipartite dynamic networks. This problem has not been treated present so far. Further on, in Section 2 we discuss some basic notions and results for minimum flow problem in general static networks. Section 3 deals with minimum flow problem in general dynamic networks. In Section 4 we present algorithms for flow problems in bipartite static network with lower bounds and in Section 5 we discuss the minimum flow problem in bipartite dynamic networks with lower bounds. Section 6 deals with an example for problem presented in Section 5.

Terminology and preliminaries.

In this section we discuss some basic notations and results used in the rest of the paper.

Let G = (N, A, l, u) be a static network with the set of nodes N = {1, . . . , n}, the set of arcs A = {a 1 , . . . , a k , . . . , a m }, a k = (i, j), i, j ∈ N , the lower bound function l : A → N, the upper bound (capacity) function u : A → N, with N the natural number set, 1 the source node and n the sink node.

For a given pair of subset X, Y of the nodes set N of a network G we use the notation:

(X, Y ) = {(i, j)|(i, j) ∈ A, i ∈ X, j ∈ Y }
and for a given function f on arcs set A we use the notation:

f (X, Y ) = (X,Y ) f (i, j)
A flow is a function f : A → N satisfying the next conditions:

f (i, N ) -f (N, i) =    v, if i = 1 0, if i = 1, n -v, if i = n (1a) l(i, j) ≤ f (i, j) ≤ u(i, j), (i, j) ∈ A (1b)
for some v ≥ 0. We refer to v as the value of the flow f . The minimum flow problem is to determine a flow f for which v is minimum.

For minimum flow problem a preflow f is a function f : A → N satisfying the next conditions:

f (i, N ) -f (N, i) ≤ 0, i ∈ N -{1, n} (2a)
l(i, j) ≤ f (i, j) ≤ u(i, j), (i, j) ∈ A (2b)

For a preflow f the deficit of each node i ∈ N is

e(i) = f (i, N ) -f (N, i) (3)
and if e(i) < 0, i ∈ N -{1, n} then we say that node i is an active node. We refer to a node i with e(i) = 0 as balanced. A preflow f satisfying the condition e(i) = 0, i ∈ N -{1, n} is a flow. Thus, a flow is a particular case of preflow.

We further assume, without loss of generality, that if (i, j)

∈ A the (j, i) ∈ A (if (j, i) / ∈ A we consider that (j, i) ∈ A with l(j, i) = u(j, i) = 0). A cut is set of arcs [X, X] = (X, X) ∪ ( X, X), (X, X) = {(i, j)|(i, j) ∈ A, i ∈ X, j ∈ X}, X ⊂ N , X = N -X, ( X, X) = {(i, j)|(i, j) ∈ A, i ∈ X, j ∈ X}.
The set (X, X) denote the set of forward arcs of the cut, and the set ( X, X) denote the set of backward arcs of the cut. We refer to a cut [X, X] as a 1 -n cut if 1 ∈ X and n ∈ X. For the minimum flow problem, the capacity c[X, X] of a 1 -n cut is c[X, X] = l(X, X) -u( X, X). We refer to a 1 -n cut whose capacity is the maximum among all 1 -n cuts as a maximum cut.

Theorem 1. The value of the minimum flow from the source node 1 to the sink node n in a network G = (N, A, l, u) equals to the capacity of the maximum 1 -n cut.

For the minimum flow problem, the residual capacity r(i, j) of any arc (i, j) ∈ A with respect to a given flow (preflow) f is given by r(i, j) = c(j, i) -f (j, i) + f (i, j) -l(i, j). The residual network is Ĝ = (n, Â, r) with  = {(i, j)|(i, j) ∈ A, r(i, j) > 0}.

The minimum flow problem in a network G = (N, A, l, u) can be solved in two phases:

(1) estabilish a feasible flow if it exists;

(2) if exists a feasible flow, estabilish a minimum flow; The solution of first fase is presented in [1]. There are three approaches for solving the minimum flow problem:

(1)using decreasing path algorithms;

(2) using preflow algorithms;

(3) minmax algorithms. Any directed path from the source nodes 1 to the sink node n in the residual network Ĝ corresponds to a decreasing path in the original network G. Thus, the decreasing path algorithms for minimum flow problem correspond to the augmenting path algorithms for maximum flow problem.

Also, the preflow algorithms for minimum flow problem correspond to the preflow algorithms for maximum flow problem.

In the next presentation we assume familiarity with minimum flow algorithms, and we omit many details. The reader interested in further details is urged to consult the paper [START_REF] Ciurea | Sequential and parallel algorithms for minimum flows[END_REF].

Minimum flows in dynamic networks

Dynamic network models arise in many problem settings, including production distribution systems, economic planning, energy systems, traffic systems, and building evacuation systems [3], [START_REF] Ford | Flow in Networks[END_REF], [START_REF] Skutella | An Introduction to Network Flows Over Time[END_REF].

Let G = (N, A, l, u) be a static network with the node set N = {1, . . . , n}, the arc set A = {a 1 , . . . , a m }, the lower bound function l, the upper bound (capacity) function u, 1 the source node and n the sink node. Let N be the natural number set and let H = {0, 1, . . . , T } be the set of periods, where T is a finite time horizon, T ∈ N. Let use state the transit time function h : A × H → N the time lower bound function e : A × H → N, the time upper bound function q : A × H → N, e (i, j; t) ≤ q(i, j; t), for all (i, j) ∈ A and for all t ∈ H. The parameter h(i, j; t) is the transit time needed to traverse an arc (i, j). The parameters e(i, j; t) and q(i, j; t) represents the minimum and respective maximum amount of flow that can travel over arc (i, j) when the flow departs from i at time t and arrives at j at time θ = t + h(i, j; t).

The minimum dynamic flow problem for T time periods is to determine a flow function g : A×H → N, which should satisfy the following conditions in dynamic network D = (N, A, h, e, q) :

T t=0 (g(1, N ; t) - τ g(N, 1; τ )) = w (4a) g(i, N ; t) - τ g(N, i; τ ) = 0, i = 1, n, t ∈ H (4b) T t=0 (g(n, N ; t) - τ g(N, n; τ )) = -w (4c)
e(i, j; t) ≤ g(i, j; t) ≤ q(i, j; t), (i, j) ∈ A , t ∈ H (5)

min w, (6)
where τ = t -h(k, i; τ ), w = T t=0 v(t), v(t) is the flow value at time t and g(i, j; t) = 0 for all t ∈ {T -h(i, j; t) + 1, . . . , T }.

Obviously, the problem of finding a minimum flow in dynamic network D = (N, A, h, e, q) is more complex than the problem of finding a minimum flow in static network G = (N, A, l, u). Happily, this complication can be resolved by rephrasing the problem in dynamic network D into a problem in static network

R 1 = (V 1 , E 1 , l 1 , u 1 ) called the reduced expanded network.
The static expanded network of dynamic network D = (N, A, h, e, q) is the network R = (V, E, l, u) with V = {i t |i ∈ N, t ∈ H}, E = {(i t , j θ )|(i, j) ∈ A, t ∈ {0, 1, . . . , T -h(i, j; t)}, θ = t+h(i, j; t), θ ∈ H}, l(i t , j θ ) = e(i, j; t), u(i t , j θ ) = q(i, j; t), (i t , j θ ) ∈ E. The number of nodes in static expanded network R is n(T + 1) and number of arcs is limited by m(T + 1) -

A • h(i, j),
where • h(i, j) = min{h(i, j; 0), . . . , h(i, j; T )}. It is easy to see that any flow in dynamic network D from the source node 1 to the sink node n is equivalent to a flow in static expanded network R from the source nodes 1 0 , 1 1 , . . . , 1 T to the sink nodes n 0 , n 1 , . . . , n T and vice versa. We can further reduce the multiple source, multiple sink problem in static expanded network R to a single source, single sink problem by introducing a supersource node 0 and a supersink node n+1 constructing static super expanded network

R 2 = (V 2 , E 2 , l 2 , u 2 ), where V 2 = V ∪ {0, n + 1}, E 2 = E ∪ {(0, 1 t )|t ∈ H} ∪ {(n t , n + 1)|t ∈ H}, l 2 (i t , j θ ) = l(i t , j θ ), u 2 (i t , j θ ) = u(i t , j θ ), (i t , j θ ) ∈ E, l 2 (0, 1 t ) = l 2 (n t , n + 1) = 0, u 2 (0, 1 t ) = u 2 (n t , n + 1) = ∞, t ∈ H.
We construct the static reduced expanded network

R 1 = (V 1 , E 1 , l 1 , u 1 ) as follows. We define the function h 2 : E 2 -→ N, with h 2 (0, 1 t ) = h 2 (n t , n + 1) = 0, t ∈ H, h 2 (i t , j θ ) = h(i, j; t), (i t , j θ ) ∈ E. Let d 2 (0, i t )
be the length of the shortest path from the source node 0 to the node i t , and d 2 (i t , n + 1) the length of the shortest path from node i t to the sink node n + 1, with respect to h 2 in network R 2 . The computation of d 2 (0, i t ) and d 2 (i t , n + 1) for all i t ∈ V are performing by means of the usual shortest path algorithms. The network

R 1 = (V 1 , E 1 , l 1 , u 1 ) have V 1 = {0, n + 1} ∪ {i t |i t ∈ V, d 2 (0, i t ) + d 2 (i t , n + 1) ≤ T }, E 1 = {(0, 1 t )|d 2 (1 t , n + 1) ≤ T, t ∈ H}∪{(i t , j θ )|(i t , j θ ) ∈ E, d 2 (0, i t )+h 2 (i t , j θ )+d 2 (j θ , n+1) ≤ T }∪{(n t , n+1)|d 2 (0, n t ) ≤ T, t ∈ H} and l 1 , u 1 are restrictions of l 2 , u 2 at E 1 .
Now, we construct the static reduced expanded network R 1 = (V 1 , E 1 , l 1 , u 1 ) using the notion of dynamic shortest path. The dynamic shortest path problem is presented in [3]. Let d(1, i; t) be the length of the dynamic shortest path at time t from the source node 1 to the node i and d(i, n; t) the length of the dynamic shortest path at time t from the node i to the sink node n, with respect to h in dynamic network D. Let as consider

H i = {t|t ∈ H, d(1, i; t) ≤ t ≤ T -d(i, n; t)}, i ∈ N , and H i,j = {t|t ∈ H, d(1, i; t) ≤ t ≤ T -h(i, j; t) -d(j, n; θ)}, (i, j) ∈ A.
The multiple source, multiple sinks static reduced expanded network R 0 = (V 0 , E 0 , l 0 , u 0 ) have

V 0 = {i t |i ∈ N, t ∈ H i }, E 0 = {(i t , j θ )|(i, j) ∈ A, t ∈ H i,j }, l 0 (i t , j θ ) = e(i, j; t), u 0 (i t , j θ ) = u 1 (i, j; t), (i t , j θ ) ∈ E 0 . The static reduced expanded network R 1 = (V 1 , E 1 , l 1 , u 1 ) is constructed from network R 0 as follows: V 1 = V 0 ∪{0, n+1}, E 1 = E 0 ∪{(0, 1 t )|1 t ∈ V 0 }∪{(n t , n+1)|n t ∈ V 0 }, l 1 (0, 1 t ) = l 1 (n t , n + 1) = 0, u 1 (0, 1 t ) = u 1 (n t , n + 1) = ∞, 1 t , n t ∈ V 0 , l 1 (i t ; j θ ) = l 0 (i t , j θ ) and u 1 (i t , j θ ) = u 0 (i t , j θ ), (i t , j θ ) ∈ E 0 .
We remark the fact that the static reduced expanded network R 1 (R 0 ) is always a partial subnetwork of static super expanded network R 2 (R). In references [START_REF] Ciurea | An algorithm for minimal dynamic flow[END_REF], [START_REF] Ford | Flow in Networks[END_REF] it is shown that a dynamic flow for T periods in the dynamic network D with e = 0 is equivalent with a static flow in a static reduced expanded network R 1 . Since an item released from a node at a specific time does not return to the location at the same or an earlier time, the static networks R, R 2 , R 0 , R 1 cannot contain any circuit, and are therefore acyclic always.

In the most general dynamic model, the parameter h (i) = 1 is waiting time at node i, and the parameter e(i; t), q(i; t) are lower bound and upper bound for flow g(i; t) that can wait at node i from time t to t + 1. This most general dynamic model is not discussed in this paper.

The maximum flow problem for T time periods in dynamic network D formulated in conditions (4), ( 5), ( 6) is equivalent with the maximum flow problem in static reduced expanded network R 0 as follows:

f 0 (i t , V 0 ) -f 0 (V 0 , i t ) =    v t , if i t = 1 t , t ∈ H 1 0, if i t = 1 t , n t , t ∈ H 1 , t ∈ H n -v t , if i t = n t , t ∈ H n (7a) l 0 (i t , j θ ) ≤ f 0 (i t , j θ ) ≤ u 0 (i t , j θ ), (i t , j θ ) ∈ E 0 (8)
min

H 1 v t , (9)
where by convention i t = 0 for t = -1 and i t = n + 1 for t = T + 1.

In stationary case the dynamic distances d(1, i; t), d(i, n; t) become static distances d(1, i), d(i, n).

Minimum flows in bipartite static networks

In this section we consider that static network G = (N, A, l, u) is bipartite static network. A bipartite network has the node set N partitioned into two subsets N 1 and N 2 , so that for each arc (i, j) ∈ A,

either i ∈ N 1 and j ∈ N 2 or i ∈ N 2 and j ∈ N 1 . Let n 1 = |N 1 | and n 2 = |N 2 |.
For minimum flow problem, without any loss of generality, we assume that n 2 ≤ n 1 . We also assume that source node 1 belongs to N 2 (if the source node 1 belonged to N 1 , then we could create a new source node 1 ∈ N 2 , and we could add an arc (1 , 1) with l(1 , 1) = 0, u(1 , 1) = ∞). A bipartite network is called unbalanced if n 2 << n 1 and balanced otherwise.

F IF Opref low n 3 n 2 2 n n 2 m + n 3 2 Highestlabel n 2 √ m n 2 n √ m n 2 m Excessscaling nm + n 2 log ū n 2 m + n 2 n log ū n 2 m + n 2 2 log ū
Table 1: Several minimum flows algorithms

In reference [START_REF] Gusfield | Fast algorithms for bipartite network flow[END_REF] the authors show that the time bounds for several minimum flow algorithms automatically improve when the algorithms are applied without modification to unbalanced networks. A careful analysis of the running times of these algorithms reveals that the worst case bounds depend on the number of arcs in the longest node simple path in the network. We denote this length by L. For general network, L ≤ n -1 and for a bipartite network L ≤ 2n 2 + 1. Hence for unbalanced bipartite network L << n. Column 3 of Table 1 summarizes these improvements for several network flow algorithms.

The authors of references [START_REF] Ciurea | Improved algorithms for minimum flows in bipartite networks[END_REF] obtain further running time improvements by modifying the algorithms. This modification applies only to preflow pull algorithm. They call it the two arcs pull rule. According to this rule, always pull flow from a node in N 2 and pull flow on two arcs at a time, in a step called a bipull, so that no excess accumulates at nodes in N 2 . Column 4 of Table 1 summarizes the improvements obtained using this approach. The reader interested in further details is urged to consult the paper [START_REF] Ciurea | Improved algorithms for minimum flows in bipartite networks[END_REF].

Minimum flows in bipartite dynamic networks

In this Section the dynamic network D = (N, A, h, e, q) is bipartite.

We construct the static reduced expanded network R 0 = (V 0 , E 0 , l 0 , u 0 ) and we remark the fact that the network R 0 is a bipartite network with

V 0 = W 1 ∪ W 2 , W 1 = {i t |i ∈ N 1 , t ∈ H}, W 2 = {i t |i ∈ N 2 , t ∈ H}. Let w 1 , w 2 , ε 0 be w 1 = |W 1 | w 2 = |W 2 |, ε 0 = |E 0 |. If n 2 << n 1 then
obvious that w 2 << w 1 . In the static bipartite network R 0 we determine a minimum flow f 0 with a generalization of bipartite FIFO preflow algorithm.

We recall that the FIFO preflow algorithm might perform several saturating pulls followed either by a nonsaturating pull or relabeled operation. We refer to this sequence of operations as a node examination. The algorithm examines active nodes in the FIFO order. The algorithm maintains the list Q of active nodes as a queue. Consequently, the algorithm select a node j from the front of Q, performs pulls from this node, and adds newly active nodes to the rear of Q. The algorithm examines node j until either it becomes inactive or it is relabeled. In the latter case, we add node j to the rear of the queue Q. The algorithm terminates when the queue Q of active nodes is empty.

The modified version of FIFO preflow algorithm for minimum flow in bipartite network is called bipartite FIFO preflow algorithm. A bipull is a pull over two consecutive admisible arcs. It moves deficit from a node j θ ∈ W 2 to another node k τ ∈ W 2 . This approach means that the algorithm moves the flow over the back path D → D = (j θ , i t , k τ ), i t ∈ W 1 , and ensures that no node in W 1 ever has any deficit. A pull of α units from node j θ to node i t decreases both e(i t ) and r 0 (i t , j θ ) by α units and increases both e(j θ ) and r 0 (j θ , i t ) by α units.

We specify that maintain the arc list E - 0 (j θ ) = {(i t , j θ )|(i t , j θ ) ∈ E 0 }. We can arrange the arcs in these lists arbitrarily, but the order, once decided, remains unchanged throughout the algorithm. Each node j θ has a current arc, which is an arc in E - 0 (i t ) and is the next candidate for admissibility testing. Initialy, the current arc of node j θ is the first arc in E - 0 (j θ ). Whenever the algorithm attempts to find an admissible arc emanating from node j θ , it tests whether the node's current arc is admissible. If not, it designates the next arc in the arc list as the current arc. The algorithm repeats this process until either to finds admissible arc or reaches the end of the arc list.

The generalizate bipartite FIFO preflow (GBFIFOP) algorithm is presented in Figure 1.

1: ALGORITHM GBFIFOP; 2: BEGIN 3: PREPROCESS 4: while Q = ∅ do 5:
select the node j θ from the front of Q;

6:

BIPULL/RELABE(j θ ); 7: end while 8: END. f 0 (j θ , n λ ) := l 0 (j θ , n λ ) and adds node j θ to the rear of Q for all (j θ , n λ ) ∈ E 0 7:

d(n λ ) := 2w 2 + 1; 8: end for 9: END. G = (N 1 ∪ N 2 , A, l, u) in O(n 2 m + n 3
2 ) how is specified in table from Section 4. We apply the generalizate bipartite FIFO preflow algorithm in static reduced expanded bipartite network R 0 . Hence the algorithm has the complexity O(w 2 ε 0 + w 3 2 ). From Section 4 we have w 2 = n 2 T and ε 0 ≤ mT . As a result the algorithm has the complexity O(n 2 mT 2 + n 3 2 T 3 ). We specify that in the first phase the feasible flow f 0 is zero flow and in the second phase the feasible flow f 0 is the feasible flow f 0 determined in the first phase.

Example

The support digraph of bipartite dynamic network is presented in Figure 2 and time horizon being set T = 5, therefore H = {0, 1, 2, 3, 4, 5}. The transit times h(i, j; t) = h(i, j), t ∈ H, the lower bounds e(i, j; t) = e(i, j) and the upper bounds (capacities) q(i, j; t) = q(i, j), t ∈ H for all arcs are indicate in Table 2 We have N 1 = {2, 3, 4, 7} and N 2 = {1, 5, 6}. Applying the GBFIFOP algorithm in the first phase and the second phase we obtain the flows f 0 (i t , j θ ), f * 0 (i t , j θ ) (feasible flow, minimum flow) which are indicated in Figure 3. We have W 1 = {2 1 , 2 2 , 3 2 , 3 3 , 4 1 , 7 4 , 7 5 } and W 2 = {1 0 , 1 1 , 5 2 , 5 3 , 6 3 , 6 4 }. A maximum (1 0 , 1 1 ) -(7 4 , 7 5 ) cut in static network R 0 is Y 0 , Ȳ0 = (Y 0 , Ȳ0 ) ( Ȳ0 , Y 0 ) with Y 0 = {1 0 , 1 1 , 2 1 , 2 2 } and Ȳ0 = {3 2 , 3 3 , 4 1 , 5 2 , 5 3 , 6 3 , 6 4 , 7 4 , 7 5 }. Hence Y 0 , Ȳ0 = {(1 0 , 3 2 ), (1 0 , 4 1 ), (1 1 , 3 3 ), (2 1 , 5 2 ), (2 1 , 6 3 ), (2 2 , 5 3 ), (2 2 , 6 4 )}∪∅. We have w0 = f * 0 (Y 0 , Ȳ0 ) -f * 0 ( Ȳ0 , Y 0 ) = 14-0 = 14 and l 0 (Y 0 , Ȳ0 )-u 0 ( Ȳ0 , Y 0 ) = 14 -0 = 14. Hence f * 0 is a minimum flow.

1 :

 1 PROCEDURE PREPROCESS; 2: BEGIN 3: f 0 is a feasible flow in R 0 ; Q := ∅; 4: compute the exact distance function d in residual network Ĝ; 5: for λ ∈ H n do 6:

Figure 2 :

 2 Figure 2: The suport digraph of network D = (N, A, h, q)

Figure 3 :

 3 Figure 3: The network R 0 = (V 0 , E 0 , u 0 ).

  .

	(i, j)	(1, 2)	(1, 3)	(1, 4)	(2, 5)	(2, 6)	(3, 6)	(4, 5)	(4, 6)	(5, 3)	(5, 7)	(6, 7)
	h(i, j)	1	2	1	1	2	1	2	3	1	1	1
	e(i, j)	5	1	2	3	2	1	1	1	2	0	3
	q(i, j)	8	3	5	6	3	5	3	3	3	7	9

Table 2 :

 2 The functions h, e, q

6:

if (i t , j θ ) is admissible arc then 7:

select the first arc (k τ , i t ) in E - 0 (i t ) with r0 (k τ , i t ) > 0;

8:

if (k τ , i t ) is admissible arc then 9:

pull α := min {-e(j θ ), r0 (i t , j θ ), r0 (k τ , i t )} 10:

units of flow over the arcs (i t , j θ ), (k τ , i t );

adds node k τ to the rear of Q;

select the next arc in E - 0 (j θ )

24:

else 25:

d(j θ ) := min d(i t ) + 1|(i t , j θ ) ∈ E - 0 (j θ ), r0 (i t , j θ ) > 0 ; end if 30: until e(j θ ) = 0 or β = 0 31: if e(j θ ) > 0 then 32: adds node j θ to the rear of Q; 33: end if 34: END.

Figure 1: The generalizate bipartite FIFO prefow (GBFIFOP) algorithm

We remark that any path in the residual network R0 = (V 0 , Ẽ0 , r 0 ) can have at most 2w 2 + 1 arcs. Therefore we set d(n λ ) := 2w 2 + 1 in PROCEDURE PREPROCES.

The correctness of the GBFIFOP algorithm results from correctness of the algorithm for minimum flow in bipartite network [START_REF] Ciurea | Improved algorithms for minimum flows in bipartite networks[END_REF].

Theorem 2. The GBFIFOP algorithm which determine a minimum flow into bipartite dynamic networks D = (N, A, h, q) has the complexity O(n 2 mT 2 + n 3 2 T 3 ). Proof. In Section 3 we specify that the minimum flow problem for T time periods in dynamic network D = (N, A, h, e, q) is equivalent with the minimum flow problem in static reduced expanded network R 0 = (V 0 , E 0 , l 0 , u 0 ). The networks D and R 0 are bipartite with


Anthony Giacoma

David Dureisseix
email: david.dureisseix@insa-lyon.fr.

Anthony Gravouil

Michel Rochette

A multiscale LATIN/FAS algorithm with time-space model reduction for frictional contact problems

Keywords: Frictional contact, LATIN, Model reduction, Proper Orthogonal Decomposition (POD), Multigrid (MG-FAS)

come     A multiscale LATIN/FAS algorithm with time-space model reduction for frictional contact problems

Introduction

Simulations of mechanical systems involving frictional contact between flexible bodies can lead to strong difficulties. Non-linear behavior for bodies, large contact zones, and somehow large friction coefficient lead to numerical issues. Frictional contact problems are characterized by non-smooth behaviors at contacting boundaries of solids (discontinuities at contact interfaces for both normal and tangential components). These strong non-linearities dare the computational mechanic community to solve efficiently such problems. Over the last decades, various resolution methods have been developed. Contact problems can be formulated as an optimization problem subject to constraints with the augmented Lagrangian formulation. Such formulation can be solved using augmented Lagrangian methods [START_REF] Fortin | Augmented Lagrangian methods: applications to the numerical solution of boundary-value problems[END_REF][START_REF] Simo | An augmented lagrangian treatment of contact problems involving friction[END_REF][START_REF] Laursen | Formulation and treatment of frictional contact problems using finite elements[END_REF][START_REF] Wriggers | Finite element algorithms for contact problems[END_REF][START_REF] Heinstein | An algorithm for the matrix-free solution of quasistatic frictional contact problems[END_REF][START_REF] Alart | A scalable multiscale latin method adapted to nonsmooth discrete media[END_REF][START_REF] Bussetta | The adapted augmented lagrangian method: a new method for the resolution of the mechanical frictional contact problem[END_REF][START_REF] Pietrzak | Large deformation frictional contact mechanics: continuum formulation and augmented lagrangian treatment[END_REF]. Non-linear Gauss-Seidel solvers are also developed and applied to problems involving multi-contact bodies (rigid or flexible) [START_REF] Jourdan | A gauss-seidel like algorithm to solve frictional contact problems[END_REF][START_REF] Jean | The non-smooth contact dynamics method[END_REF][START_REF] Dostál | Implementation of the fixed point method in contact problems with coulomb friction based on a dual splitting type technique[END_REF]. But these last methods are not the most efficient to solve large systems (poor rate of convergence). Gradient methods are developed and adapted [START_REF] Barboteu | A domain decomposition strategy for nonclassical frictional multi-contact problems[END_REF][START_REF] Dureisseix | A numerically scalable domain decomposition method for the solution of frictionless contact problems[END_REF][START_REF] Renouf | Conjugate gradient type algorithms for frictional multicontact problems: applications to granular materials[END_REF][START_REF] Dostál | Minimizing quadratic functions subject to bound constraints with the rate of convergence and finite termination[END_REF] (with projection of the descent direction on feasible region for contact laws or active-set violation control strategy ...) in order to reproduce well-known performances of such methods for linear systems. But the main drawback of them is the difficulty to provide a convergence proof.

In spite of the non-smooth property of contact problems, widely used Newton's methods (generalized Newton's methods) are studied [START_REF] Curnier | A generalized newton method for contact problems with friction[END_REF][START_REF] Alart | Méthode de newton généralisée en mécanique du contact[END_REF][START_REF] Li | Smoothing newton and quasi-newton methods for mixed complementarity problems[END_REF][START_REF]Generalized newtons methods for the approximation and resolution of frictional contact problems in elasticity[END_REF]. Such methods proved their efficiency and robustness in large scale elasticity problems with contact and friction. Once again and due non-smooth property of contact problems, convergence proof of generalized Newton's methods are hard to provide. Another method consists in casting the contact problems into a linear complementarity problem (LCP) and using specific solvers such as active-set methods, Lemke's algorithm, projected successive over relaxation (PSOR). In [START_REF] Anitescu | Formulating dynamic multi-rigid-body contact problems with friction as solvable linear complementarity problems[END_REF], an LCP formulation is used to solve a frictional contact problem by faceting the Coulomb's cone. Nevertheless, this approximation is tough and not efficient (the problem to solve becomes larger). Nowadays, LCP formulations are suited to frictionless problems.

Generally speaking and even within a quasi-static context, all these nonlinear solvers can lead to prohibitive time of computations. Acceleration strategies based on multigrid methods were proposed [START_REF] Alart | Solution of frictional contact problems using ilu and coarse/fine preconditioners[END_REF][START_REF] Lebon | Multigrid methods for unilateral contact problems with friction[END_REF]. Computational costs spur recent and intensive works on efficient model reduction techniques in various fields [START_REF] Kerfriden | Bridging proper orthogonal decomposition methods and augmented newton-krylov algorithms: An adaptive model order reduction for highly nonlinear mechanical problems[END_REF][START_REF] Kunisch | Galerkin proper orthogonal decomposition methods for parabolic problems[END_REF][START_REF] Ryckelynck | On the a priori model reduction: Overview and recent developments[END_REF][START_REF] Amsallem | Nonlinear model order reduction based on local reduced-order bases[END_REF]. But and due to the non-differentiable nature of frictional contact, application of such methods seems to be heretic.

The aim of this work is to propose a strategy accelerating solution for frictional contact problems in the finite element framework embedding model reduction techniques also well-suited for parametric studies. Linear elastic, homogeneous, isotropic behavior and quasi-static evolution are assumed for the bodies. Frictional contact laws involve a non-linear and non-smooth behavior at the boundary of the body. To solve this mechanical problem, the nonincremental LATIN (LArge Time INcrement) method is used. This method is well-known for its ability to solve difficult non-linear large problems (nonlinear material, contact problems...) [START_REF] Ladevèze | Nonlinear Computational Structural Methods: New Approaches and Non-Incremental Methods of Calculation[END_REF][START_REF] Champaney | Large scale applications on parallel computers of a mixed domain decomposition method[END_REF] with a global time-space approach. This method is closed to augmented Lagrangian methods. Its great advantages are non re-factorization of matrices (stiffness matrix remains constant through LATIN iterations) and explicit subsequent iterations (no iterations to handle the non-linear behavior solved only on the contacting boundary). So the cost of a LATIN iteration is low even if the number of iterations can be elevated as augmented Lagrangian methods. We will try to exploit scale separability of contact problems as observed in [START_REF] Galland | A global model reduction approach for 3D fatigue crack growth with confined plasticity[END_REF][START_REF] Galland | An adaptive model reduction approach for 3D fatigue crack growth in small scale yielding conditions[END_REF] for fracture mechanics. The set of vector solutions can be compressed into a small number of vectors. Moreover, each of them describes a particular scale of the problem (e.g. scale of the global structure vs. scale of the crack). This phenomenon fosters multiscale strategies. At first, we propose to illustrate multiscale contents of solutions for three different contact problems. Then, some recalls about the LATIN method will be given in order to introduce the specific approach we propose. Subsequently, applications with convergence results on those three cases will be performed in order to highlight the feasibility of our strategy and its aptitude in the parametric studies framework.

Ω ∂ 2 Ω ∂ 1 Ω ∂ 3 Ω n g e 2 e 1

Reference problem

A flexible body is contacting a rigid basement (figure 1). We define Ω, the domain occupied by this solid with an open boundary (denoted by ∂Ω) at time t ∈ [0, T ] ⊂ R. This boundary can be split into three parts:

• ∂ 1 Ω: the part of the boundary where displacements u p are prescribed;

• ∂ 2 Ω = ∂Ω \ ∂ 1 Ω: the part of the boundary where external loads f ext are prescribed;

• ∂ 3 Ω: the part of the boundary where contact conditions occurs (contacting interface).

Strong formulation

Linear hypothesis are assumed (small displacements, homogeneous-isotropicelastic material). External loads are time dependent but quasi-static regime is assumed. The problem can be set as follows: find the displacement field and the Cauchy's stress field denoted respectively by u(x, t) and σ(x, t) satisfying:

• kinematic admissibility:

u(x, t) ∈ U [0,T ] with U [0,T ] = u has a finite energy, ∀(x, t) ∈ ∂ 1 Ω × [0, T ], u(x, t) = u p (1)
and the strain-displacement (in small perturbations) relationship reads:

∀(x, t) ∈ Ω × [0, T ] : ε(u) = 1 2 ∇u + ∇ T u (2)
• static admissibility: σ(x, t) is balanced with external force, i.e.

∀(x, t) ∈ ∂ 2 Ω × [0, T ] : σ n = f ext ∀(x, t) ∈ Ω × [0, T ] : div σ = 0 (3)
• constitutive laws:

-Hooke law for elasticity :

∀(x, t) ∈ Ω × [0, T ] : σ = Kε(u) ( 4
)
where K is the linear elasticity operator.

-Coulomb frictional law for contact. With v and λ, the traces of displacement field and contact force field (both normal and tangential) on the contacting interface ∂ 3 Ω, these conditions can be written formally as:

R(v, λ) = 0 (5)
To be more precise for the frictional behavior, at the contacting interface ∂ 3 Ω, we denote the local outward normal vector n. Then, the two fields can be split into a normal and a tangential part as:

∀(x, t) ∈ ∂ 3 Ω × [0, T ] : v(x, t) = u N n + u T λ(x, t) = λ N n + λ T (6)
R governs the contact behavior at the interface ∀(x, t) ∈ ∂ 3 Ω × [0, T ], according to Signorini's conditions and Coulomb's law on the displacement field and contact force field:

• Normal contact or Signorini's conditions given in [START_REF] Signorini | Questioni di elasticità non linearizzata e semilinearizzata. Rendiconti di Matematica e delle sue applicazioni[END_REF],

u N = u • n g Non-penetration condition λ N = (σ n) • n 0 Compressive contact force λ N (g -u N ) = 0 Complementarity condition ( 7
)
where g is the initial gap.

• Tangential contact or Coulomb's friction law [START_REF] Coulomb | Théorie des machines simples[END_REF]:

λ T 2 = -µλ N sliding: uT 2 = 0 and uT = -ρλ T with ρ > 0 λ T 2 -µλ N sticking: uT 2 = 0 ( 8
)
where µ ∈ [0 , +∞[ is the friction coefficient. Assuming a radial evolution of the displacement field per time step, the velocity can be replaced with the increment of displacement in the previous law.

Weak formulation

We define the following set of homogeneous kinematically admissible functions as:

u(x, t) ∈ U [0,T ] 0 with U [0,T ] 0 = {u has a finite energy and ∀(x, t) ∈ ∂ 1 Ω × [0, T ], u(x, t) = 0} (9)
Equilibrium ( 3) is equivalent to the following integral formulation with compatibility of contact force field and displacement field to contact conditions [START_REF] Bussetta | The adapted augmented lagrangian method: a new method for the resolution of the mechanical frictional contact problem[END_REF] and ( 8):

∀u * ∈ U [0,T ] 0 : Ω σ : ε(u * )dV - ∂2Ω f ext • u * dS - ∂3Ω λ • u * dS = 0 (10)

Discretized weak formulation

Using Hooke law for σ = Kε and the finite element approximation for displacement field yield to the so-called equilibrium equation with respect to frictional contact conditions at each node:

∀u * : T 0 u * T Ku dt = T 0 u * T f ext + f ctc dt with      v = Bu f ctc = B T λ R(v, λ) = 0 (
11) K is the stiffness matrix, f ext and λ are generalized forces corresponding to f ext and λ, and B is a boolean matrix mapping the global vector of nodal values to the values on contact boundary nodes. All in all, we have to find displacement field u and contact force field f ctc verifying [START_REF] Dostál | Implementation of the fixed point method in contact problems with coulomb friction based on a dual splitting type technique[END_REF] for all weight functions u * . Time interval is discretized into a subsequent of time steps t 0 k m describing a regular time stepping, i.e. t k+1 = t k + ∆t.

Reduced weak formulation and a posteriori analysis

Assuming u(x, t) being the solution of the discretized reference problem depicted in section 2 defined over the whole space (including boundaries) and time domains, the corresponding stress field is denoted by σ(x, t) and the contact force field by λ(x, t).

Discretized nodal displacements (with n degrees of freedom) are known for all time steps t 0 k m and are cast into a snapshot matrix defined as follows:

U =      u(x 1 , t 0 ) u(x 1 , t 1 ) • • • u(x 1 , t m ) u(x 2 , t 0 ) u(x 2 , t 1 ) • • • u(x 2 , t m ) . . . . . . . . . . . . u(x n , t 0 ) u(x n , t 1 ) • • • u(x n , t m )      = w(t 0 ) w(t 1 ) • • • w(t m ) (12)
U ∈ R n×m is a real rectangular matrix. According to the singular value theorem [START_REF] Eckart | The approximation of one matrix by another of lower rank[END_REF], U can be factorized using singular value decomposition (SVD):

U = ΥΣΦ T = Υ 1 • • • Υ n         σ 1 0 • • • 0 0 σ 2 . . . . . . . . . 0 0 • • • 0 σ r 0 • • • 0 0             Φ T 1 . . . Φ T m     (13)
with r = min(n, m), Υ ∈ R n×n is a unitary matrix containing left-singular vectors, Φ ∈ R m×m a unitary matrix containing right-singular vectors and Σ ∈ R n×m contains positive singular values σ i in increasing order. Moreover, left-singular vectors Υ i are space vectors which are combined with time functions Φ i scaled with the singular value σ i . Consequently, displacement field can be written as a sum of radial solutions:

U = Υ 1 σ 1 Φ T 1 + Υ 2 σ 2 Φ T 2 + • • • + Υ i σ i Φ T i i th radial solution + • • • + Υ r σ r Φ T r ( 14
)
This decomposition is unique up to an arbitrary sign for the couples of vectors (Υ i , Φ i ).

Given U and its SVD, and taking into account only K ≤ r first singular values, allows to define an approximation Ǔ such that:

Ǔ = Υ 1 σ 1 Φ T 1 + Υ 2 σ 2 Φ T 2 + • • • + Υ K σ K Φ T K (15)
According to the Eckart-Young's low rank approximation theorem [START_REF] Eckart | The approximation of one matrix by another of lower rank[END_REF], Ǔ is the best approximation of rank K of U according to the Frobenius' norm. It corresponds to the discrete POD of U. Moreover, if K = r then Ǔ = U.

Consequently, we are able to build a basis of relevant functions from Υ containing the K first left-singular space vectors. Moreover, the first vectors are the most energetic ones (i.e. they describe globally the solution and contain a great fraction of the elastic energy), as it will be exemplified in the following.

The elastic energy of the displacement field solution integrated over the studied time interval is given by

E = 1 2 T 0 w(t) T Kw(t) dt (16)
It can be distributed among radial solutions as E = r i=1 E i where

E i = σ i 2 2 Υ T i KΥ i T 0 Φ T i Φ i dt ( 17
)
corresponds to the elastic energy contained in the i th radial solution.

3 Application to 2D frictional contact problems suited to approximate it. The following examples are aimed to exemplify this on 2D frictional contact problems once their space-time solutions have been obtained. The first test model is a single contact case, the second one represents a larger structure when a part of it encounters contact conditions while the last one uses a more complex loading which is known to be less easily represented on a reduced radial base. Each time, we discuss the contribution of the different SVD modes and of their scale lengths. Linear 3-node triangular finite elements are used for the spatial discretization.

Material properties, simulation characteristics and parameters are given in table 1, for the upcoming three cases.

Hemicylindrical body contacting a rigid basement (case

A)

The first model (figure 2) consists in an elastic hemicylindrical body contacting a frictional rigid basement. We prescribe displacement on the upper boundary such that a large contact zone occurs. The scale of the structure and the scale of the contact zone are similar. Initially, the solid is not contacting. From time-step t 1 to t 40 , it approaches and becomes pressed on the rigid basement. Meanwhile, a tangential displacement is prescribed to activate a global frictional reaction at the contacting zone. We apply SVD to displacement field in order to generate singular vectors in space (figure 4) and in time (figure 3). We notice that for this case a few couple of modes (4 or 5) are needed to render 99.9% of the elastic energy of the highly reducible solution, which is therefore provided that we consider the space-time evolution in a while.

On figure 4, we notice that the three first modes refer to a large scale contribution acting on the whole body. Whereas after the fourth one, space modes are remarkably localized in the vicinity of contact area. After the 6 th mode, the characteristic length scale of the mode get closer and closer to the discretization length scale (i.e. the mesh size). Similarly, respective time functions have a global / local behavior on the time interval.

Ω ∂ 3 Ω R u p ∂ 1 Ω intial gap

Cantilevered indenter contacting a rigid basement (case B)

The second problem (figure 5) consists in an elastic clamped cantilevered hemicylindrical indenter contacting a frictional rigid basement. External loads are designed in order to get different status at the contacting nodes during the studied time interval (no contact, sliding, sticking). The size of the effective contact area is smaller than the scale of the structure.

In the same way, by applying the SVD to the displacement field, we get time and space functions. A few couple of modes (4 or 5) are also needed to render 99.9% of the elastic energy of the simulation (figure 6). The solution of this second problem is also highly reducible.

From figures 6 and 7, same conclusions as for the first simulation case can be drawn. Status switching seems to not affect reducibility of the solution. For this problem, the length scale of the global structure is more separated from the one of the contact area, and this is also depicted in the succession of spatial modes.

Cantilevered indenter contacting a rigid basement with a non radial external load (case C)

The last problem (figure 8) is the same as the second excepting for the external loads. We design a mobile half-sine lobe load moving along the upper boundary of the body with an increasing amplitude. It is well-known that this kind of solicitation is not efficiently described by a finite set of radial functions [START_REF] Dostál | Minimizing quadratic functions subject to bound constraints with the rate of convergence and finite termination[END_REF]. A soft contact occur at the bottom of the hemicylindrical part.

In figure 8, we introduce the space variable localizing he position of the

√ σ i . √ σ 2 Space mode Υ 3 √ σ 3 Space mode Υ 4 √ σ 4 Space mode Υ 5 √ σ 5 Space mode Υ 6 √ σ 6 0 0.55 1.1 0 0.55 1.1
R Ω ∂ 1 Ω : u p = 0 (clamped) ∂ 2 Ω f ext ∂ 3 Ω
intial gap R

Ω ∂ 1 Ω : u p = 0 (clamped) 0 L f ext = -a × lobe( )e 2 ∂ 3 Ω
intial gap lobe on the upper boundary. The evolution of according the time variable is prescribed as follows:

         t ∈ [0, T 4 ] : = 4L t T t ∈ [ T 4 , T 2 ] : = 4L(1 -t T ) t ∈ [ T 2 , 3T 4 ] : = 4L t T t ∈ [ 3T 4 , T ] : = 4L(1 -t T ) (18)
According to figure 9, the reduction of the solution is clearly less efficient than other problems. Twenty modes are needed (instead of five) to contain 99.9% of the elastic energy. Space modes on figure 10 are concerned with three scales of the problem: the body length scale, the contact zone length and the load length scale (the wavelength of the lobe) on the upper boundary where it is applied.

SVD functions and multiscale contents

Time-space model reduction emphasizes scale separability phenomenon. At least, two scales of the problem can be distinguished: the global scale of the structure and a local scale where contact occurs. First radial solutions obtained by SVD describe globally the displacement field (both in space but also in time) whereas higher order modes bring corrections enhancing precision at contact interface.

We can notice that time-space modes can be dominated by specific solicitations occurring during the deformation process. It is clear in the model A : the first and the second modes are linked to the compression on the rigid basement whereas the third one is linked to shearing effects involved by friction: the appropriate linear combination of these three modes for each time step covers roughly 95% of the elastic energy involved. Nevertheless and due to orthogonality properties and energetic considerations, SVD modes are not prescribed to be characterized by load components. Consequently, we have to distinguish SVD decomposition and linear superposition principle even though both span the exact solution.

Efficient resolution using reduced order model

Scale separability observation encourages the deployment of multiscale strategies [START_REF] Brandt | Multiscale solvers and systematic upscaling in computational physics[END_REF]. Generally speaking, they consist in distinguishing a global behavior occurring on the whole structure and a behavior localized on specific area. For instance, in [START_REF] Pierres | A two-scale extended finite element method for modelling 3D crack growth with interfacial contact[END_REF], a cracked body is studied taking into account friction contacts between crack faces. Behavior at interfaces (crack faces) is treated as an autonomous problem (own discretization, constitutive law, internal variables) and then linked in a weak sense to the global problem. In [START_REF] Kerfriden | Local/global model order reduction strategy for the simulation of quasi-brittle fracture[END_REF], a local/global strategy is used by solving the global problem in a reduced space (obtained by POD) while the local problem (where localized damage phenomenon occurs) is fully resolved. Then, both of the global and local problems are coupled. Other efficient approaches are multigrid solvers suitable for various non-linear problems [START_REF] Venner | Multilevel Methods in Lubrication[END_REF][START_REF] Cavin | Efficient FEM calculation with predefined precision through automatic grid refinement[END_REF][START_REF] Rannou | Prise en compte d'effets d'échelle en mécanique de la rupture tridimensionnelle par une approche X-FEM multigrille localisée nonlinéaire[END_REF][START_REF] Biotteau | Three dimensional automatic refinement method for transient small strain elastoplastic finite element computations[END_REF]. For frictional contact problems, multigrid methods were proposed in [START_REF] Lebon | Multigrid methods for unilateral contact problems with friction[END_REF]. Various iterative solvers (called "smoothers") have a remarkable smoothing effect on high-frequency component of the error (relatively to the discretization). In contrast, low-frequency of the error is slowly damped by these solvers. Consequently, coarser discretization are introduced (coarse grids) to bring corrections to fine grids to increase convergence rate.

We can draw an analogy between the previous observations and multigrid solvers. We propose to consider global scale modes as coarse grid components. In the same way as multigrid methods, we try to bring time-space corrections from a subspace spanned with these first space modes obtained by an a posteriori analysis. Corrections are quickly computed on a reduced order model (defined by a reduced basis built on left-singular Υ vectors of SVD). Then, these corrections are transferred on the full order model (full dimensional model) to increase convergence rate of the iterative LATIN solver [START_REF] Ladevèze | Nonlinear Computational Structural Methods: New Approaches and Non-Incremental Methods of Calculation[END_REF] (considered as the smoother).

In this study we rely on the SVD modes that are precomputed to exemplify these effects on the solvers. Of course, such modes, taking into account their scale properties could have been computed by other means (for instance, coarse grid computations, analytical local solutions ...) but this is out of the scope of this article.

LATIN method

One of the major points of our strategy is the non-incremental LATIN method introduced by Ladevèze and designed for time-dependent problems. Champaney et al. in [START_REF] Champaney | Une nouvelle approche modulaire pour l'analyse d'assemblages de structures tridimensionnelles[END_REF][START_REF] Champaney | Large scale applications on parallel computers of a mixed domain decomposition method[END_REF][START_REF] Champaney | Modular analysis of assemblages of three-dimensional structures with unilateral contact conditions[END_REF] used this method to solve problems involving contact between substructures. In [START_REF] Ladevèze | A modular approach to 3-D impact computation with frictional contact[END_REF][START_REF] Lemoussu | A 3D shock computational strategy for real assembly and shock attenuator[END_REF], LATIN method is used to solve 3D impacts with frictional contacts with plasticity interfaces. Herein, we give a short recall of the LATIN method directly formulated for frictional contact problems with a two-field weak formulation. Stabilized formulations (three-field weak formulation with stabilization operator) can be used as in [START_REF] Trollé | Optimization of a stabilized X-FEM formulation for frictional cracks[END_REF].

LATIN method is based on three principles illustrated on figure 11 (P3) Radial approximation. The unknowns are represented as a sum of products between a time function and a space functions. This principle is not use herein.

To specify the equations to solve at global and local stage, we use a displacement formulation. Velocity formulation well-suited for dynamical framework could also be used [START_REF] Ladevèze | A modular approach to 3-D impact computation with frictional contact[END_REF].

Global stage. This stage consists in finding a displacement field v(x, t) and a contact force field λ(x, t) defined over ∂ 3 Ω × [0, T ] verifying the balance of the solid with the search direction [START_REF]Generalized newtons methods for the approximation and resolution of frictional contact problems in elasticity[END_REF].

(E -) : λ -λ = -k(v -v) ⇔ λ T λ N -λT λN = - k T I d 0 0 k N I d v T v N -vT vN ( 19
)
I d is the identity matrix of size equal to the number of contacting nodes. k N and k T are parameters of the method defined at each node. It is similar to the penalty parameter in an augmented Lagrangian formulation. Different values for the normal and tangential problem can be chosen. These values influence only on the convergence rate and optimal values are related to mechanical properties of the studied body [START_REF] Ladevèze | Nonlinear Computational Structural Methods: New Approaches and Non-Incremental Methods of Calculation[END_REF][START_REF] Champaney | Modular analysis of assemblages of three-dimensional structures with unilateral contact conditions[END_REF][START_REF] Boucard | A suitable computational strategy for the parametric analysis of problems with multiple contact[END_REF][START_REF] Trollé | Optimization of a stabilized X-FEM formulation for frictional cracks[END_REF]. Generally speaking and according to [START_REF] Boucard | A suitable computational strategy for the parametric analysis of problems with multiple contact[END_REF], optimum values are searched as follows:

Γ A s s i = (v i , λ i ) ŝi = (v i , λi ) E + s i+1 E - ŝi+1 E +
k = EL c ( 20
)
with E the Young's modulus and L c a contact characteristic dimension. k has the dimension of a stiffness. A parametric study was carried out for case A (figure 12) to identify the optimal value of k for both LATIN and FAS/LATIN algorithm (algorithm proposed). Displacement field and contact force field verifying equilibrium are given by [START_REF] Alart | Solution of frictional contact problems using ilu and coarse/fine preconditioners[END_REF].

     Ku = f ext + B T λ v = Bu λ = λ -k(v -v) (21)
Local stage. The local stage is an updating stage for the contacting interface. Its solution is explicit and is deduced from search direction [START_REF] Lebon | Multigrid methods for unilateral contact problems with friction[END_REF] and non-linear contact conditions. We can remark that FAS/LATIN doesn't modify the optimal value and spreads the "optimal zone" Normal components

λN = c N - vN = j- 1 k N c N + Tangential components Sliding: c T 2 > µ| λN | Sticking: c T 2 µ| λN | t = c T / c T 2 λT = µ| λN |t vk T = vk-1 T - 1 k T c T 2 -µ| λN | t λT = c T vk T = vk-1 T Table 2: Solutions of the local stage. (E + ) : λ -λ = k(v -v) ⇔ λT λN - λ T λ N = k T I d 0 0 k N I d vT vN - v T v N (22)
Equivalently, search directions can be rewritten for both tangential and normal contact. For each contacting node, the search direction yields to equations [START_REF] Kerfriden | Bridging proper orthogonal decomposition methods and augmented newton-krylov algorithms: An adaptive model order reduction for highly nonlinear mechanical problems[END_REF] and [START_REF] Kunisch | Galerkin proper orthogonal decomposition methods for parabolic problems[END_REF]. v k T refers to the nodal displacement on the contacting interface in the tangential plane at time step t k .

c N = λ N -k N (v N -j) = λN -k N (v N -j) (23)
c T = λ T -k T (v k T -vk-1 T ) = λT -k T (v k T -vk-1 T ) (24)
Solutions of the local stage are given in table 2 and have to be computed for each node belonging to the contacting interface and for each time step.

Initialization and stopping criterion. To initiate the iterative process, a linear contactless prediction is computed for s 0 . We use stopping criterion [START_REF] Ryckelynck | On the a priori model reduction: Overview and recent developments[END_REF] introduced in [START_REF] Ribeaucourt | A new fatigue frictional contact crack propagation model with the coupled X-FEM/LATIN method[END_REF]. Iterative process is stopped once I < ε.

I = max s N -ŝN 2 ∞ s N 2 ∞ + ŝN 2 ∞ , s T -ŝT 2 ∞ s T 2 ∞ + ŝT 2 ∞ with s N,T 2
∞ = max x,t 1 k λ 2 N,T + ku 2 N,T (25)
where k is the search direction parameter. Note that subscripts N and T refer respectively to normal and tangential components at a given node belonging to the contacting interface. Criterion ( 25) is a very tough convergence indicator

I G = T 0 w(t) -w ref (t) T K w(t) -w ref (t) dt T 0 w ref (t) T Kw ref (t)dt ( 26
)
where w ref is the reference solution for displacement field. As it proceeds in a global space-time approach implying a specific error appreciation (e.g. convergence criterion, constitutive law error...), LATIN method cannot be directly compared to usual methods. The LATIN method is summarized up in the algorithm 1 and convergence indicators are plotted for all described cases on figure 13.

Algorithm proposal

The full approximation scheme (FAS) multigrid method [START_REF] Venner | Multilevel Methods in Lubrication[END_REF] is a well-suited solver for non-linear problems. Difficulty of a such a solver is to define a coarse problem where low-cost corrections (computationally speaking) are calculated. For that purpose, we propose the following algorithm:

8 Solve K + B T kB u = f ext + B T ( λ + kv) 9 Extract displacement interface v = Bu 10 Search direction λ = λ -k(v -v) 11 end 12
Compute convergence criterion I(v, λ, v, λ)

13 end 14 if i = i max then No convergence
We have to define operators to transfer quantities between the full order space and the reduced order one. These are analog to intergrid operators (restriction and interpolation operators). Thus, we define:

• restriction operator projecting dual quantities from FOM to ROM: R (e.g. R = Υ T )

• restriction operator projecting primal quantities from FOM to ROM: R (e.g. R = Υ T )

• prolongation operator projection primal quantities from ROM to FOM: P = Υ (e.g. P = Υ)

As mentioned in [START_REF] Rannou | Prise en compte d'effets d'échelle en mécanique de la rupture tridimensionnelle par une approche X-FEM multigrille localisée nonlinéaire[END_REF], prolongation and restriction operators have to fulfill energy conservation condition between ROM and FOM. Consequently, we have to verify :

R = P T ( 27
)
Fine problem (FOM). LATIN iterations are considered as relaxations on the FOM (or the finest description of the problem). The solution of the FOM is the solution of the reference problem. Consequently, convergence indicator and errors have to be evaluated at this stage. We chose to end relaxations by a local stage that produces a solution ( u, λ) satisfying the contact behavior but not admissibility conditions. There is therefore a balance residual defined as:

res = f ext + B T λ -K u ( 28
)
Coarse problem (ROM). We have to define it in such a way that we get relevant corrections for both displacement field and contact force field ( u, λ) obtained from previous relaxation stage. For this purpose, we propose a mixed version of the LATIN solver and FAS coarse solver: a two-stage iterative process on corrective increments to satisfy contact conditions and reducing equilibrium residual.

• First stage. If one gets a first approximation of the correction ( ∆u, ∆λ), we search for a better approximation (∆u, ∆λ) trying to reduce the balance residual res which is:

res = f ext + B T λ + ∆λ -K u + ∆u (29)
using an additional search direction:

∆λ = ∆λ -k(∆v -∆v) (30)
one gets

res = res + B T ∆λ + k ∆v -K + B T kB ∆u (31)
The FAS scheme nullifies the coarse component of this residual (i.e. Rres = 0) with the representation:

∆u = P Ψ -R u (32)
Therefore it leads to solve:

R K + B T kB P (Ψ -R u) = R res + B T ∆λ + k ∆v ⇔ R K + B T kB P Ψ = R res + B T ∆λ + k ∆v + K + B T kB PR u (33)
Once Ψ is obtained, the correction (∆u, ∆λ) is easily computed and the corresponding full fields are ( u + ∆u, λ + ∆λ).

• Second stage. It consists in obtaining a new correction ( ∆u, ∆λ) to satisfy contact conditions for the full fields. Here, we simply reuse the local stage to get full fields ( u, λ) and corrective increments ( ∆u, ∆λ).

The pseudocode is described in algorithm 2.

Application

Knowing a posteriori a reduced basis defining the coarse problem, we compare the proposed strategy and the former non-linear LATIN solver in terms of number of relaxations for each three problems presented before. Sets of parameters used for the application of the proposed algorithm are described on table 3. The value of search direction parameters are set to optimal values (the less number of LATIN iterations are needed to solve the problem). 1) , λ (j-1) ) for ν 1 relaxations do

7 if j = 1 then (v, λ) ← (Bu (j-
8
for each time step do Solve local stage to get (v, λ) satisfying R(v, λ) = 0 for each time step do

9 Solve K + B T kB u = f ext + B T ( λ + kv) 10 Extract displacement interface v = Bu 11 Search direction λ = λ -k(v -v) 12 end 13 end 14 Solve local stage to get (v, λ) satisfying R(v, λ) = 0 15 Unbalanced solution assembly ( u, λ) ← u Ω\∂3Ω u ∂3Ω = v , λ 16
Create or update operators: R,R,P As no analytic solution exists, the reference is a solution evaluated with a precision of ε = 10 -12 . The number of SVD vectors defining the reduced order model is studied.

Parameters Case A Case B Case C µ 0.15 ν 1 1 1 3 ν 2 1 1 3 i C,max 15 15 15 ε C -1 -1 -1

Discussion

Results on figure 14 were obtained by enforcing coarse step iterating 15 times for each cycle (one cycle = ν 1 relaxations + coarse step + ν 2 relaxations). The threshold ε C is therefore set to a negative value to prescribe a given number of iterations for the coarse correction. In this study, we chose a rather large number of iterations for the coarse step to get an "exact" coarse correction. Indeed, it would be possible to make it only run until a threshold fixed by a level of accuracy (convergence criterion ε C,max > 0). Generally speaking, the proposed algorithm converges faster than the LATIN method (number of relaxations can be divided from 2 to 10 depending on parameters and desired precision).

On figure 14, the first comment concerns the content of the reduced basis and the definition of operators P, R, R. A reduced order model built exclusively on a posteriori SVD vectors leads to a weak acceleration of convergence even if the size of the basis is large (for case A, the 48 SVD vectors correspond to all vectors associated to non-null singular values). As a consequence, other vectors have to be introduced to define the coarse problem. Indeed, at each cycle, an approximated solution is built which could be not properly represented in the ROM (SVD vectors are built from the solution of the problem and not from the iterates). As a consequence, relevant coarse corrections can not be brought to the FOM from ROM. For that purpose, we have to add other vectors spanning the orthogonal residual (to SVD basis) of iterated solution. Arbitrary, we complete SVD vectors with 15 other vectors spanning the residual. At each cycle, these vectors are updated and the size of the coarse problem remains the same. On figure 14, we study the influence of the number of SVD vectors fixing the number of "residual vectors" to 15.

Generally speaking, the more SVD vectors there are, the larger is the number of used SVD vectors, the smaeller is the number of iterations to converge, but the larger is the coarse problem. Consequently, we have to make a compromise between the size and the relevancy of the coarse model. A too few number of vectors defining the coarse problem makes the coarse step inefficient (no accelerating effect) or, even worse, tends to struggle convergence of the algorithm (coarse correction becomes irrelevant). The additional residual vectors allow to counterbalance a possible lake of relevancy of SVD vectors and ensure convergence.

FAS/LATIN convergence plots exhibit different regimes in the converging process: an acceleration of the convergence is clearly seen during a first group of iterations, then the convergence rate of the former classical LATIN method is recovered. The larger the reduced basis is, the longer the accelerating regime pertains. This can be explained by the multiscale content of the reduced basis as it was exemplified in section 3. During the first iterations, global scale of the problem is captured (and contained in the first modes of the SVD) involving a steep slope of the error. Then, the slope reduces depending on the content of the coarse problem. During this second regime, local scales are solved and higher order SVD vectors are used to compute coarse correction. The last regime shows a convergence rate roughly similar to the former LATIN method convergence rate. Only relaxation steps are useful and the coarse step does not bring corrections anymore.

Besides the definition of the coarse problem, other parameters can influence the behavior of the proposed algorithm, namely: the number of relaxations before and after the coarse step. This parameter allows to give a prominent role to the FOM or the ROM. Generally speaking, more relaxations steps renders the algorithm slower but more robust as shown on figure 15: a too few relaxations makes fail the convergence of the problem (up to a certain level of accuracy: for instance on figure 15, ν 1 = 1 and ν 2 = 0 stops to converge at 10 -4 ). This may also depend on the number of SVD vectors that are selected. According to results on figure 15, the less (but a sufficient) number of relaxations is the optimal configuration. Note that the more relaxations there are for ν 1 and ν 2 , the closest to the LATIN method our strategy becomes (the less the coarse step occurs).

Parametric study

Product design requires several optimization phases. In particular, physical parameters has to be choose among discrete or continuous interval values in order to minimize a cost function. To achieve this goal, response surface methodology (using e.g. basic interpolations or kriging) can be used to identify the optimum set of parameters. A such methodology needs several evaluations of the cost function (a sample of set of parameters) which is a very expensive phase. Thus, acceleration strategies might be a very useful tool to save up computational cost. In this respect, we propose to deploy our strategy. Originally, the LATIN method is very well-suited for such parametric studies. The main idea consists in re-using a previous solution for a given set of parameters to start LATIN iterations for a new set of parameters. This way, a large amount of relaxations (LATIN iterations on the full order model) can be spared as it is described in [START_REF] Boucard | A suitable computational strategy for the parametric analysis of problems with multiple contact[END_REF][START_REF] Odièvre | A parallel, multiscale domain decomposition method for the transient dynamic analysis of assemblies with friction[END_REF][START_REF] Soulier | A multiparametric strategy for the two step optimization of structural assemblies[END_REF].

In the other hand, reduced basis methods are also used to accelerate solutions [START_REF] Bolzon | An effective computational tool for parametric studies and identification problems in materials mechanics[END_REF][START_REF] Deparis | Reduced basis method for multi-parameter-dependent steady navier-stokes equations: Applications to natural convection in a cavity[END_REF][START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF]. Generally speaking, these consist in solving once the full-order problem associated to a first set of parameters and extract a reduced basis from computed solution (an expensive first phase). Then, the reduced basis is used (and eventually enriched) to accelerate solutions of other sets of parameters.

We wish study different values for the friction coefficient µ ∈ {µ m = 0, 0.01, 0.5, 0.7, µ M = 1} focusing on the case A, knowing a first result for µ = 0.15.

First, we propose a reduced basis approach using our strategy. Then, a combination between reusing previous solution to compute new ones (as proposed in [START_REF] Boucard | A suitable computational strategy for the parametric analysis of problems with multiple contact[END_REF][START_REF] Odièvre | A parallel, multiscale domain decomposition method for the transient dynamic analysis of assemblies with friction[END_REF]) and our strategy is tested. All these approaches are compared to a brute force method.

Brute force method

For each µ, a calculation is carried out with the non-linear full-order LATIN method (figure 16).

Reduced basis methods

Two approaches can be used [START_REF] Heyberger | Multiparametric analysis within the proper generalized decomposition framework[END_REF]:

(I) First ROM-POD approach. A reduced basis is built from the solution with µ = 0.15. Then, this basis is used to solve all other problems associated to the other friction coefficients.

(II) Second ROM-POD approach. A reduced basis is constructed from resolutions of problem associated to extreme values for the varying parameter (µ = 0 and µ = 1). Then, we use this basis to solve all problems associated to all other µ. SVD reduced basis is built by concatenating and orthogonalizing the first 20 SVD vectors of the solution for µ = 0 and the first 20 SVD vectors of the solution for µ = 1.

For both ROM-POD approaches, we update a set of vectors of the reduced basis on the fly without a posteriori knowledge, but we keep the same size of the reduced basis.

Results obtained using the proposed algorithm are depicted on figure 16. Clearly, using our strategy to solve similar problems with varying parameters allows to reduce the number of relaxations in comparison to brute force method. Bringing coarse corrections on reduced basis (global scale of the structure) is relevant for all problems differing only on a localized scale (contacting boundary). Generally speaking, in both approaches, simulations with small coefficients of friction (µ = 0 and µ = 0.01) are highly accelerated. But for higher ones, accelerating effect is weaker for both approaches. That can be explained by the toughness of tangential problem which is difficult for large coefficient of friction requiring particular local modes. Moreover the SVD basis is built with a problem where the friction coefficient is relatively small (µ = 0.15). Approach (II) is globally more efficient than approach (I). Indeed, coarse model of approach (II) is built from a larger panel of solutions and is spanned with more vectors.

Enriching the collection of SVD vectors with vectors spanning the unprojected residual adapts the definition the coarse stage to a specific problem to solve. It can be seen as an a priori method which tries to maintain the relevancy of the subspace spanned by the the coarse model.

Reusing previous solutions

In this paper, we initialize the LATIN and the FAS/LATIN method with the linear elastic solution (kinematically and statically admissible) without contact conditions. Now, as in [START_REF] Boucard | A suitable computational strategy for the parametric analysis of problems with multiple contact[END_REF][START_REF] Odièvre | A parallel, multiscale domain decomposition method for the transient dynamic analysis of assemblies with friction[END_REF], we initialize the LATIN method or FAS/LATIN with the previous solution computed (figure 17). We compare the number of relaxations done to reach a precision of 10 -4 for all approaches (figure 18). For the FAS/LATIN, the first reduced basis approach is used (basis computed for µ = 0.15).

We firstly note that the solution reusing doesn't reduce the number of relation for LATIN method. Indeed, the linear elastic solution (used originally) of the studied case is already a good and relevant first guess of the low-frictional solutions for case A. In case of µ = 0 (frictionless problem), the solution for µ = 0.15 is irrelevant (mid-frictional problem) and the number of relaxations increases. Nevertheless, for high-frictional problems (µ = 0.5 or 0.7), initializing with a (high-)frictional solution is this time a relevant choice. A quite obvious conclusion can be drawn, reusing solutions is suited for initializing LATIN to solve quite similar problems. In addition to this observation, the FAS/LATIN algorithm provides a very interesting gain and shows its ability to perform such parametric problems.

Conclusion

This article shows that a frictional contact problem can be reduced and efficiently described with a set of space-time radial functions (radial solutions) provided that a time-space model reduction is used. Analyses of each radial solution reveal a clear scale separability of the solution. This observation justifies the use of multiscale approach. Then, an algorithm based on the non-linear FAS multigrid solver was proposed. The relaxation step is based on the non-linear LATIN solver and the non-linear coarse step is based on a reduced basis strategy. However, a sole POD basis is insufficient to increase efficiently the convergence rate even if the solution of one problem is highly reducible. This is why a few more vectors have to complete the POD basis: the iterates of solutions built by the non-linear iterative solver have to be taken into account. Consequently, a possible enrichment of the reduced basis is a key point to increase efficiently the convergence rate of the multiscale approach. Moreover, it has been shown that such an approach is also well-suited to parametric studies.

In the proposed approach, enrichment of POD basis vectors with residual vectors can be compared with an a priori method. This enrichment can be integrated by considering the so-called radial approximation (or separated representation). The efficiency of this approach has already been emphasized for thermal and mechanical problems [START_REF] Ladevèze | The LATIN multiscale computational method and the proper generalized decomposition[END_REF][START_REF] Heyberger | Multiparametric analysis within the proper generalized decomposition framework[END_REF][START_REF] Lemoussu | A 3D shock computational strategy for real assembly and shock attenuator[END_REF] which is close to PGD (Proper Generalized Decomposition) methods [START_REF] Chinesta | Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models[END_REF][START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF]. PGD methods operate without a priori knowledge of the solution (contrary to POD approaches). Many works underline the potential of such a priori approaches [START_REF] Heyberger | Multiparametric analysis within the proper generalized decomposition framework[END_REF][START_REF] Boucinha | Space-time proper generalized decompositions for the resolution of transient elastodynamic models[END_REF]. That point is currently being investigated.

Figure 1 :

 1 Figure 1: Problem set.

Figure 2 :

 2 Figure 2: (a) Model -(b) Prescribed displacement evolution -(c) Shape of the deformed body at t 50 .

Figure 3 :

 3 Figure 3: (a) The a posteriori analysis of the first 50 singular values and the cumulative elastic energy of the modes -(b) The a posteriori analysis of the first 6 right-singular time modes Φ i scaled by √ σ i .

Figure 4 :

 4 Figure 4: Shapes of deformed body for the left-singular space modes scaled by respective √ σ and the grayscale color map depicts the gradient of the amplitude of the displacement field.

Figure 5 :√Figure 6 : 6 0

 566 Figure 5: (a) Model -(b) External loads evolutions. Initially, the solid is not contacting. -(c) Shape of the deformed body at t 50 .

Figure 7 :

 7 Figure 7: Shapes of deformed body for the left-singular space modes scaled by respective √ σ and the grayscale color map depicts the gradient of the amplitude of the displacement field.

Figure 8 :

 8 Figure 8: Model C -(a) Model with a mobile external load. The load is a half-sine moving twice from = 0 to = L and back to = 0. Its amplitude varies linearly according to . At = 0 external load are null -(b) Shape of the deformed body when = L.

√Figure 9 : 6 0

 96 Figure 9: (a) The a posteriori analysis of the first 50 singular values and the cumulative elastic energy through modes -(b) The a posteriori analysis of the first 6 right-singular time modes Φ i scaled by √ σ i .

Figure 10 :

 10 Figure 10: Shapes of deformed body depicts the left-singular space modes scaled by respective √ σ and the grayscale color map depicts the gradient of the amplitude of the displacement field.

  : (P1) Separation of the linear and non-linear behavior. A denotes the set of solutions s = (v, λ) satisfying linear constitutive law, kinematic admissibility and static admissibility. These are defined on the whole time-space domains Ω × [0, T ]. Γ denotes the set of solutions ŝ = (v, λ) verifying frictional contact conditions and are defined locally at the contacting interface and on the whole time interval ∂ 3 Ω × [0, T ]. The solution of the problem is s ∈ A ∩ Γ. (P2) A two-staged iterative algorithm. The solution of the problem is searched with the construction of two sequences of approximations belonging alternatively to A and then Γ. At the i th iteration, global stage consists in finding s i+1 = (v i+1 , λ i+1 ) ∈ A with a search direction (s i+1ŝi ) = (v i+1vi , λ i+1 -λi ) ∈ E -. Note that ŝi = (v i , λi ) is known at this stage. Then, local stage consists in finding ŝi+1 = (v i+1 , λi+1 ) ∈ Γ with another search direction (ŝ i+1s i+1 ) = (v i+1v i+1 , λi+1λ i+1 ) ∈ E + . And s i+1 = (v i+1 , λi+1 ) is known from the previous stage.

Figure 11 :

 11 Figure 11: LATIN non-linear iterative solver.

Figure 12 :

 12 Figure12: Number of LATIN iterations to reach a precision of 10 -3 according different value for search directions (normal k N and tangential k T problem). We can remark that FAS/LATIN doesn't modify the optimal value and spreads the "optimal zone" 19

Figure 13 :

 13 Figure 13: Convergence rates according to the local convergence criterion I and the global one I G for cases A, B and C. providing a possible low convergence rate. Another far less restrictive criterion (a global energetic one) can be used:

Algorithm 1 : 6 i ← i + 1 7

 161 1. a few iterations (ν 1 relaxations) on the full order model (FOM) with the non-linear LATIN solver; 2. compute a coarse correction on a subspace spanned by space vectors (on a reduced order model or ROM); 3. a few iterations (ν 2 relaxations) on the full order model with the non-linear LATIN solver. LATIN solver for frictional contact Input: Numerical model Output: Admissible displacement field and contact force field 1 Initialization: i = 0 2 Initialization with a contactless linear elastic prediction. Solve Ku = f ext 3 Extract displacement interface v = Bu 4 Set contact force λ = 0 5 while I > ε and i < i max do for each time step do Solve local stage to get (v, λ) satisfying R(v, λ) = 0 for each time step do

Algorithm 2 : 6 j

 26 LATIN/FAS Input: Numerical model, reduced basis Υ Output: Admissible displacement field and contact force field 1 Initialization with a contactless linear elastic prediction. Solve Ku = f ext 2 Extract displacement interface v = Bu 3 Set contact force λ = 0 4 Initialization j = 0 5 while I > ε or j < j max do ← j + 1

17 20 i C ← i C + 1 21if i C = 1 then 22 for 24 B∆u 26 Corrected displacement field u = u + ∆u 27 Search 30 for ν 2 relaxations do 31 for 36 end 37 ComputeLATINLATINLATINFigure 14 :

 201222426273031363714 Figure 14: Convergence rates of the algorithm proposed vs. LATIN method: (a) Case A -(b) Case B -(c) Case C. Convergence diagrams are plotted with errors evaluated at each relaxation (coarse problem is considered costless). As no analytic solution exists, the reference is a solution evaluated with a precision of ε = 10 -12 . The number of SVD vectors defining the reduced order model is studied. 25

ν1 = ν2 = 1 ν1 = ν2 = 2 ν1 = ν2 = 3 ν1 = ν2 = 4 Figure 15 :

 123415 Figure15: Convergence rates of the algorithm proposed vs. LATIN method for case A. Convergence diagrams are plotted with errors evaluated at each relaxation (coarse problem is considered costless). The number of relaxations before (ν 1 ) and after (ν 2 ) coarse stage is studied.

Table 1 :

 1 Simulation parameters.

	Though frictional contact problems lead to non-smooth solutions, the represen-
	tation of a solution in the space-time domain allows to define a reduced basis

Table 3 :

 3 Algorithm parameters used for previous different cases.

Acknowledgement We gratefully acknowledge the French National Association for Research and Technology (ANRT, CIFRE grant number 2011/0705).