Rohit Pethe

Thomas Heuzé

Laurent Stainier

Variational h-adaption for coupled thermomechanical problems

Keywords: mesh adaption, variational approach, thermo-mechanics, nonlinear coupled problems, multi-physics

Introduction

In a number of transient problems (purely thermal, mechanical or thermo-mechanical), zones of high gradients of fields of interest evolve with time and loading. It is therefore interesting to use a dynamical mesh adaption algorithm to capture the solution in zones of high gradients in order to maintain the required precision. Many methods of mesh adaption have been proposed in the literature based on error-estimation. In these methods, the strategy is to adapt the mesh to minimize an error estimate, typically an upper bound, among all meshes of fixed size; or by recursive application of local refinement steps [START_REF] Verfürth | A Review of A Priori Error Estimation and Adaptive Mesh-refinement Techniques[END_REF] [START_REF] Ainsworth | A Posterior Error Estimation in Finite Element Analysis[END_REF]. But these methods have certain limitations. Rigorous estimates can be derived for linear constitutive models (for example elasticity), but it becomes more complex when non-linear constitutive models are used. Moreover, admissible fields need to be reconstructed [START_REF] Ladeveze | Error estimation and mesh optimization for classical finite elements[END_REF] [START_REF] Zienkiewicz | A simple error estimator and adaptive procedure for practical engineering analysis[END_REF]. In addition, standard error bounds require a certain regularity of the solution for their validity [START_REF] Ciarlet | Mathematical Elasticity[END_REF]). Therefore, it can be difficult and costly to use this approach for complex problems involving non-linear constitutive models and/or large deformation. In addition, methods based on global remeshing of the domain of interest require to transfer internal variables between meshes, which can lead to artificial diffusion of the latter unless specific methods are used [START_REF] Barthold | Error indicators and mesh refinements for finite-element-computations of elastoplastic deformations[END_REF] [START_REF] Brancherie | On a consistent field transfer in non linear inelastic analysis and ultimate load computation[END_REF].

Variational formulations allow us to express finite element problems as problems of minimization (or maximization) of an energy-like potential. This holds true for non-linear problems as well [START_REF] Dacorogna | Direct Methods in the Calculus of Variations[END_REF]. In some cases, the energy functional is evident, whereas, in some cases it needs careful formulation. For instance, in inelastic problems and dynamical problems, minimum principles can be obtained by careful time discretization [START_REF] Radovitzky | Error estimation and adaptive meshing in strongly nonlinear dynamic problems[END_REF] [START_REF] Ortiz | The variational formulation of viscoplastic constitutive updates[END_REF] [START_REF] Yang | A variational formulation of the coupled thermomechanical boundary-value problem for general dissipative solids[END_REF] [START_REF] Ortiz | Nonconvex energy minimisation and dislocation in ductile single crystals[END_REF]. In these cases, the energylike functional is incremental and incorporates the free energy, inertia and kinetics of material.

An alternative approach of mesh adaption for purely mechanical problems was recently proposed [START_REF] Mosler | Variational h-adaptation in finite deformation elasticity and plasticity[END_REF] [START_REF] Mosler | An error-estimate-free and remapping-free variational mesh refinement and coarsening method for dissipative solids at finite strains[END_REF], based on the variational approach of [START_REF] Ortiz | The variational formulation of viscoplastic constitutive updates[END_REF]. This technique uses an error indicator rather than an error estimator. In a variational approach, an energy-like potential is to be minimized (or maximized), and the gain in this scalar value associated to a given mesh adaption indicates the level of approximation, following the minimum (or maximum) criterion. No error estimates are then used at any stage of the algorithm. It allows mesh adaption in presence of large deformations and non-linear constitutive behavior. In addition, it was shown in [START_REF] Mosler | Variational h-adaptation in finite deformation elasticity and plasticity[END_REF]) that variational h-adaption could be combined with variational r-adaption, at least for hyperelastic behavior. Indeed, r-adaption would involve remapping in the presence of internal variables, and was not considered by these authors for dissipative behaviors. Recently, hp-adaptive energy minimisation has also been treated in [START_REF] Houston | Adaptive energy minimisation for hp-finite element methods[END_REF] for linear problems, the h-part of the adaptive process being close to former ones [START_REF] Mosler | Variational h-adaptation in finite deformation elasticity and plasticity[END_REF] [START_REF] Mosler | An error-estimate-free and remapping-free variational mesh refinement and coarsening method for dissipative solids at finite strains[END_REF]. In these, the authors addressed isothermal, steady state mechanical problems. The extension of the algorithm to thermo-mechanical problems raises some additional difficulties, which are addressed in the present work.

In this work, an h-adaption algorithm for problems in multi-physics is developed. The variational energy-like potential value is used to construct an error indicator and the variational principle itself drives mesh refinement and coarsening. Similar to [START_REF] Mosler | Variational h-adaptation in finite deformation elasticity and plasticity[END_REF], the algorithm presented is also based on a variational approach [START_REF] Ortiz | The variational formulation of viscoplastic constitutive updates[END_REF] but with an extended functional that admits a saddle point [START_REF] Yang | A variational formulation of the coupled thermomechanical boundary-value problem for general dissipative solids[END_REF][START_REF] Stainier | A Variational Approach to Modeling Coupled Thermo-Mechanical Nonlinear Dissipative Behaviors[END_REF]. This algorithm now accounts for heat conduction, transient thermal and thermo-mechanical coupling effects. For problems in multiphysics, the different physics have different temporal and spatial scales. Different meshes are used for each physics in order to account for its own spatial scale, these multiple meshes are adapted sequentially, the solution scheme relying on staggered algorithms [START_REF] Farhat | An unconditionally stable staggered algorithm for transient finite element analysis of coupled thermoelastic problems[END_REF] [START_REF] Armero | A new unconditionally stable fractional step method for nonlinear coupled thermomechanical problems[END_REF]. This permits to accurately capture the different spatial scales associated with each physics while maintaining the cost effectiveness of the approach. For constitutive behaviors involving internal variables, these are stored at integration points of each mesh, and no transfer (remapping) is necessary. Transfers between meshes only involve interpolation of external fields (displacement and temperature in this case). Though this paper focuses on 1-D problems in order to better illustrate and analyze the method, this algorithm can easily be extended to 2-D and 3-D problems using different subdivision schemes such as longest edge propagation path (LEPP) bisection algorithm of Rivara [START_REF] Rivara | Local modification of meshes for adaptive and/or multigrid finite-element methods[END_REF] [START_REF] Rivara | refinement algorithm suitable for adaptive and multi-grid techniques[END_REF] [START_REF] Bänsch | Local mesh refinement in 2 and 3 dimensions[END_REF] [START_REF] Rivara | Using longest-side bisection techniques for the automatic refinement of delaunay triangulations[END_REF] [START_REF] Rivara | Using longest-side bisection techniques for the automatic refinement of delaunay triangulations[END_REF]. In order to demonstrate this, first results from our ongoing work in 2-D is shown on a transient thermal test case. The edge bisection technique is used for refinement and coarsening the mesh. It enables to keep the same mesh topology so that a simple interpolation can be used between data located at integration points of the two meshes, avoiding a costly and diffusive projection of fields.

The structure of the article is as follows. In the second section, the variational formulation which is the base for our mesh adaption strategies is first presented. Then the mesh adaption algorithm is introduced. For the sake of clarity, the algorithm is first explained for steady state (thermal) problems followed by transient (thermal) problems and finally for strongly coupled (thermo-mechanical) problems. In the third section, studies on different test cases are presented. For each test case, cost and parametric analyses of the variational mesh adaption algorithm are carried out and comparison is made with respect to a uniform refinement mesh technique. At the end of this section, first results of a 2-D test case are presented to emphasize the extendability of the presented approach to 2D (and 3D) problems.

Formulations and methods

In this section, we recall the outlines of the variational formulation of coupled thermo-mechanical problems initially proposed in [START_REF] Yang | A variational formulation of the coupled thermomechanical boundary-value problem for general dissipative solids[END_REF]) (see also [START_REF] Stainier | A Variational Approach to Modeling Coupled Thermo-Mechanical Nonlinear Dissipative Behaviors[END_REF]). The presentation is specialized to the 1D case for simplicity, but the framework is actually very general, and is not limited to this case.

Steady state thermal problem

Assuming homogeneous properties, the local uni-dimensional thermal equilibrium reads

k d 2 T dx 2 + r = 0 ∀x ∈]0, L[, (1)
where k is the thermal conductivity, r is some external heat source and T is the temperature field.

In addition, let us consider the Dirichlet boundary conditions:

T (x = 0) = T 0 T (x = L) = T 1 (2)
This can be reformulated as a variational problem. The solution of the steady thermal problem is then found by minimizing the following convex potential:

φ(T (x)) = 1 2 L 0 k dT dx 2 dx - L 0 rT dx (3)
Indeed, the first variation of equation (3) yields

< DΦ, δT >= L 0 k dT dx d(δT ) dx dx - L 0 rδT dx = 0 ∀δT : δT (0) = δT (L) = 0, (4)
which is the weak form of equation ( 1). Since Φ is convex, any approximated, discretized temperature field T h will lead to the following inequality:

Φ(T h ) ≥ Φ(T analytical ) (5)
The variational mesh adaption algorithm presented here exploits directly this property plus the fact that, as a scalar energy, Φ is the sum of all the elementary contributions Φ(T h ) = e Φ e (T e h ) in a finite element discretisation framework. Indeed, the latter property allows us to define local patches of elements, refine them to decrease their values of Φ e , then add them to those of other patches. One local patch of elements may be isolated of the rest of the mesh by fixing the known temperature field on its boundary during the refinement procedure. Assuming a 1D thermal problem and patch made of one linear element, the mesh refinement procedure consists first in fixing the temperature at its end nodes, and in adding an additional node at the middle of this element. The procedure then amounts to minimize the value of energy-like potential Φ e over the patch. If the local improvement in the potential is considered as significant, the node is added to the global mesh. Otherwise, the mesh (in the local patch) is kept unchanged. This procedure is shown in figure 1. After looping over all the elements, a new global mesh is defined. However, since the boundary temperature of all the patches has been fixed for the refinement procedure based on the previous known temperature field, a global thermal problem needs to be solved on the refined mesh. This defines an iterative refinement procedure until the required precision is obtained. Thus, a given local patch may be refined or derefined, depending on the energy difference associated with mesh subdivision, the energy being itself related to the local regularity of the solution. Two tolerance parameters T ol r and T ol d associated with the refinement and derefinement bounds respectively are then introduced in the variational mesh adaption algorithm. A third tolerance parameter T ol 0 enables to stop the overall iterative procedure of adapted meshes. In one iteration of the procedure, two successive loops on patches are carried out in the sequence for the refinement and derefinement, followed by a new solution on the new mesh.

1: We begin with an arbitrary coarse mesh and solve our problem on that mesh. Get energylike potential φ G1 .

2: while

φ G 2 -φ G 1 φ G 1 ≤ T ol 0 do 3: φ G1 = φ G2 4:
Division of our full geometry Ω into different patches Ω i .

5:

for Ω i = F irst P atch to Ω i = Last P atch do 6:

Refine the current patch locally and solve a local problem on this small patch with the temperature field on the boundary of the patch imposed (given by the complete solution we calculated in earlier iteration).

7:

Calculate the local energy-like potential φ L2 .

8:

if φ L 2 -φ L 1 φ L 1 < T ol r then 9:
Add patch to list of unrefined patches. Refine the patch in global mesh.

12:

Add the patch to the set of refined patches. for first unrefined patch to last unrefined patch do 16:

Locally derefine the mesh on the patch or locally merge the patch with adjacent unrefined patch.

17:

Calculate the values of fields on deleted nodes by interpolation.

18:

Calculate local values of energy-like potentials with interpolated values φ L2 and calculated values φ L1 .

19:

if φ L 2 -φ L 1 φ L 1 < T ol d then 20:
Derefine the patch in the global mesh (or merge two unrefined patches in global mesh).

21:

end if

22:

end for

23:

Solve the problem again on the new mesh thus obtained and calculate φ G2 . 24: end while

Transient thermal problem

Transient thermal problems are governed by the heat equation, written in 1D as

c ∂T ∂t + ∂ ∂x -k ∂T ∂x = r ∀x ∈]0, L[ ∀t ∈ [0, t end ], (6)
plus a set of initial and boundary conditions. Here, c is the heat capacity per unit volume, k is the thermal conductivity, r is internal heat source density, and T is the temperature field. Assuming that k is homogeneous in Ω, we get

c ∂T ∂t -k ∂ 2 T ∂x 2 = r. (7)
The above equation can be written in a discretized time setting (backward-Euler) as follows

c ∆T ∆t -k ∂ 2 T n+1 ∂x 2 = r n+1 . (8)
Let's now define T ref as the reference temperature, and θ as the variation around it so that T = T ref + θ. Putting this in above equation, one gets

c T ref ∆θ ∆t - k T ref ∂ 2 θ ∂x 2 = r T ref . ( 9
)
The Helmholtz free energy can be defined as

W (θ) = - c T ref θ 2 2 , (10)
so that the entropy reads

η(θ) = - dW dθ = c θ T ref . ( 11
)
Introducing equations ( 10) and (11) in equation ( 9), one gets

∆η ∆t - k T ref ∂ 2 θ ∂x 2 = r T ref . ( 12
)
Let's now define the following incremental energy-like potential:

Φ(θ n+1 ) = L 0 (W (θ n+1 ) -W (θ n ) + η n (θ n+1 -θ n ))dx -∆t L 0 1 2 k T ref ∂θ n+1 ∂x 2 dx + ∆t L 0 r T ref (θ n+1 -θ n )dx (13)
It is important to note that the first variation of equation ( 13) gives the weak form of equation ( 12). The incremental variational principle reads:

θ n+1 = arg sup θ Φ(θ n+1 ) (14)
Although presented here in a heuristic fashion, it can easily be verified that this variational principle is a specific case of the more general variational formulation derived in [START_REF] Yang | A variational formulation of the coupled thermomechanical boundary-value problem for general dissipative solids[END_REF]).

In the case of transient problems, the mesh is adapted at the first time step as explained in section 2.1. From the second time step onward, the mesh adaption procedure is started with the final adapted mesh obtained at the previous time step as an initial mesh. This avoids a complete re-meshing of the domain at each time step. Two techniques may also be used in order to maintain a low cost of the algorithm without compromising the precision of the solution obtained.

The first one is to avoid to adapt the mesh at every time step. If domains of high gradients of fields of interest do not change much between two time steps, the mesh does not need to be adapted. After the refinement loop over the patches of the procedure, if no major changes have been found in the mesh, adaption at that time step can be bypassed.

The second technique pertains to the values of tolerance parameters in our algorithm. The current solution uses the final adapted mesh from the previous time step as initial mesh. Therefore, well-adapted meshes are required at the early time steps. Accordingly, small values of all of our parameters (T ol r ,T ol d ,T ol 0 and T ol u ) are hence used at the first time steps. Then, these values can increase with time steps till a given time, after which the values can be set constant.

Thermo-elasticity problem

Time discretization

Following [START_REF] Armero | A new unconditionally stable fractional step method for nonlinear coupled thermomechanical problems[END_REF], the thermoelasticity coupled equations read:

ρ ∂ 2 u ∂t 2 = ∂ ∂x (E(ε -αθ)) + ρb (15) c ∂θ ∂t = ∂ ∂x k ∂θ ∂x -αE ∂ε ∂t + r T ref ( 16
)
where ρ is the mass density, u denotes the displacement field, E is the Young's modulus, ε is the linearized strain tensor, α is the coefficient of thermal expansion, θ is the temperature variation such that T = T ref + θ. External loads ρb and r represent the body force vector and the external heat source density respectively. One also sets c = c T ref and k = k T ref , where c is the heat capacity and k is the thermal conductivity. An incremental energy-like potential, convex with respect to the displacement u and concave with respect to the temperature θ, defined at time step n + 1 is by:

Φ(u n+1 , θ n+1 ) = L 0 ρ 2β ∆u ∆t 2 - ρ β v n + 1 2 -β ∆ta n ∆u ∆t + W ( n+1 , θ n+1 ) -W ( n , θ n ) + η(θ n+1 -θ n ) -∆tψ(θ n+1 ) +∆t r(t n+1 ) T ref (θ n+1 -θ n ) -b(t n+1 )(u n+1 -u n ) dx -f (t n+1 )(u n+1 -u n ) x∈Γt -∆t q(t n+1 ) θ n+1 -θ n T ref x∈Γq (17)
as proposed in [START_REF] Yang | A variational formulation of the coupled thermomechanical boundary-value problem for general dissipative solids[END_REF] and [START_REF] Radovitzky | Error estimation and adaptive meshing in strongly nonlinear dynamic problems[END_REF]. The incremental problem then amounts to solving the saddle-point problem

u n+1 , θ n+1 = arg inf u sup θ Φ(u, θ) . (18)
Here, β and γ are the Newmark time integration parameters, u, v, a are the displacement, velocity and acceleration fields respectively, ψ(∇θ) is the heat conduction (Biot) potential, η the entropy, and W (ε, θ) the Helmholtz free energy. The boundary on which external forces f are applied is Γ t , whereas, Γ q is the boundary on which the heat flux q is imposed. The acceleration a and velocity v fields can be obtained from the displacement field u using Newmark's formulas:

a n+1 = u n+1 -u n β∆t 2 - v n β∆t - 1 2 -β β a n v n+1 = v n + (1 -γ)∆ta n + γ∆ta n+1 (19)
The Newmark's parameters are set at β = 1 4 and γ = 1 2 so that scheme can be unconditionally stable with respect to the time step size ∆t. The Helmholtz free energy can be given as follows:

W (ε, θ) = 1 2 Eε 2 -αEεθ - 1 2 c T ref θ 2 (20)
so that, the entropy η reads

η = - ∂W ∂θ = cθ T ref + εEα. (21)
In addition, the heat conduction potential can be given as

ψ(∇θ) = 1 2 k T ref ∂θ ∂x 2 .
(22)

Algorithm

An adiabatic staggered algorithm is used to solve the problem. This partition of the thermomechanical operator is known to preserve numerical stability [START_REF] Armero | A new unconditionally stable fractional step method for nonlinear coupled thermomechanical problems[END_REF].The mechanical part is solved by

u n+1 = arg inf u sup θ ad Φ ad (u, θ ad ), (23)
where Φ ad is built from the potential Φ, defined by equation ( 17), by removing the conduction (ψ(∇θ)) and the prescribed heat flux ( Γq ) terms. The stationary condition of u leads to the following system of equations

1 β∆t 2 [M ] + [E] [B] [B] T -[ C] {u n+1 } {θ ad } = [M ] β∆t 2 {u n } + [M ] β∆t {v n } + ( 1 2 -β)∆t{a n } + {b} + { f } -[ C]{θ n } + [B] T {u n } (24) Here, [M ], [E], [B],
[ C] and [ k] are the classical mass, elastic stiffness, coupling, normalized capacity and normalized conduction matrices respectively. Similarly, the thermal part is solved by

θ n+1 = arg sup θ Φ(u n+1 , θ). (25)
Therefore, one gets

-[ C] -∆t[ k] {θ n+1 } = -[ C]{θ n } + [B] T {u n } - {r} T ref + ∆t T ref {q} -[B] T {u n+1 } (26)
This staggered algorithm consists first in solving the mechanical problem by equation ( 24) assuming adiabatic conditions. Then, a thermal step at constant geometry is performed solving equation ( 26).

It has been shown [START_REF] Armero | A new unconditionally stable fractional step method for nonlinear coupled thermomechanical problems[END_REF]) that this staggered algorithm is unconditionally stable, provided a fixed time step [START_REF] Adam | Modélisation du comportement thermo-élasto-viscoplastique des métaux soumis à grandes déformations. Application au formage superplastique[END_REF]. This separation into two steps also allows us to use different meshes for the mechanical and the thermal fields. Therefore, in one time step, first, the mechanical mesh is adapted according to the adaption procedure explained in section 2.2 to get the mesh that best describes the displacement solution field using a minimum number of elements. During this, thermal fields are interpolated onto the mechanical mesh. Second, the thermal mesh is also adapted to get a mesh that best describes temperature field using a minimum number of elements. During this, the mechanical fields are interpolated on the thermal mesh. The adapted mechanical and thermal meshes at a given time step serve as initial meshes for the mechanical and thermal steps for the next time step.

For unidimentional problems, doing interpolation or projection of fields associated with one mesh onto another mesh won't make big difference in terms of cost and precision. However, doing interpolation rather than projection becomes important for 2D problems (and even more for 3D ones) since the gain of computational cost is significant, while it preserves a good accuracy.

Mathematical analysis

The variational functional Φ in equation ( 17) is equivalent to an H 1 norm. This is shown below for transient purely thermal problem for simplicity. Equation ( 13) can first be rewritten as follows

Φ(θ n+1 ) = L 0 - c T ref θ 2 n+1 2 + cθ n+1 θ n T ref - cθ 2 n 2T ref -∆t k 2T ref ∇θ • ∇θ dx. ( 27
)
Defining the following dimensionless field and variable θ = θ T ref , x = x∆tk c , and the functional J = -Φ( θ) cT ref , one gets

J( θ) = - Φ( θ) cT ref = L 0 1 2 ( θ -θn ) 2 + 1 2 ∇ θ • ∇ θ dx. (28)
The functional J( θ) can be put on the form

J( θ) = 1 2 a( θ, θ) -l( θ) + c 1 , (29)
where

a( θ, θ) = θ 2 1 = Ω( θ2 + ∇ θ 2 )dΩ (30) l( θ) = Ω θ θn dΩ (31) c 1 = Ω 1 2 θ2 n dΩ (32)
denote the H 1 (squared) norm, a linear form and a constant depending on the solution at the previous time step, respectively. The associated weak form is obtained by minimizing the functional

J( θ): a( θ, δ θ) -l(δ θ) = 0 ∀δ θ ∈ V, ( 33
)
where V is the space of continuous real valued functions in the domain Ω, satisfying homogeneous Dirichlet boundary conditions. Let θh be the finite element solution and θ the exact one. The error in the potential J( θ) due to the finite element discretization reads

J( θ) -J( θh ) = 1 2 a( θ, θ) -a( θh , θh ) + l( θh ) -l( θ). ( 34
)
Provided the H 1 (squared) norm of the error reads

1 2 || θ -θh || 2 1 = 1 2 a( θ, θ) + 1 2 a( θh , θh ) -a( θ, θh ), (35)
the combination of equations ( 35), ( 34) and ( 33) yields:

1 2 || θ -θh || 1 = J( θh ) -J( θ). ( 36
)
The difference in the variational potential is directly related to the H 1 (squared) norm of the interpolation error. This result is also found in [START_REF] Houston | Adaptive energy minimisation for hp-finite element methods[END_REF], although shown differently.

Results and discussion

4.1 Steady state

Analytical solution

Consider the following particular boundary conditions and heat source density:

r = x m T 0 = 0 T 1 = 0 ( 37
)
where m is a constant. Therefore, the analytical solution of the problem for a bar of length L reads

T = L m+1 x (m + 1)(m + 2) - x m+2 (m + 1)(m + 2) . ( 38
)
The energy potential is given by:

Φ(T ) = k 2 L 2m+3 1 (m + 1) 2 (2m + 3) - 1 (m + 1) 2 (m + 2) 2 -L 2m+3 1 (m + 1)(m + 2) 2 - 1 (m + 1)(m + 2)(2m + 3)
.

(39)

Numerical solution

The algorithm is started with a coarse initial mesh of two elements. Figure 2 shows the analytical solution as well as the variationally adapted temperature solutions on each iterated mesh, computed with P1-finite elements. These solutions are obtained with m = 51, L = 10 m, the thermal conductivity and the cross-section area are set to unity. After few iterations and with less than 50 elements, the sharp temperature gradient generated by the large value of m is well captured. Refined elements have been introduced close to that sharp gradient, and few elements are sufficient to represent the remaining part of the solution.

Cost Analysis

In order to assess the usefulness of this algorithm, the error between the computed solution and the analytical one can be plotted as a function of the number of nodes of the mesh. Three cases are considered. In the first case, the plot is made for uniform refinement of the mesh. This can be used as a reference. In the second case, the error at each refinement iteration in the variational adaptive mesh algorithm is plotted with respect to the number of nodes of the mesh. However, since the mesh adaption is done in several iterations, a consistent comparison between a uniform refinement and the variational one should account for path of refinement followed during mesh adaption. One way to accomplish this is to account for a cumulated number of nodes associated with all the calculations performed during the mesh adaption process. Therefore, the error at each refinement iteration is also plotted with respect to the cumulative number of nodes. In the third case, a comparison is performed with the Superconvergent Patch Recovery method (Zienkiewicz and Zhu 1992a;Zienkiewicz and Zhu 1992b) (also denoted ZZ2 ). The above method provides the global error estimator

e q 2 = N i=1 e q 2 i , (40)
computed from elementary contributions defined in 1D as

e q 2 i = Ωi e 2 q k dx, (41)
provided

e q = q * -q h , (42)
where q h = -k dT h dx is the finite element reconstruction of the heat flux vector, which is here constant elementwise, and q * is the nodally reconstructed ZZ2 solution from patches (see (Zienkiewicz and Zhu 1992a)). Elements are split into two if the relative value of the local estimator e q 2 i / q h 2 i is lower than some tolerance, set here identical to that of the variational mesh adaption algorithm.

Figure 3 shows the adapted temperature solutions on each iterated mesh adapted with the ZZ2 error estimator as well as the analytical solution. A close path of mesh adaption is observed between the variational (figure 2) and the ZZ2 (figure 3) ones, though the final meshes differ. Figure 4a shows the L 2 norm of the error on the temperature field, and figure 4b shows some energy norm of the error. More precisely, the potential energy Φ (3) is computed with the difference between the approximated and the analytical temperature fields. From figures 4a and 4b, several points can be emphasized. First, the curves of adaptive meshing algorithms are below that linked to uniform mesh refinement, which is quite expectable. Second, the curves of adaptive meshing techniques plotted as a function of the cumulated number of nodes cross from above the uniform mesh one, showing as expected that there is a number of nodes beyond which adaptive remeshing techniques are more performant and more cost effective than a uniform mesh technique. At last, the variational mesh adaption technique appears slightly more performant than the ZZ2 one on the range of error computed, both for cumulated and non-cumulated nodes, and both for L 2 error and the energy norm. Remark 1. In the variational mesh adaption approach, the cost of the refinement step of the mesh essentially consists of the computation of local solutions on refined patches. These problems are very small and cheap ones, because they consist in evaluating the solution at a single node. For thermal analysis, it yields one linear equation per patch to be solved since the temperature is fixed on the patch boundary. The computation of the error indicator for the whole domain thus have a complexity of O(N ), where N is the number of nodes. Once the mesh has been adapted, a new solution is computed on the updated mesh, whose complexity is O(N 3 ) if a direct solver is used. Hence, the cost of the refinement step remains far smaller than the computation of the solution.

Remark 2. The nodally reconstructed heat flux vector q * in the ZZ2 approach requires the solution of a small system of linear equations on each patch, whose dimension equals that of the polynomial basis used for the reconstruction. Though the dimension of the basis is usually small, the solution of a linear system on each patch yields a higher complexity of the refinement step than that achieved by the variational mesh adaption approach.

Parametric Analysis

The algorithm explained in section 2 exploits the additive property of the energy-like potential Φ by summing its local values over all elements. Improving the local value of Φ on a patch allows to reduce its global value on the mesh, and hence leads to a reduction of error. Three tolerance parameters have been introduced in the algorithm which influence its performance. However, the parameters T ol 0 ,T ol r and T ol d are not independent. In order to study the sole effect of each parameter, others are set to a constant value. The results of this parametric study are shown below.

Effect of T ol 0 : Figures 5 and6 show the influence of the tolerance parameter T ol 0 , while fixing T ol r and T ol d to a particular value. The parameter T ol 0 allows to decide when to stop the algorithm. It doesn't have any effect on the path followed. Therefore, one can observe that, as it decreases, the number of iterations increases. Therefore, the parameter T ol 0 should be selected such that the algorithm stops when a solution of a required precision (with respect to the current T ol r ) has been obtained. For example as shown in figure 6, since the T ol r and T ol d parameters are set at 0.5, a value of T ol 0 ranging between 10 3 and 10 2 is enough.

Effect of T ol r : Figures 7 and8 show the influence of the tolerance parameter T ol r , while fixing T ol 0 and T ol d to a particular value. The parameter T ol r drives the precision of the converged solution. As shown in figure 7, when the T ol r ranges between 5 × 10 -3 and 10 -2 the error of the converged solution is of the order of 10 -3 . In the graphs, one can observe that the algorithm carries out a few more iterations after convergence. This occurs because the value of T ol 0 has been set to a constant value in order to study the sole effect of T ol r . Whereas in normal circumstances, the value of T ol 0 is changed according to that of T ol r , so that the algorithm stops immediately after reaching convergence.

Effect of T ol d :

The results are shown in figures 9 and 10. Convergence and stability of the algorithm depend on the parameter T ol d . The latter should be less than or equal to T ol r , otherwise the algorithm will keep on refining and derefining the same patch entering in an unending loop. For example in figure 10, T ol r is fixed to 10 2 . All the curves that correspond to the values of T ol d less than or equal to 10 2 converge to the solution, whereas all the other curves diverge. This effect can also be observed in one of the curves in figure 8.

Improved Algorithm

In problems involving sharp gradients of the main field, many iterations of this iterative adaption process may be performed before convergence occurs, particularly if the initial mesh is coarse. Hence, it could be interesting to accelerate the refinement procedure by dividing an element in more than two elements. An application of this idea is shown in figures 11 and 12. The refinement procedure is so that:

1: if φ L 2 -φ L 1 φ L 1 > T ol u then 2:
Subdivide 1 element in 4 elements, where T ol u > T ol r .

3: else if T ol u > φ L 2 -φ L 1 φ L 1 > T ol r then 4:
Subdivide 1 element in 2 elements. 5: else if

φ L 2 -φ L 1 φ L 1 < T ol d then 6:
Consider derefinement. Parametric Analysis: Figures 13 and14 show the results of a parametric analysis carried out for the parameter T ol u while keeping T ol r constant at 10 4 . All algorithms give equivalent results after convergence and at the beginning. However, there is a big difference in the path followed to reach the converged state.

Thermo-elasticity

Consider a bar with homogeneous Dirichlet thermal and mechanical boundary conditions at its two ends:

T (0, t) = T (L, t) = 0; u(0, t) = u(L, t) = 0 ∀t (43)
along with a sinusoidal initial velocity:

u(x, 0) = 0; v(x, 0) = sin πx L ; T (x, 0) = 0 ∀x ∈]0, L[ (44)
This test case has been introduced in [START_REF] Armero | A new unconditionally stable fractional step method for nonlinear coupled thermomechanical problems[END_REF]. With these conditions, the bar is expected to vibrate, though damped through thermal dissipation.

Numerical solution fields

An adiabatic staggered scheme is used for the solution, as well as the algorithm of mesh adaption explained in section 2.3.2. The time step is set at 1 second in this test case. The problem is solved on a very fine mesh (4097 nodes) to obtain a reference solution, which has also been compared with the results obtained in [START_REF] Armero | A new unconditionally stable fractional step method for nonlinear coupled thermomechanical problems[END_REF] to ensure correctness. 15a, 15b, 15c, 15d, 15e and 15f show the displacement and temperature fields at times 1, 50 and 301 seconds respectively. In some parts of the bar, the algorithm has instroduced more nodes even where the solution field does not vary much. Two items may explain this behavior. First, our criterion for the mesh refinement is not directly related to the smoothness of the solution profile, but to the value of the energy-like potential. Recall that its value also allows to account for the variation of the solution field with respect to time, hence in a sense it generalizes the quasi-static approach of the ZZ2 error estimator. Second as explained in section 2.2, the mesh is not adapted at each time step in order to achieve a better cost effectiveness. Therefore, more nodes are required to capture the solution at different time steps for which the same mesh is used. This maintains the accuracy of the solution field and also the cost-effectiveness of the algorithm. However, it is evident from these figures that a very good solution field is captured at all the time steps.

Cost analysis

A cost analysis is performed on both thermal and mechanical meshes. The L 2 error of the displacement field and of the temperature field are computed on the mechanical and the thermal meshes respectively. The results are shown in figures 16a, 16b, 16c, 16d, 16e, and 16f. It is evident that the introduced mesh adaption algorithm is almost always more cost-effective with respect to a simple uniform mesh for both thermal and mechanical meshes. Recall that obtaining cost-effectiveness and a good accuracy of the solution at the first time step of the calculation is crucial, because the following adapted meshes directly depend on the previous ones. For example, at time steps 50 and 301, no mesh adaption is carried out because the meshes used at the previous time steps are good enough to represent the solution.

In this test case, the solution fields do not vary sharply with respect to time and space. Therefore, the algorithm drives the mesh adaption at quite few time steps, from time 1 to 301 seconds. The mechanical mesh is adapted at only 4 time steps, whereas the thermal mesh is adapted only at 3 time steps. Heuzé et al. [START_REF] Heuzé | Benchmark tests based on the couette viscometer -II: Thermo-elasto-plastic solid behaviour in small and large strains[END_REF]) have extended the well-known viscometer test case to thermoelastic-plastic solid behaviors in small and large strains. In this test case, the sole mechanical part acts on the thermal part, so that the mechanical problem is solved independently of the thermal problem. The mechanical problem is first solved, followed by the thermal one taking into account the effect of the mechanical solution. The geometry of the problem is shown in figure 17. The gap between the two cylinders is discretized by a radial 1D mesh. Zero displacement is prescribed on the inner cylinder while a driven rotation is prescribed on the outer cylinder. Temperature of external and internal cylinders are fixed to zero. Therefore, the boundary conditions of the problem can be stated as:

Thermo-elasto-plasticity

u θ (r = a) = 0 u θ (r = b) = u θ (b) T (r = a) = 0 T (r = b) = 0 (45)
where a and b denote the inner and outer radii respectively, and u θ the curved arc length swept since the finite strain framework is assumed. However the following differences arise between the test case given in (ibid.) and the present one. A hyperelastic-plastic constitutive law is considered here, whereas a hypo-elastic-plastic constitutive law was used in (ibid.). The analytical solution developed in (ibid.) relies on certain assumptions, that is: dilatation effects are neglected, thermal and mechanical parameters are fixed independently of the temperature and additional terms linked to the objective derivative are neglected. The solution developed in small strains is extended to the large strains in a straightforward manner, but its validity remains bounded from above when the rotations and hence the objective derivative become important. We solve the problem on a very fine mesh (5000 elements) with the numerical data of [START_REF] Heuzé | Benchmark tests based on the couette viscometer -II: Thermo-elasto-plastic solid behaviour in small and large strains[END_REF] The reference solutions in displacement and equivalent plastic strains are plotted in figures 18 and 19 respectively. According to (ibid.), the thermal solution is valid once the viscometer is completely elastic-plastic. As seen in figure 19, at rotation of θ = 3 • of the outer cylinder, the viscometer is completely elastic-plastic. Therefore, a coupled mechanical problem is solved starting at a rotation of outer cylinder of 3 degrees, provided the initial temperature being given by the analytical solution at that rotation, and a rotation evolution of the outer cylinder prescribed so that the plastic crown radius varies exponentially in time (eq.( 35) of (ibid.)), consistently with eq.( 24) of (ibid.). The reference solution in temperature is shown in figure 20. As seen from figure 18, the displacement field does not vary much but the thermal field presents the interest of a strong temperature gradient close to the inner cylinder. Therefore, it will not be very interesting to use adaptive meshing technique on the displacement mesh. Therefore, we solve our problem by adapting only the thermal mesh and keeping the mechanical mesh constant. Mech Mesh,θ = 3 One can observe that the numerical solution is very close to the analytical solution. The small differences between the numerical solution and the analytical one can be attributed to the different formulations of the mechanical constitutive models in large strains adopted in these two solutions. However, it is harmless for the mesh adaption purpose we are interested in here.

θ =1 • θ =3 • θ =6 • θ =9 •
p eq θ =1 • θ =3 • θ =6 • θ =9 •
θ =3 • . θ =4 • . θ =5 • . θ =6 • . θ =9 • .
• . Ther Mesh,θ =3 • . Analytical,θ = 3 • . Mech Mesh,θ = 6 • . Ther Mesh,θ =6 • . Analytical,θ = 6 • . Mech Mesh,θ = 9 • . Ther Mesh,θ =9 • . Analytical,θ = 9 • .

Analysis

As seen in figure 20, the domains of interest of the solutions field (domains with high gradients of temperature) do not evolve much in time. Therefore, our algorithm adapts the mesh only at the first time step, and then decides to use the same mesh for the following time steps. 23a, 23b, 23c, and 23d show the L 2 error of the adapted meshes. One can observe that mesh adaption has taken place only at a rotation of 4 degrees, all other time steps use the same mesh. The adaptive meshing still appears more economical than a uniform mesh.

Bidimensional transient thermal test case

In order to demonstrate the extendibility of this algorithm to 2D problems, the results of one test case obtained from our ongoing work in 2D is here presented. Let's consider a rectangular geometry whose boundary temperature is prescribed to zero. An external heat source is introduced in the domain, which follows a circular path in time centered within the rectangle. The heating area at one instant is generated by an arc length of 1 degree and a length of 1m in the radial direction. Figures 24a,24b, 24c, 24d, 24e and 24f show the solution fields on the adapted meshes at different time steps. A strong mesh adaption is performed in this test case because the location of strong temperature gradients moves with the prescribed heat source. The mesh coarsening upstream from the heat source appears as efficient as the mesh refinement where the heat source is located. It is therefore evident that this adaption strategy works well and is more cost effective than using a simple uniform mesh.

Conclusion

In this paper, a strategy for mesh adaption based on a variational approach for multiphysics problems has been proposed, in particular attention has been paid to thermo-mechanical problems. The variational approach uses an error indicator to adapt the mesh, based on the optimality property of the solution. The geometry is divided into patches and according to the level of improvement of the local value of an energy-like potential, refinement or de-refinement of the patch is performed. This strategy was first tested on simple steady state and transient thermal problems, for which a complete parametric analysis was performed. The effect of each parameter was studied and the strategies of selection of these parameters were discussed. Then the algorithm was applied to a strongly coupled problem of thermo-elasticity, using an adiabatic staggered algorithm and two different meshes for the thermal and the mechanical solution fields. Finally, the strategy was successfully tested on a more complicated thermo-elasto-plasticity benchmark test case. In this weakly coupled problem, the mechanical solution field was calculated on a fixed mesh, whereas the thermal mesh was adapted.

In all these test cases, it has been demonstrated that the developed strategy is reliable, economical and more effective than using a simple uniform mesh. The first perspective of this method in current progress is to extend this strategy to 2-D and 3-D problems, using different subdivision schemes, for example, the mesh subdivision scheme by [START_REF] Rivara | Using longest-side bisection techniques for the automatic refinement of delaunay triangulations[END_REF] and [START_REF] Rivara | Using longest-side bisection techniques for the automatic refinement of delaunay triangulations[END_REF].

Figure 1 :

 1 Figure 1: Mesh adaption algorithm on 1D problem.

Figure 2 :

 2 Figure 2: Analytical solution and variationally adapted temperature fields.

Figure 3 :

 3 Figure 3: Analytical solution and adapted temperature fields with the ZZ2 error estimator.

  , nodes non-cumulated ZZ2 adaptative algorithm, nodes non-cumulated Variational adaptative algorithm, nodes cumulated ZZ2 adaptative algorithm, nodes cumulated Uniform Refinement (b) Energy norm

Figure 4 :

 4 Figure 4: Error norms of the energy-like potential as a function of the number of nodes.

Figure 5 :Figure 6 :

 56 Figure 5: Effect of the parameter T ol0, T olr and T ol d being fixed to 10 -4 (figure shows L2 error vs. the number of nodes).

Figure 7 :

 7 Figure 7: Effect of the parameter T olr, T ol0 is set at 10 -2 and T ol d at 10 -4 (figure shows L2 error vs. the number of nodes).

  Figure 8: Effect of the parameter T olr, T ol0 and T ol d being fixed to 10 -2 (figure shows L2 error vs. the number of nodes).

Figure 9 :

 9 Figure 9: Effect of the parameter T ol d , T ol0 is set at 10 -2 and T olr at 10 -3 (figure shows L2 error vs. the number of nodes).

Figure 10 :

 10 Figure 10: Effect of the parameter T ol d , T ol0 is set at 10 -2 and T olr at 10 -2 (figure shows L2 error vs. the number of nodes).

Figure 11 :

 11 Figure 11: Comparison of the L2 error of the temperature field between the original and improved algorithms.

Figure 12 :Figure 13 :Figure 14 :

 121314 Figure 12: Comparison of the energy error between the original and improved algorithms.

  Temperature field at time = 301 seconds.

Figure 15 :

 15 Figure 15: Displacement and temperature fields at different times.

Figures

  Figures15a, 15b, 15c, 15d, 15e and 15f show the displacement and temperature fields at times 1, 50 and 301 seconds respectively. In some parts of the bar, the algorithm has instroduced more nodes even where the solution field does not vary much. Two items may explain this behavior. First, our criterion for the mesh refinement is not directly related to the smoothness of the solution profile, but to the value of the energy-like potential. Recall that its value also allows to account for the variation of the solution field with respect to time, hence in a sense it generalizes the quasi-static

  Cost analysis of the thermal mesh at time=301 seconds.

Figure 16 :

 16 Figure 16: Cost analysis of the mechanical and thermal meshes different times

Figure 17 :

 17 Figure17: Geometry of the test case[START_REF] Heuzé | Benchmark tests based on the couette viscometer -II: Thermo-elasto-plastic solid behaviour in small and large strains[END_REF]).

Figure 18 :

 18 Figure 18: Displacement solution (Heuzé et al. 2014).

Figure 19 :

 19 Figure19: Equivalent plastic strain reference solution[START_REF] Heuzé | Benchmark tests based on the couette viscometer -II: Thermo-elasto-plastic solid behaviour in small and large strains[END_REF]).

Figure 20 :

 20 Figure 20: Temperature reference solution.

Figure 21 :

 21 Figure 21: Superposed plot of the analytical (Heuzé et al. 2014) and numerical equivalent plastic strain.

Figure 22 :

 22 Figure 22: Solution on the adapted mesh for different rotations θ of the outer cylinder.

Figures

  Figures22a, 22b, 22c and 22d show solutions on adapted mesh at rotations θ = 4 • , 5 • , 6 • , 9 • of the outer cylinder respectively. Figure21shows analytical and numerical plastic strain distributions. One can observe that the numerical solution is very close to the analytical solution. The small differences between the numerical solution and the analytical one can be attributed to the different formulations of the mechanical constitutive models in large strains adopted in these two solutions. However, it is harmless for the mesh adaption purpose we are interested in here.

  21

  Rotation of θ = 4 • of the outer cylinder.

  Rotation of θ = 5 • of the outer cylinder.

  Rotation of θ = 6 • of the outer cylinder.

  Rotation of θ = 9 • of the outer cylinder.

Figure 23 :

 23 Figure 23: L2 error analysis at different rotations θ of the outer cylinder.

Figures

  Figures23a, 23b, 23c, and 23d show the L 2 error of the adapted meshes. One can observe that mesh adaption has taken place only at a rotation of 4 degrees, all other time steps use the same mesh. The adaptive meshing still appears more economical than a uniform mesh.

Figure 24 :

 24 Figure 24: Temperature field at different time steps.

  and use that solution as a reference one.

	0.16
	0.14
	0.12
	0.10
	0.08 u
	0.06
	0.04
	0.02
	0.5 0.00	0.6	0.7	r	0.8	0.9	1.0


Cesar Franck

Marle Franck

A Network-based Contract Management Framework for Prioritizing Contractual Performance Progress Axis

Keywords: Contract Management, network-based modeling, network-based analysis, network-based decision, decision aiding

Contracts play a major role in today's projects and operations. Contractual performance is defined as the result of the contract and the process and organization that support the contract lifecycle. It is determined by several interdependent factors, albeit they are generally managed without considering these interactions. A DSM is built to first identify potential relationships, either causeeffect or correlation, between contractual performance factors. This will allow for an easier and better decision about progress axis priorities in the field of Contract Management.

Introduction

Contracts play a major role in today's projects and operations. The way they are designed, awarded, executed and closed may strongly influence the performance of both clients and contractors.

Contract

Some examples of definitions are given as follows.

"An agreement between private parties creating mutual obligations enforceable by law. The basic elements required for the agreement to be a legally enforceable contract are: mutual assent, expressed by a valid offer and acceptance; adequate consideration; capacity; and legality. Possible remedies for breach of contract include general damages, consequential damages, reliance damages, and specific performance" [START_REF] Kim | Contract | Wex Legal Dictionary / Encyclopedia[END_REF].

"A contract is an agreement giving rise to obligations which are enforced or recognized by law. The factor which distinguishes contractual from other legal obligations is that they are based on the agreement of the contracting parties" [START_REF] Mckendrick | Contract Law: Text, Cases, and Materials[END_REF].

"A contract is an agreement with the intent of creating one or several juridical obligations. It is a voluntary engagement, formal or informal, between two or more parties and recognized by law" [START_REF] Barnett | Contracts: Cases and Doctrine, Sixth Edition[END_REF].

Contract is thus a deliverable (document), designed, written, signed, executed, managed and closed by organizational entities / actors, following a specific process. At least two A network-based Contract Management Framework for prioritizing contractual performance progress axis DSM 2019 entities/actors are involved in a contract, often called Client and Contractor. More broadly, the contract implies the notion of stakeholders and their satisfaction at different moments of its lifecycle. A contract may be a unique, isolated operation, or may be integrated in a bigger one involving multiple contracts. In all cases, contracts have to be managed.

Contract Management

According to the Chartered Institute of Purchasing and Supply, Contract Management (CM) is "the process of systematically and efficiently managing contract creation, execution, and analysis for the purpose of maximizing financial and operational performance and minimizing risk" (CIPS 2007).

It is important to emphasize that CM is an organizational process and, for this reason, can be different from one organization to another. However, the need for CM is shared by any kind of organization, public and private, in order to reduce costs, deal with regulatory requirements, and more broadly improve performance.

Contractual Performance

If we consider contracts as the tangible output of an organizational process (as seen in the previous paragraph), "contractual performance" is defined by WillBe Group as the combination of three elements:

-Performance of the output (contract): the ability of the contract to satisfy the parties involved, -Performance of the Contract Management process: the way of obtaining the output (efficiency), -Performance of the organization that performs the Contract Management process.

The aim of a Contract Management Framework (CMF) is to provide a clear vision of where and how much progress has to be done, in the aim of improving (or maintaining or restoring) global and long-term contractual performance of the entity. As described in the rest of the paper, the first part of our research work consists of building a performanceoriented CMF. A perspective will be to relate this to classical risk management, considering the risk of non-contractual performance.

The remainder of the paper is as follows. Section 2 briefly introduces a survey about existing (and available) Contract Management Frameworks. This helps to build the research question, exposed in Section 3. Then, Section 4 introduces the research approach, which is executed in Section 5 with a first model of interdependent contractual performance factors. Section 6 will discuss about this first version and will draw some conclusions and perspectives.

Background

Contract Management Frameworks (CMF) are diverse and numerous, albeit with similarities, notably in terms of nature of organizational entity that provides them. It is obvious that other CMFs exist without being public, so our literature review is restricted to available information. CMFs have been found for mainly 4 types of organizations: Most of them are process-oriented, describing either a sequential process, an iterative process or a cycle. They generally distinguish activities depending on their position before or after the contract signature:

- The next Section introduces their limits concerning complexity and introduces the academic question.

Limits and Research Problem

Complexity is recognized as influencing the way contracts are managed and their results [START_REF] Caldwell | Procuring Complex Performance in Construction: London Heathrow Terminal 5 and a Private Finance Initiative Hospital[END_REF][START_REF] Jarkas | Contractors' Perspective of Construction Project Complexity: Definitions, Principles, and Relevant Contributors[END_REF][START_REF] Perona | Complexity Management and Supply Chain Performance Assessment. A Field Study and a Conceptual Framework[END_REF][START_REF] Sivadasan | Advances on Measuring the Operational Complexity of Supplier-Customer Systems[END_REF]. For instance, the CIPS (Chartered Institute of Procurement and Supply) cites increasing contract volumes and complexity as reasons to recognize the importance and benefits of Contract Management (CIPS 2007).

In some existing CMFs, efforts have been done to link contract management elements, often processes and often with time-related links (sequential links, phases, with or without A network-based Contract Management Framework for prioritizing contractual performance progress axis DSM 2019 overlap, cycles, …). Some consequences of contract-related complexity are well-known, like for instance operational uncertainty, directly linked to financial risk [START_REF] Beg | Complexity, Contract Design and Incentive Design in the Construction Management Industry[END_REF]. However, some indirect effects, with possible propagations, or amplification, or retroaction, or variation (either of magnitude or even direction) of the links, are not well studied yet. It is thus a challenge to integrate the complexity of the contract management system, in order to improve contractual performance or increase risk of contractual performance.

The objective of the paper is thus to model interdependencies between factors of an existing Contract Management Framework, developed by the WillBe Group French consulting company. The justification of this question is that this model could be used for analyzing the CM factors network, in order to make easier and more reliable decisions about where to put priorities for performance improvement. The definition of these progress axis is today generally based on:

-models where factors are often considered as independent, -analysis based on an individual assessment of performance or maturity, -a decision-making process which does not simultaneously consider multiple parameters, like existing performance and importance of the factor, and required effort of the progress action.

All these reasons, operational need and literature gap, have justified the start of this research program, the first step being the network modeling project, presented here.

Research method

In this field, it is quite mandatory to be as close as possible of fieldwork, this is why we chose to use a constructivist approach [START_REF] Mäkinen | A Strategic Framework for Business Impact Analysis and Its Usage in New A network-based Contract Management Framework for prioritizing contractual performance progress axis DSM 2019 Product Development[END_REF], notably a Design Research Methodology-based approach, with alternation of prescriptive and descriptive study [START_REF] Blessing | DRM, a Design Research Methodology[END_REF][START_REF] Laurel | Design Research : Methods and Perspectives[END_REF].

The target is an overall decision aiding process inspired by [START_REF] Marle | A Structured Process to Managing Complex Interactions between Project Risks[END_REF], with a 3-step process, consisting of: 1/ modeling, 2/ analyzing a network, then 3/ making decisions considering specific complexity-oriented information obtained from previous steps. This paper focuses on step 1, the modeling of contractual performance factors network.

This step is based on the well-known DSM modeling technique, which emerged around thirty years ago and has spread, both in terms of volume of use and of nature of elements [START_REF] Browning | Design Structure Matrix Extensions and Innovations : A Survey and New Opportunities[END_REF][START_REF] Browning | Applying the Design Structure Matrix to System Decomposition and Integration Problems: A Review and New Directions[END_REF][START_REF] Eppinger | Design Structure Matrix Methods and Applications[END_REF][START_REF] Marle | Potential Applications of DSM Principles in Project Risk Management[END_REF][START_REF] Steward | The Design Structure System: A Method for Managing the Design of Complex Systems[END_REF]. Classically, the DSM offers a compact structure for visualization of complex networks. Moreover, it allows for systematic exploration of potential interdependencies by navigating through rows and columns. The main hypothesis was that, to our knowledge, no similar network currently exists, and we argue that DSM is more suitable to create this network between such immaterial and intangible factors. The assumption was more about appropriation by users, both for building and manipulating this network.

This modeling project relies on four components:

-1: internal experts' judgment (done) -2: literature review (done), -3: missions (ongoing), -4: clients' interviews (perspective at the moment of paper submission).

Two of them depend on external parameters and actors, with associated uncertainty (actual timeline and scope of missions, relevance of missions regarding the research question, number, adequacy and response time of interviewees).

One mission is currently ongoing with a possible link with the research. Authors hope that at the moment of the final submission of the paper, some additional information will have been provided by this fieldwork. Moreover, interviews will be structured and conducted with a refined version of the model, obtained by combining literature, internal experts' judgment and confrontation to the ongoing mission. This is planned for the end of 2019start of 2020. After this initialization phase, a continuous change management process will be applied.

Proposal of a first model of interdependencies between contractual performance factors

This Section introduces, first the initial version of the CMF developed by WillBe Group, second the modeling principles for interdependence identification, and third the refined network-based version of the CMF.

The initial 19-factor framework

Based on a CMF literature review and expertise of internal consultants, numerous factors had been selected, renamed, merged or broken down in order to get the initial 19-factor framework developed by WillBe Group: F1-Contract type, F2-Contractor choice, F3-Risk allocation, F4-Governance quality, F5-Quality of execution control, F6-Abusive claims from the contractor, F7-Client opportunism, F8-Respect of mutual obligations, F9-Interdependence between contractors, F10-Difference in clauses interpretation, F11-Transaction Costs, F12-Input data, F13-Economic equilibrium, F14-Deliberate underestimation of constraints, F15-Incompleteness of information, F16-Evolution of needs, F17-Renegotiation ex-post, F18-Alignment to evolution (market, technology), F19-Relational dimension.

The following paragraph illustrates how interdependencies are identified.

The seek for interdependencies between these factors

Interdependency exists when actions or events in one system affects or may affect elements in another system [START_REF] Marle | Managing Complex, High-Risk Projects[END_REF]. Several authors introduced different types of interdependencies, like for instance:

-Klir and Valach with linear, circular and recursive interdependencies, in terms of causality of phenomena [START_REF] Klir | Cybernetic Modeling[END_REF]. -Thompson and later Browning with the parallel, the sequential, the coupled and the contingent dependencies [START_REF] Browning | Applying the Design Structure Matrix to System Decomposition and Integration Problems: A Review and New Directions[END_REF][START_REF] Thompson | Organizations in Action[END_REF]. The notion of A network-based Contract Management Framework for prioritizing contractual performance progress axis DSM 2019 temporal relative position between the two elements is still present, so as the notion of one-way or reciprocal dependency and the notion of sure or potential dependency. -Le Moigne with interdependencies between three sub-systems: the decision system, the information system and the operational system (Le Moigne 1990). -Marle with the hierarchical link (typically found in WBS), the contribution link (meaning that one element contributes to the advancement of the other one), the sequential link, the influence link (if a decision or a change in element 1 may involve a change in element 2), and the exchange link if the two elements have an information flow [START_REF] Marle | Modèles d'information et Méthodes Pour Aider à La Prise de Décision En Management de Projets[END_REF]). -Worren with the commitment, the governance, the activity, the resource and the social interdependencies [START_REF] Worren | Organisation Design: Re-Defining Complex Systems[END_REF]).

An interdependency is thus described like in [START_REF] Marle | Managing Complex, High-Risk Projects[END_REF] through the characterization of which attribute(s) change(s), with a level of certainty (potential or sure) and of intensity (the number of attributes that change and the amount of change for each attribute). Attributes depend on the nature of elements, but classical attributes in contracts are related to cost, time, performance, and broadly speaking it describes the state of the element. This notion of attribute is not treated in this paper.

The contractual performance factors network v1.0

Based on initial list of factors and on previously described strategy for interdependency identification, the initial version of the contractual performance DSM has been built.

While looking for interdependencies, we noticed that a factor, F19, was linked with several ones, but in fact the links were more on a hierarchical nature (in terms of belonging to the same branch of a tree). We then considered this factor as a family of factors, this is why this updated version contains 18 factors. Moreover, in some cases, it is possible to give a sense to the positive or negative nature of a link. For instance, a good "Respect of mutual obligations" may decrease the necessity to run numerous controls and meetings, and then decrease "Transaction costs". However, other factors like "Contract type" and "Contractor choice" may be linked, albeit without any sense to an amplifying of retroactive effect. This is why for the moment we stay with positive values only, with this additional information on some specific cells.

Figure 2 hereunder represents the combined DSM, where each cell is the product of likelihood of interdependence and impact in case of interdependence. There are values from 1 to 9, 9 meaning that the interdependency has the strongest likelihood and the strongest impact. This development has been done right before submission of revised paper. Analysis is still ongoing, but for now it seems that there are enough levels to discriminate interdependencies. The reading convention for (i,j) cell is: factor Fi has an influence on factor Fj.

DSM 2019

Discussion and perspectives

Results of this first DSM-based modeling step are promising. The first lesson learnt is that factors list had to be adapted while considering interdependence between them. Questions arise about incorporating new factors, breaking existing ones down, merging other ones, and renaming some factors (for instance to improve clarity of definition and avoid overlaps). Another question arising from this initial work was the meaning of the cell. It is now after an iteration related to occurrence and impact of the interdependency.

Internal experts (from the consulting company) and current panel of external experts (from the client of the ongoing mission) recognize: 1/ the potential use of this model as an assistance to the prioritization of contract management progress actions; this may be helped by identification of key interdependences (strongest values in Figure 2), and possible chains and loops; 2/ the need to keep refining the network-based model and making it more robust, simultaneously generic but with the capacity to be quickly declined in specific contexts.

Future steps of this modeling project will consist in using maturity-based assessment of real organizations in order to: 1/ improve robustness of the network model by studying possible correlations between factors, depending on their maturity status; 2/ display factors list (and DSM) according to different sequences, the chronological vision being one of them (since there are factors that are clearly before or after contract signature, and some factors may be present all along the contract lifecycle); 3/ assess priorities for improvement, by combining maturity, importance and amount of work. This is inspired by previous work in automotive development context [START_REF] González | Measuring Project Maturity: Example in a French Automotive Organization[END_REF][START_REF] Gonzalez Ramirez | Assessing Project Maturity: A Case Study[END_REF]; 4/ group factors and the group actors according to the factors they are associated with and the links between these factors. This is inspired by other existing clustering works in nuclear construction and oil and gas sectors [START_REF] Marle | Organizing Project Actors for Collective Decision-Making about Interdependent Risks[END_REF][START_REF] Ventroux | Organizational Reshuffling to Facilitate Coordinated Decisions in Complex Projects[END_REF]. 5/ finally, we realized that the link between factors and activities should be more explicit, and an MDM-based approach could help that. We will also try to relate performance factors to outcomes and results of both contract and contract management (we may have one time a good result but with poor contract management and in fact poor efficiency and high risk of failure on a next contract). More important for us is the relation between network consideration and final performance, meaning that our main assumption in complex systems (and DSM modeling) is that if they are not managed with a mental representation and understanding of their complexity, then decisions may be inadequate and involve failures or at least major waste (of time, money, energy, trust, …).

  Upstream (or pre-award) activities: Preparing the business case and securing management approval; Assembling the project team; Developing contract strategy; Risk Assessment; Developing contract exit strategy; Developing a contract management plan; Drafting specifications and requirements; Establishing the form of contract; Establishing the pre-qualification, qualification & tendering procedures; Appraising suppliers; Drafting ITT (Invitation To Tender) documents; Evaluating tenders; Negotiating; Awarding the contract -Downstream (or post-award) activities: Change management; Service delivery management; Relationship management; Contract administration; Risk reassessment; Purchasing organization's performance and effectiveness review; Contract closure

Figure 2 :

 2 Figure 2: version 1 of the network-based CMF Next Section will discuss about results at date and next steps to reinforce the robustness of the model.

  Department of Housing and Public Works 2017; Municipal Association of Victoria 2018; Purchasing and Procurement Center of Malaysia 2017; Renewal SA 2015): the process is constrained by public regulations, so at least the award process should be public, but in some cases the entity publishes its overall process. -Consulting companies, that aim at assisting either clients or contractors, like Deloitte (Anderson and Mendel 2017) or PriceWaterHouseCoopers (PriceWaterHouseCoopers 2010). -Worldwide institutions, that aim at capturing best practices or just providing guidance and standards (World Bank 2017).

	CESAR Franck and MARLE Franck
	-Public organizations, at a national or local level, that audit procurement or
	purchasing-related activities (Department of Finance and Administration of
	Australian Government 2007; Department of Treasury and Finance Victoria State
	Government 2013; National Audit Office UK 2016; National Treasury
	Department of South Africa 2010).
	-Public organizations, at a national or local level, that are clients for contracts with
	generally private contractors (Department of Health of Government of Western
	Australia 2012;
	DSM 2019

DSM 2019

Contact: F. MARLE, Laboratoire Genie Industriel, CentraleSupélec, Université Paris-Saclay, 3 rue Joliot-Curie, 91170 Gif-sur-Yvette, FRANCE


Meriam Kilani

Franck Marle

Robustness of Analysis of Simplified Networks

Keywords: networks, network-based analysis, graph simplification, robustness, analysis sensitivity 1. Introduction

come

the same: does the simplified network behavior remain close to the complete one's? This is a question of robustness. The remainder of the paper is as follows. Section 2 introduces related work on network modeling, analysis and robustness. Section 3 presents research question and our approach to attempt to respond to it. Section 4 describes two robustness analyses to compare complete and simplified networks. Section 5 draws some conclusions and perspectives.

Related work

Network-based analysis

The composition and structure of the network may give information about its potential behavior. This is a static analysis, based on network snapshot, and the complementary analysis is the propagation analysis. It consists on dynamic simulation of an event, either the occurrence or disappearing of something in the network, and its consequences on the rest of the network (Koh, Caldwell, and Clarkson 2009;Oduncuoglu and Grebici 2010). This can be used for simulating the occurrence of desired or undesired events. In our case, it is a desired, albeit fictitious event, which is the suppression of one edge. This corresponds to network reduction ambition, in order to reduce calculation time (or prevent it to diverge to infinite) and increase human capacity to visualize and understand network behavior.

Network reduction

Possible strategies are to group nodes or edges, for different reasons [START_REF] Arenas | Size Reduction of Complex Networks Preserving Modularity[END_REF][START_REF] Blondel | Fast Unfolding of Communities in Large Networks[END_REF]. The groups are renamed or called super nodes or super edges. Elements may be grouped because they are strongly connected, or because they are connected to similar elements. However, we are not in this case. We do not want to reduce the number of nodes, but to test the suppression of some edges. It is more related to importance measures of edges, or robustness of network behavior.

Network robustness (to edge suppression, not node suppression)

Some works have been done to study network disruption, notably in transportation (Matisziw, Murray, and Grubesic 2007), telecommunications (O'Kelly and Kim 2007) or energy (Zio, Golea, and Sansavini 2012) networks. They analyze the impact of the loss of one edge (Jenelius, Petersen, and Mattsson 2006), and sometimes propose indicators, like Network Trip Robustness (Sullivan et al. 2010), or Network Robustness Index, or k-edge survivability index (Myung and Kim 2004), to characterize the capacity of the network to resist to the disruption [START_REF] Chen | Capacity Reliability of a Road Network: An Assessment Methodology and Numerical Results[END_REF]. We analyze here the sensitivity of removing edges to the occurrence of nodes knowing interdependencies between them. Removing edges will have an impact, the question will be how much.

Research question and approach

for instance. Calculation times may be exponential and saving for instance half of the edges may permit, either to run more simulations, or to expand network size. Moreover, it is important for human decision-makers to keep the possibility to visualize and understand the behavior of the network, even with the assistance of previous calculations and simulations. This paper thus tests the possibility to remove some elements in networks, mainly weak edges, in order to know if and how much precision is lost in terms of analysis. The research question can be formulated as follows: "Under what circumstances and conditions is it worthy to simplify network by removing weakest links?" This implies several sub-questions: 1/ What is the definition of weak link, and is it absolute or relative? 2/ What does "worthy" mean, or what is an acceptable loss of precision in the network analysis? 3/ How to characterize networks in order to know a priori if simplification under these conditions is reasonable or not? Our research work does initially aim at answering to these questions in familiar contexts, meaning that our history of studied networks delimitated the contexts. Basically, we considered two past projects with a direct focus on risk interdependency modeling. Several tests on real past modeled networks are made to illustrate this research. We are then in the first and second step of a classical Design Research Methodology [START_REF] Blessing | DRM, a Design Research Methodology[END_REF]Laurel 2003), formalizing questions from descriptive studies, and proposing and testing a new concept on these past case studies. We are not at the prescriptive study step yet, with an ongoing analysis and this simplification question at the moment when decision-makers are present. This could be done in 2020 or 2021 with a new decision-decision network modeling project recently started with an Oil and Gas company.

Comparison of risk rankings with complete network, simplified network and absence of network modeling

Two cases have been studied, one in the field of urban transportation system delivery, and one in the field of nuclear component production plant design and construction. In these contexts, project risks were initially identified, assessed and analyzed without considering potential interactions. This is classical risk management, called here BA for Basic Analysis.

Previous work had introduced a risk network modeling and analysis, called AA for Advanced Analysis. Simplified Advanced Analysis (sAA) is defined as the network-based analysis obtained with a simplified network, with interdependency values superior to a certain threshold.

The risk interdependencies were estimated by experts, using a discrete scale from 0 (impossible) to 10 (certain). For numerical calculation, expert judgments were transformed using a truncated negative exponential curve.

The tramway project

It was a complete turnkey project, aiming at delivering civil work, infrastructure, equipment and rolling stock. Budget was around several hundreds of million dollars, but the amount at stake was higher, since one secondary target was to prove competence of the company in order to get the 30-year operations and maintenance contract. In the initial Basic Analysis, 56 risks were identified, 13 of them being judged as critical. This was not an assistance to prioritization of risk mitigation actions. Some changes had been highlighted in the analysis of the risk network compared to the analysis of the risk list. 95 interactions have been identified, with interactions values from 1 to 10 (Fig. 1). These 95 interactions concerned 51 risks, 5 remaining isolated (in the limits of our model, not in reality of course). In AA, two significant changes occurred. First, rankings changed for 3 risks, R2, R7 and R10. The two most important changes were for R2, which became the most critical risk while considering network effect, and R10 which started from a "Negligible" category and became considered as an "Average" risk (27 th on 56). Second, even if the same rankings, the values changed and so gaps between risks. This means that there are fewer equal values, which may facilitate decision-making. For instance, 13 risks initially had the 0.381 value. After that, only 4 risks kept an equal value, those which did not change and remained at 0.381. The other ones changed, and thus the question becomes to know whether the gap is significant compared to the error or not. For R37, with an AA value of 0.416, this may not be the case, but for R43 with a value of 1.422 the gap is huge. This means that, with confidence, which was not the case with BA, decision-makers can put priority on R43 and R55 compared to R3, R16, R19 and R29.

In this study, we decided to filter values inferior to 7, which implies that only 26 nodes were connected by 27 links. 30 nodes became isolated. Figure 2 represents simplified network.

Fig. 2: simplified project risk network, with only edges superior to 7

Using the same analysis methodology, we were able to calculate a new refined criticality. This permitted to compare the latter (sAA) to the initial ones (from BA) and the refined ones (from AA). The first question is: knowing what should be detected (the two risks that see their ranking dramatically change, and the higher gaps between values which were initially close or even equal), does the simplified analysis capture this information or not? The analysis of the difference between the ranking with AA and the ranking with sAA, gives two information.

In the first place, the main ranking changes are captured by sAA, albeit it is partially done for R10 (which progresses of 7 places with sAA instead of 17 places with AA). We see the loss of the simplification, the fact that the change is in progress but not finished. Secondly, some minor changes are made for the first part of the ranking, meaning that biggest risks oscillate in terms of ranking before stabilizing to the final AA-based ranking. This does not

Part IV: New Methods and Algorithms DSM 2019 142 appear significant in terms of error. Thirdly, the biggest change is the fact that R2 becomes the most critical risk, and this is captured by sAA. So, it is good for R2 and R7, in progress but not so bad for R10, and minor changes for some of the biggest risks. The global answer to the first question is YES, but it seems to depend on the structure of the network. If strongest links are connected together and form the heart of the network, also called kernel, then there is little chance to miss a critical information. On the other hand, if two strong parts of the network are connected with multiple weak links, then this could not be captured by the simplified analysis. This was not the case here, but this remains as an interrogation about the network structure. The second question is: does the simplified network bring opportunity to decision-maker to visualize and understand the main drivers of the risk network? In this case, comparison between Figures 1 and2 gave a clear YES. We ran in parallel another test with another network, similar in nature (project risk network), but very different in structure, as described in following paragraph.

The nuclear component production plant

The second project is about the design, procurement, construction and ramp-up of a production plant in the nuclear area. Many risks are critical by themselves, but also critical by the impact they may have on further phases of the plant lifecycle. Just as in previous paragraph, basic and advanced analysis had already been done, and the object of the study was to compare them with the intermediary simplified analysis. From a 77-risk list with 12 possible likelihood values uses (from 0.1 to 0.9, with some intermediary values like 0.25, 0.35 and 0.75), we obtained a network with 74 connected nodes and 495 links, as illustrated in Figure 3. This network is quite difficult to read and to understand, however network-based calculations highlighted some risks due to their collective importance in the network behavior. The simplification of low values gave a network with 71 connected risks and 254 links (Figure 4). What is interesting here is that the network kept its global structure, meaning that high and low interaction values were quite distributed. This implies that it is not really easier to analyze the simplified network but may have the positive consequence to keep the same conclusions about network behavior. Secondly, this means that the highest interaction values (above 7 on a 10-level scale) represented half of the interactions. Decision-maker questioned this but kept the estimations since it is the nature of this kind of project to have strongly connected components. She decided not to change estimates but to assume the highly complex nature of this project. In terms of values and rankings, an additional indicator has been used, which is the fact that the risk ranking class changed or not. Due to specific management policy of the company, a classification by tens of risks was used, to get the top 10, top 20, and so on. So, we counted the number of risks that saw their class change. There are four cases:

-Case 1: A slight change in ranking (inferior to 5) and the same class -Case 2: A big change in ranking (superior to 6) and the same class -Case 3: A slight change in ranking (inferior to 5) and a class change -Case 4: A big change in ranking (superior to 6) and a class change 65% of risks were in case 1, 14% in case 2, 11% in case 3 and 9% in case 4. All risks in case 4 were not in Top 30, meaning that biggest changes were for low value risks only.

Only 3 risks of case 3 went out of Top 10 to integrate Top 20 (respectively 11 th , 12 th and 13 th places), meaning that, if final decision is made to allocate budget on Top 10 only, then we lose some information. However, Top 10 of the simplified analysis is still closer than the initial one obtained with basic analysis. The conclusion after this second test is that, if the structure of the network is quite distributed, then there is no significant loss of information when simplifying network. However, it remains sensitive to management practices and policies like for instance managing Top 10, Top 20, Top 30 with decreasing priorities (prioritizing blindly the 10 th risk above the 11 th ).

Conclusions

Regarding the initial question, "Under what circumstances and conditions is it worthy to simplify network by removing weakest links?", first results show that it is possible in some situations to make a simplification by removing lowest interdependency values. The simplified network analysis is closer to the complete network analysis than the basic analysis, often made without network consideration. This precision is still sensitive to the structure of the initial network and the position of weak links.

1/ What is the definition of weak link, and is it absolute or relative? Since we have as initial estimates discrete values, with a lot of ex-aequo risks, it is easier to use absolute thresholds. However, it could be combined if decision-maker wants to remove a significant percentage of edges (moving the absolute threshold until the relative percentage is obtained). Another perspective arose, which is to assess the location/nature of these weak links. Are these links buried with main blocks (or clusters) in a DSM or are they off diagonal (and outside the blocks)? Do these links represent feed forward or feedback information dependencies? 2/ What does "worthy" mean, or what is an acceptable loss of precision in the network analysis? Once again, this depends on decision-maker's intention and on company policies.

The two examples introduced here showed different practices with different perceptions of what is acceptable or not. In the first case, relative gaps and rankings were important, albeit in second case there could be different to be at the 10 th or 11 th place. Moreover, a sensitivity analysis could be done, in order to test different threshold values. It is equivalent to creating new intermediary snapshots of the network, with Figures 1 and2 (for the tramway project) as starting and ending points. The intermediary versions of the network, in this case, would keep the same structure. However, in a generic case, it is impossible to know a priori when and how the nature of structure may change. It is a good idea to test different values, also to help decision-maker to know what an acceptable threshold is (tradeoff between amount and precision of information about network behavior). This is a perspective for further work (not asked in this case).

3/ How to characterize networks in order to know a priori if simplification under these conditions is reasonable or not? This is more a perspective after these two initial descriptive tests. Many works exist to characterize network structure and a combination of real and fictitious situations could be used to know if Simplified Advanced Analysis is always closer to Advanced Analysis than to Basic Analysis. The idea is to identify patterns of networks and location of weak links in order to know whether weak links simplification has an influence on network behavior, i.e. on its structure. For instance, in this paper, the first case "Tramway project" has a funnel structure, with many sources leading to final effects. Removing weak links isolated a part of the network, while keeping its funnel structure. Similarly, for the second case "nuclear plant project", the network was very dense and connected, and kept this kind of structure. This means that we could (but further work will be done to confirm this assumption) assume that in such simple situations, the nature of the network is kept. However, the question is more difficult when a network cumulates different types of patterns, with a very connected part, another one with a funnel or cascade structure, and another one with multiple loops.

A last perspective is to analyze the sensitivity of proposed risk clusters to the removing of weak edges. A work is ongoing based on two other past case studies. Depending on the advancement of the study, some results may be presented during the conference presentation.

Fig. 1 :

 1 Fig. 1: initial project risk network (node labels are risk Ids, edge labels are interdependency strength)

Fig. 3 :

 3 Fig. 3: complete risk network of the nuclear production plant project

Fig. 4 :

 4 Fig. 4: simplified risk network of the project with half of the initial number of links

Modeled networks may contain hundreds or thousands of nodes and edges, which may decrease capacity to understand, analyze and make decisions. This is notably true when attempting to make calculations, either optimization-or heuristics-based, like clustering

Engineering Management 62(4). Jaber, H., F. Marle, L. A. Vidal, and L. Didiez. 2017. "Criticality and Propagation Analysis of Impacts between Project Deliverables." Research in Engineering Design. Jenelius, Erik, Tom Kristian Alex Petersen, and Lars-Göran Mattsson. 2006. "Importance and Exposure in Road Network Vulnerability Analysis." Transportation Research Part A: Policy and Practice 40(7):537-60. Retrieved (https://econpapers.repec.org/article/eeetransa/v_3a40_3ay_3a2006_3ai_3a7_3ap_3a537-560.htm). Koh, Edwin C. Y., Nicholas H. M. Caldwell, and P. John Clarkson. 2009. "Using a Matrix-Based Approach to Model Change Propagation." in International Design Structure Matrix Conference. Laurel, Brenda. 2003. Design Research : Methods and Perspectives. MIT Press. Retrieved (http://oskicat.berkeley.edu/record=b15960285~S1). Marle, F. 2014. "A Structured Process to Managing Complex Interactions between Project Risks." International Journal of Project Organisation and Management 6(1-2). Marle, F. and L. A. Vidal. 2008. "Potential Applications of DSM Principles in Project Risk Management." in Proceedings of the 10th International DSM Conference. Matisziw, Timothy C., Alan T. Murray, and Tony H. Grubesic. 2007. "Evaluating Vulnerability and Risk in Interstate Highway Operation." in Transportation Research Board. Washington DC, USA. Retrieved (https://trid.trb.org/view/801045). Myung, Young-Soo and Hyun-joon Kim. 2004. "A Cutting Plane Algorithm for Computing K-Edge Survivability of a Network." European Journal of Operational Research 156(3):579-89. Retrieved (https://www.sciencedirect.com/science/article/abs/pii/S0377221703001358). O'Kelly, Morton E. and Hyun Kim. 2007. "Survivability of Commercial Backbones with Peering: A Case Study of Korean Networks." Pp. 107-28 in Critical Infrastructure. Berlin, Heidelberg: Springer Berlin Heidelberg. Retrieved (http://link.springer.com/10.1007/978-3-540-68056-7_6). Oduncuoglu, A. and K. Grebici. 2010. "Quantitative Assessment Framework for Product Value and Change Risk Analysis in Early Design Process." Pp. 1-10 in Proceedings of the ASME 2010 International Design Engineering Technical Conferences & Computers and Information in Engineering Conference IDETC/CIE 2010. Retrieved May 29, 2012 (http://scholar.google.com/scholar?hl=en&btnG=Search&q=intitle:QUANTITATIVE+ASS ESSMENT+FRAMEWORK+FOR+PRODUCT+VALUE+AND+CHANGE+RISK+ANAL YSIS+IN+EARLY+DESIGN+PROCESS#0). Pollack, J., K. Costello, and S. Sankaran. 2013. "Applying Actor-Network Theory as a Sensemaking Framework for Complex Organisational Change Programs." International Journal of Project Management. Retrieved March 4, 2013 (http://linkinghub.elsevier.com/retrieve/pii/S0263786312001834). Sullivan, J. L., D. C. Novak, L. Aultman-Hall, and D. M. Scott. 2010. "Identifying Critical Road Segments and Measuring System-Wide Robustness in Transportation Networks with Isolating Links: A Link-Based Capacity-Reduction Approach." Transportation Research Part A: Policy and Practice 44(5):323-36. Retrieved (https://econpapers.repec.org/article/eeetransa/v_3a44_3ay_3a2010_3ai_3a5_3ap_3a323-336.htm). Zio, E., L. R. Golea, and G. Sansavini. 2012. "Optimizing Protections against Cascades in Network Systems: A Modified Binary Differential Evolution Algorithm." Reliability Engineering & System Safety 103:72-83. Retrieved October 11, 2012 (http://linkinghub.elsevier.com/retrieve/pii/S0951832012000415).

Contact: F. MARLE, Laboratoire Genie Industriel, CentraleSupélec, Université Paris-Saclay, 3 rue Joliot-Curie, 91170 Gif-sur-Yvette, FRANCE, franck.marle@centralesupelec.fr