Identification of Contact Pressures in Two and Three-Dimensional Solid Bodies from Cauchy Data

This note deals with the identification of contact pressures in two and three-dimensional elastic bodies via two approaches relying on domain decomposition using electrostatic measurements. These approaches consist in recasting the problem in terms of primal or dual Steklov Poincar ́e equations. The numerical performances of these formulations are compared. The proposed methods are applied to some inverse problems: the first application deals with the identification of a Hertizian contact pressures distribution, the second deals with the identification of the indentation pressure of a heterogeneous solid, and the third one with the identification of boundary data at the interface of a bonded structure. value decomposition to solve the same problem numerically. A related inverse problem which allows for interior displacement measurements and inter-facial crack detection has been investigated by Huang and Shih [7]. Weikel et al. [8] have proposed an alternating iterative algorithm in order to reconstruct an internal planar crack laying on an a priori known internal surface inside a three-dimensional elastic body from over determined electrostatic boundary data on the outer surface. Numerical experiments for the identification of internal cracks in a three-dimensional elastic body using the primal and dual formulations of the Steklov-Poincar ́e equation are recently carried out and their numerical performances are compared [9]. A general framework for the different approaches (primal, dual and mixed) is presented presentad by Azaiez et al. [10]. In this work, the Steklov-Poincar ́e method is applied to the linear elastic data completion problem. In section 2, the Cauchy problem is presented in the context of linear elasticity. In section 3 this problem is recast in condensed form that we will refer to as the Cauchy-SteklovPoincar ́e problem, which leads to the Cauchy-Steklov-Poincar ́e equation acting on the boundary of the unknowns. In section 4 we present the Dirichlet-to-Neumann algorithm and we show that it can be interpreted as a preconditioned Richardson procedure for the CauchySteklov-Poincar ́e equation. The numerical procedures are presented in sections 3 and 5 and the results obtained by FEM discretization of the problems are presented in sections 7 and 8. The methods are used to solve applications borrowed from engineering mechanics in 2D and 3D frameworks: the identification of contact pressure between two elastic bodies, identification of the indentation pressure in a two layer solid and the identification of boundary data at the interface of a bonded structure. The Cauchy Problem in Linear Elasticity Let Ω denote a bounded domain in R2 or R3 with regular boundary Journal of Applied Mechanical Engineering J o u r n al o f A pp lied ical Eninee r i n g


Introduction
Physical phenomena are often governed by partial differential equations, which need an essential set of data to solve them. In linear elasticity, these data are: the geometry of the solid, the mechanical properties of the materials and the boundary conditions. However, in many industrial applications, some of these data are unknown and have to be identified. This leads to an inverse problem whose resolution requires over specified measured data. In this paper we focus on problems of boundary condition identification in linear elasticity. In this case, data measured on part of the easily accessible border are often available. However, contrary to the direct problem, two kinds of boundary conditions (e.g. displacements and tractions) are imposed on the same part of the boundary while no information exists on the remaining part of it. Hence, data completion consists in reconstructing the boundary conditions for the whole boundary of a domain by using the partially over specified measurements. This is the well-known Cauchy problem, which is ill-posed. The ill-posedness of inverse problems may concern the existence and/or the uniqueness of the solution, but their most critical feature is their instability: the solution, whenever a problem has one, is not continuous with respect to the data, i. e. small measurement errors in the data may dramatically amplify the errors in the solution. This is ill-posedness in the Hadamard sense [1]. The Cauchy problem pertains to this kind of inverse problem. Therefore suitable regularizing algorithms that are exempt from this ill-posedness phenomenon are required in order to solve the inverse problem correctly.
For inverse problems in elasticity, we refer to Bonnet et al. overview paper [2] and the huge amount of references therein. The Cauchy problem in linear elasticity was first studied by Yeih et al. [3]. In this paper, the existence and uniqueness of the solution are analyzed as well as the continuity of the solution with respect to the data. Other authors have proposed an alternative regularization procedure, namely the indirect fictitious boundary method, which is based on the simple or the double layer potential theory. The numerical implementation of the aforementioned method has been carried out by Koya et al. [4] who used the BEM and the Nystrom method for discretizing the integrals involved. Marin et al. [5] have determined the approximate solutions of the Cauchy problem in linear elasticity using an alternating iterative BEM that reduces the problem to solving a sequence of wellposed boundary value problems. Marin et al. [6] have used singular value decomposition to solve the same problem numerically. A related inverse problem which allows for interior displacement measurements and inter-facial crack detection has been investigated by Huang and Shih [7]. Weikel et al. [8] have proposed an alternating iterative algorithm in order to reconstruct an internal planar crack laying on an a priori known internal surface inside a three-dimensional elastic body from over determined electrostatic boundary data on the outer surface. Numerical experiments for the identification of internal cracks in a three-dimensional elastic body using the primal and dual formulations of the Steklov-Poincar´e equation are recently carried out and their numerical performances are compared [9]. A general framework for the different approaches (primal, dual and mixed) is presented presentad by Azaiez et al. [10].
In this work, the Steklov-Poincar´e method is applied to the linear elastic data completion problem. In section 2, the Cauchy problem is presented in the context of linear elasticity. In section 3 this problem is recast in condensed form that we will refer to as the Cauchy-Steklov-Poincar´e problem, which leads to the Cauchy-Steklov-Poincar´e equation acting on the boundary of the unknowns. In section 4 we present the Dirichlet-to-Neumann algorithm and we show that it can be interpreted as a preconditioned Richardson procedure for the Cauchy-Steklov-Poincar´e equation. The numerical procedures are presented in sections 3 and 5 and the results obtained by FEM discretization of the problems are presented in sections 7 and 8. The methods are used to solve applications borrowed from engineering mechanics in 2D and 3D frameworks: the identification of contact pressure between two elastic bodies, identification of the indentation pressure in a two layer solid and the identification of boundary data at the interface of a bonded structure.
Γ=∂Ω. The whole domain is assumed to be filled with a homogeneous linear elastic isotropic medium. It is assumed that Γ is splitted into two open subsets Γ c , and Γ i , Γ=Γ c ∪Γ i where Γ c , Γ i ≠∅ and Γ c ∩ Γ i =∅. In what follows, u(x) denotes the displacements field on Ω.
The local equilibrium equation is given by where σ is the stress tensor and f the volume forces. The strain tensor ε is given by These tensors are related by the Hooke's constitutive law, which is where λ and µ are the Lam´e constants of the material and I is the identity tensor.
Let n(x) be the outward normal vector at Γ and t(x) be the traction vector at a point x ∈ Γ defined by In the well-posed direct problem formulation, the knowledge of the displacement on a part of the boundary and traction vectors on another part of the boundary enables us to determine the displacement vector in domain Ω. Then, the strain tensor ε can be calculated from kinematic relation (2) and the stress tensor is determined by constitutive law (3).
If a part of the boundary Γ i is inaccessible and if it is possible to measure both the displacement and traction vectors on the remaining part of boundary Γ c , this leads to the mathematical formulation of a direct problem expressed as follows: where u  and t  are prescribed vector valued functions. This problem is ill-posed because of the formulation of its boundary conditions (4). It can be seen that boundary Γ c is over specified by prescribing both the displacement |Γ =  c u u and the traction |Γ =  c t t vectors, while boundary Γ i is underspecified since both the displacement |Γ = c u u and the traction |Γ = c t t are unknown and have to be determined. Then, this problem can be stated as follows: find (u, ) t such that a displacement field u(x) exists that satisfies: (u(x))n = ( ) in .
This problem, known as the Cauchy problem, is ill-posed in the sense that the dependence of u(x), and consequently of ( , t) u t is not continuous. Although the problem may have a unique solution, it is well-known that this solution is unstable with respect to small perturbations in the data on Γc. In this paper we propose to recover the lacking data by using the Steklov-Poincar´e algorithm introduced by Ben Belgacem et al. in the steady state thermal case in [11]. However, let us first introduce an operator acting on the boundary where data are unknown: the Steklov-Poincar´e operator which is very familiar in domain decomposition and recently introduced for the Cauchy problem of the Laplace equation by Andrieux et al. [13] and by Ben Belgacem et al. [11,12].

The Cauchy-Steklov-Poincar´e equation
To keep the notational complexity to a minimum let us remove x from the notations. Let λ denote the unknown displacement vector on Γ i . We consider both Dirichlet and Neumann elliptic problems obtained by duplicating the solution u into a couple of vectors (u N , u D ). The Cauchy problem (5) is then split into: Using the following notations: Equation (8) is called the Steklov-Poincar´e interface equation and S is the Steklov-Poincar´e operator. It is familiar in the domain decomposition framework [14] for the direct boundary value problem. More precisely, things happen as if vectors u D and u N were defined on two different domains with common boundary Γ i . In this case, the equation (8) expresses the Neumann transmission condition, but the (-) sign in S would be (+) in the domain decomposition formulation [14]. The (-) sign which appears in S is at the origin of the ill-posedness of the Cauchy problem. From the discrete point of view, the finite element discretization of S leads to the Schur complement matrix [14]. It corresponds to having all interior nodes eliminated by static condensation [15]. A numerical study of the Cauchy-Poisson problem, based on the Steklov-Poincar´e formulation is performed by Andrieux and Baranger [16].
We now continue with the analogy with domain decomposition and show how the Cauchy-Steklov-Poincar´e equation can be expressed, as in domain decomposition, in terms of the Dirichlet-to-Neumann problem.

Remark
A one dimensional example [17]: To illustrate how the illposedness of the Cauchy problem occurs in the Steklov-Poincare equation, we consider the problem of reconvering the end conditions (u,f) of a pre-tensioned string lying on a Winkler-type foundation, the end conditions at the other extremity ( , t)   u being given ( Figure 1). Denoting by F the tension of the string and by K the spring stiffness density of the foundation, the Cauchy problem can be written: find (u,f) such that there exists a vertical displacement field v verifying : . In this case S is a simple scalar function of k and L. It is easy to show that = − shkL chkL S chkL shkL which vanishes monotonically with respect to kL. This means that, as expected, the ill-posedness of the problem becomes more and stronger when the length of the string or stiffness of the foundation increases and when the pretension of the string decreases [17].

Interpretation of the steklov-Poincar´e equation:
Solving the steklov-Poincar´e equation is equivalent to the optimality condition of the first order associated to the energy-like error functional introduced in [18]. The proof is analogous to that done by Koslov and Maz'ya [19] for the Cauchy-Stokes problem.
The Dirichlet-to-Neumann algorithm, also borrowed from domain decomposition and introduced first by Ben Fatma et al. [20] to solve the Cauchy-Poisson problem, can be interpreted as a Richardson procedure applied to the Steklov-Poincar´e equation preconditioned by SD. The proof is the same than that used in domain decomposition [14] and that used [21] for Cauchy-Helmoltz problem.

The Dirichlet-to-Neumann solver for the Cauchy problem
When describing the Dirichlet-to-Neumann approach it should be noted that when the complete data are available on Γ, we have an over specied boundary value problem

u on n t u u on
This problem can be split into two well-posed sub problems with different boundary conditions. For one of them (Neumann/Dirichlet) while this is reversed for the other and (Dirichlet/Neumann) conditions are imposed on (Γ c , Γ i ) Solving the Cauchy system (5) is achieved when extension ( , ) u t makes û and u  coincide so the solution is Basically, the iterative method proposed for the Cauchy-Poisson problem and studied by Koslov et al. [20], is derived from these observations: starting from an arbitrary prediction of the Dirichlet condition (here the displacement vector u ) on Γ i , we add several corrections by solving alternately a Dirichlet on Γ c /Neumann on Γ i problem and a Neumann on Γ c /Dirichlet on Γ i problem, where at each iteration the appropriate boundary data are inferred from the solution computed in the previous step. More specifically, we construct a sequence of a pair of vectors from the following recurrence: given (0) D u , the following systems are solved for each k ≥ 0: The convergence of the alternating method toward the solution of the Cauchy problem and its stabilizing properties are established by Koslov et al. [20] for the steady state thermal case. In the linear elastic framework, no convergence result has been proved till now but the result of convergence established by Koslov et al. may be applied for any elliptic operator. When convergence is achieved, we may obtain ( , ) u t =(σ(u)n, u) on Γ i . By using straightforward computations, it can be established that the Dirichlet-to-Neumann scheme can be interpreted as a pre-conditioned Richardson procedure for the Cauchy-Steklov-Poincar´e equation. For this purpose, the Dirichlet-to-Neumann algorithm is rewritten, using the previous notations, as follows: Given λ 0 , the following systems are solved for each k ≥ 0

Remark
When implementing this approach, one has to keep in mind that the solution of the second PDE system in (10) is defined up to a rigid body motion. In order to relieve this problem, one can devote a wellchosen part of Γ c to set the Dirichlet boundary condition  u . This approach does not alter the generality of the work presented.

Algorithmic Strategy
The approximation of the problems 8 and 11 leads to the illconditioned linear systems: S h λ h =χ h for the SPP algorithm and D h µ h =Ψ h for the SPD one. S h and D h are the discrete Schur complement matrices. These linear systems are ill-posed since they are inverse problems. The vectors λ h and µ h are respectively the discretized unknown displacement and stress on Γ i . The Schur complement is too expensive to compute and ill-conditioned, an iterative procedure is henceforth in order. We use the GMRes algorithm, which is a popular iterative method for the solution of large linear systems. When starting vector λ 0 is zero, it generates a sequence of iterates λ 1 , λ 2 , ... such that λ k is in the Krylov subspaces K m (S h , r 0 )=span{r 0 , ..., S m−1 r 0 }, where r 0= χ h -S h λ 0 . Applying a Krylov method to S h λ h =χ h has been shown to have a regularizing effect [22]. In fact the Krylov subspace K m gives an approximation of the subspace generated by the m eigenvectors associated to the largest eigenvalues of S h [22], thus GMRes iterate λ k can be considered as an approximation of the truncated singular value decomposition solution [23].
The advantage of using the GMRES Algorithm is that it does not calculate the approximation of the solution at each iteration; only the final approximation λ m is computed.
To terminate the iterations of the GMRES algorithm, we use the L-curve criterion. The L-curve is the graph obtained by connecting consecutive points in the sequence [22,23]: In [23] Hansen shows that for discrete ill-posed problems it turns out that the L-curve always has a characteristic L-shaped appearance with a distinct corner separating the vertical and the horizontal parts of the curve. Calvetti et al. [24] shows that the optimal stopping criterion for the GMRes algorithm applied to liner ill-conditoned system corresponds to the vertex of the L-curve.
In this way, the L-curve clearly displays the trade-off between minimizing the residual norm and the solution norm.

Contact Pressures Identification
The example concerns contact identification on an inaccessible contact area. The data of the inverse problem are generated by solving the finite elements discretization of the direct problem. The thicknesses of meshes used for that are finer then that we used for solving the inverse problem.
Domain Ω is a square plate (1. × 1.) with a circular hole (R 1 =0.20225), where a fixed rigid disc (R 2 =0.2) is placed. Figure 2 shows the geometry   Table 1. When tractions are applied on the plate, it comes into contact with the lower part of the disc (Figure 2). The problem is to identify the contact pressure distribution and the displacements on the interface between the plate and the disc, by using over specified data provided for the external boundary. It should be noted here that the contact problem is nonlinear. However the data completion problem must be posed for that part of the domain where linear elasticity exists. The over specified data are generated by solving a direct problem using Hertz's analytical contact law. Here, we consider a frictionless contact so that only normal pressure is taken into account. Moreover, plane strain hypothesis is assumed (Figure 3).
The finite elements discretization of the equations (8) and (11) leads to ill-posed linear systems. This was expected as they are inverse problems. Since the Schur complement is too expensive to compute, we use the GMRes algorithm described above.
The results obtained by solving the corresponding Cauchy problem are the normal stress components and the displacements field on Γ i . Hence, the contact zone is the part of the boundary where the normal stress components are not null. When carrying out identification based on measurements, it must be kept in mind that measured data are subject to noise whose effects have to be studied. In this case, the data are synthetic, and therefore suffer from some errors (approximation error, roundoff error . . . etc). We added a noise generated by a MATLAB routine (randn) to the computational noise. The displacement measurements are polluted by a noise level at 5% and 10%. For the ND and SPP algorithms, the identification was carried out

Disc and plate Aluminium
Modulus of elasticity E=70000 M P A Poison coefficient ν=0.31

Friction coefficient µ=0
Load applied to the plate F=1e + 7 N/m   using over-specified data on the upper and lateral boundaries. It can be seen that the reconstructed stress does not approximate well the exact one: this is a consequence of the derivation operation. Although not presented here, we obtained an accurate numerical solution when the over-determined boundary Γ c and Γ i are complementary over ∂Ω (i.e ∂Ω=Γ c ∪ Γ i ). However, in practice, we cannot have stress measurements on the dirichlet boundary real-world problems. Figure 3 shows how the Dirichlet-to-Neumann algorithm is equivalent to a pre-conditioned Richardson procedure for the Cauchy-Steklov-Poincar´e equation as we have mentioned in section 4.
For SPD algorithm, the identification was carried out using over specified data only on the upper boundary which was discretized onto 65 nodes. We note that the numerical solution presented in Figures 5  and 6 approaches the exact solution as the noise level decreases and that the numerical results obtained by the SPD method are accurate and convergent with respect to decreasing the levels of noise added into the input data. We note that in order to preserve the stability of the method it is necessary to use a stopping criterion to cease the iterative process before the point where the errors start to increase due to the added noise. Figure 7 shows the error as a function of the number of iterations. It can be seen that the error e 2 decreases rapidly over the first few iterations but the rate of convergence decreases as the number of iterations increases. Thus the iterative process has to be stopped at a point where the error e 2 , obtained by comparing the calculated solution with the exact solution, stops decreasing. Moreover, if the number of iterations is large, due to the accumulation of the numerical noise, the error e2 starts to increase, this property is called semi-convergence.
As expected, displacements reconstruction is better than that for the stresses, particularly when the data are noisy. The reason is that the stresses are homogeneous with the displacements gradient and it is well known that the derivation is an ill-posed operation (the influence of noise is considerable). The identification is very satisfactory for free noisy data. For noisy data, the contact zone is well localized and the contact pressures are recovered correctly. However, some fairly significant oscillations appear on the free boundary. The table 2 provides relative error on displacements and number of iterations as a function of the number of degrees of freedom on Γ c . Γ c is the upper boundary of the plate (boundary subjected to the traction load).
We discretize Γ c respectively to 17, 33, 65, 129 and 257 nodes, the table 2 provides for each of the selected meshes, the number of iterations necessary to reach the solution for a free noisy data. The main conclusion we draw to is that provided that the data information on Γ c is sufficiently rich, the performances of the SPD algorithm are weakly sensitive to the amount of that information. Indeed, the quality of reconstruction process remains almost unchanged; no improvement is achieved by adding finite element nodes on Γ c .
The following point to investigate is how the SPD solver depends on the geometrical characteristics of Γ c : length and position. For this purpose we run five numerical experiments for different types of over specified boundary's (γ 1 : the upper and lateral boundaries of the plate, γ 2 : the upper boundary, γ 3 : the middle half of the upper boundary, γ 4 : the right half of the upper boundary and γ 5 : the left half of the upper boundary).
Two relevant indicators are recoding in the table 3: the number of iterations to reach the solution and the relative error on displacements. We note that the quality of reconstruction clearly suffer from decreasing the measure of Γ c , besides keeping the size of Γ c unchanged the position of it is very important regarding the quality of the information that can contain; we can see that γ 3 and γ 5 lead to nearly the same identification but better than that obtained with γ 2 (although of the same size) this is due to the information concerning the singularity in the stress tensor at the corner which is directly contained in γ 4 and γ 5 . The SPD method appears to be powerful and economic in comparison to the SPP one. In fact the primal schur complement matrix results from the static condensation on Γ i of the rigidity matrix, so it has the same distribution of singular values and thus inherits its initial ill-posedness, whereas the dual complement schur matrix can be seen as the inverse of the primal one so its condition number is obviously better. Therefore SPD is suitable for dealing with more complex data completion problems (e.g. 3D).

Dual Steklov-Poincare Algorithm for the resolution of the Cauchy problem in 3D linear elasticity
We investigate now the performance and accuracy of The SPD method through numerical three-dimensional examples where the measured data are extracted from the results of the direct problems: the identification of the indentation pressure of a heterogeneous solid and boundary data completion at the interface between two bonded elastic bodies.

Identification of indentation pressure on a two-layer solid
The identification of the indentation pressure on a composite solid and the stress on the interface is a major challenge in the design of composites. The example to be dealt with concerns the indentation of a two-layer solid composed of two different materials (steel and titanium) (Figure 8). By symmetry, the problem is reduced to that of a quarter of the solid. The problem is stated as follows: • Over specified data: displacements and surface traction are known on the rectangular area denoted by Γ m .
• Boundary Γ b gathering the four side faces: two have Dirichlet boundary conditions of symmetry, the two other are free from surface tractions.
• Boundary Γ u where the data are unknown includes the remaining faces (the bottom face and the remaining area of the top face).
Our aim is to identify the indentation area and the pressure distribution on the solid boundaries including the interface and the cross section by using the over specified data measured on Γ m . The measured data are generated numerically by solving the direct problem where the pressurized area (disk with radius R p ) and the expression of the pressure distribution as a function of the radius r are given by: The direct and inverse problems have been carried out using MATLAB software [25]. The geometry and mesh of the problem are displayed on Figure 8. The number of nodes on Γ u is 347, whereas the number of nodes on Γ m is 142. Hence, there are 345 × 6=2082 variables to identify for only 142 × 6=852 over specified known data. The measured data used are altered with noise that is a function of the mesh refinement, since it is extracted from finite element results of the direct indentation problem. But, in order to assess the efficiency of the method we add also some noise to the measured data.  Figure 9 show the map of the displacement norm ∥ u ∥L 2 obtained from the exact and the identified solutions with noisy data. Figure 10 show the map of the Von Mises equivalent stress obtained from the exact and the identified solutions. Figures 11 and   12 show the profiles of the Von Mises equivalent stress σ eq and the stress component σ zz on edges 1 and 2 and on the diagonal line; good reconstruction of the exact fields can be seen.           , distributed on the interface between the two materials of the solid. Here, too, a good correlation exists between the exact and identified results with reasonable accuracy, even for the interface tangential stress.

Boundary data identification of three-dimensional bonded structure
This application considers a bonded structure i.e. two bodies bonded along their common interface, by a thin adhesive layer. In the simplified models, the adhesive disappears, replaced by an interface transmission condition [26]. The problem has two planes of symmetry. Hence, only one quarter is modeled. Figure 18 shows the deformed shape and the distribution of the displacement field issued from the direct problem, defined as follows: • The first solid is a cylinder with radius r=0.248 m and a length L=3 m. The material is steel with a Young's Modulus E s =2.1 × 1011Pa and a Poisson coefficient ν s =0.34.
• The second solid is a rectangular box of dimension 1 × 1 × 2 m, which contains a cylindrical hole of radius 0.250 m. The   material is aluminum with a Young's Modulus E a =7 × 1010Pa and a Poisson coefficient ν a =0.27.
• The face x=0 is a plane of symmetry and is fixed in the x-direction: u=0.
• The face z=0 is a plane of symmetry and is fixed in the z-direction: w=0.
• On the circular face of the cylinder at z=1.5 m, a displacement is prescribed in the y-direction: D imp =−0.1 m.
A static finite element analysis was carried out using finite elements with the MATLAB Software. The result is displayed in Figure 18. Measured data (displacements and forces) are extracted on Γ m for use in the identification problem. Figure 18 shows the geometry used in the identification problem: the cylinder is ignored. The mesh used in this case has 1203 nodes, 4514 elements, 313 nodes on Γ u . We identify the stress and displacements using different localizations of the overspecified data Figure 18: • Γ ms : the top face of the solid with 82 nodes.
• Γ ml : the lateral side of the solid with 225 nodes.
• Γ l : the left lateral side of the solid with 118 nodes.
To stress the efficiency of the SPD algorithm we use also incomplete data: only tangential displacements on the over specified boundaries are used for the identification problem, the normal displacement is left unknown. Figures 19 and 20 show the identified displacement u y and stress σ yy along the edge defined by x=0 and y=0.25. We can see that the displacements are perfectly identified even when we use few and incomplete over specified data. The identified stresses are as well sensitive to the quantity of over specified data as to their localization; however, an acceptable identification remains possible.

Conclusion
In this work we presented three numerical methods for solving the Cauchy problem in the framework of linear elasticity. The methods proposed were applied in practical situations taken from engineering mechanics: contact pressure recovery, identification of the indentation pressure of a heterogeneous solid and boundary data completion in a bonded structure. The numerical results also suggest that the SPD algorithm is an accurate and reliable numerical technique for the identification of variables in elasticity in both two and threedimensional domains even if the area of the measured data is smaller than the area where the data are identified. The application of the methods to identify material properties and shape of the boundaries will be the subject of a forthcoming paper.