Cross-scale calculation method for contact pressure and contact area of rough surface

Abstract

The invention discloses a cross-scale calculation method for contact pressure and contact area of a rough surface, and the method comprises the steps: constructing a contact stiffness model of a surface local region based on a microcosmic contact mechanics theory, and building a global weak form equation of a cross-scale contact problem according to a variational principle; secondly, establishing a discretization calculation model containing roughness contribution through a finite element method and a contact judgment criterion; and finally, a contact nonlinear iterative solution algorithm is established based on a Newton-Raphson algorithm and an improved Aikine acceleration method, and contact pressure distribution, a real contact area, a contact object displacement field and the like can be calculated. According to the method, the influence of surface microstructure parameters, contact object macroconfiguration and material parameters is considered, a conventional finite element method is expanded, contact cross-scale calculation between rough surfaces is realized, and theoretical and technical bases are provided for frictional wear, sealing analysis and the like of contact parts of a mechanical system.

Description

Technical Field

[0001] This invention relates to the field of multi-scale contact calculation methods for interfaces, specifically a cross-scale contact calculation method based on the combination of the finite element method and microscopic contact mechanics theory to calculate the contact pressure and actual contact area of ​​rough surfaces. Background Technology

[0002] Rough surface contact problems are widespread in modern engineering fields, such as bolted connections, wheel-rail contact, gear-rack mating, and precision drive systems. The safety, efficiency, and vibration characteristics of mechanical systems largely depend on the matching accuracy and contact stiffness of the contact pairs. Furthermore, the actual contact area is influenced by normal loads and surface morphology, thus determining various characteristics of the system, including friction, wear, sealing, and conductivity. Therefore, considering the microscopic morphology of the interface while simulating the static / dynamic mechanical behavior of macroscopic engineering components is an important cross-scale problem.

[0003] A common approach to handling rough contact problems is deterministic methods, such as the boundary element method (BEM) and the finite element method (FE). The BEM solves for the displacement and stress fields within the solid using boundary integral equations, but it struggles with nonlinearities (such as material and interface nonlinearities) and finite geometry problems. The FE, on the other hand, transforms variational contact inequalities into equality constraint problems using penalty methods or the Lagrange multiplier method. While deterministic methods require only basic solid mechanics assumptions and offer high accuracy, simulating rough contact is extremely time-consuming and computationally expensive; therefore, in engineering practice, the contact surface is often treated as a nominally smooth surface. Another approach is based on statistical microscopic contact mechanics. These models simplify complex rough surfaces into multiple non-interacting micro-protrusions with the same radius of curvature, assuming these micro-protrusions follow a Gaussian distribution in height, and using classical Hertzian theory to characterize the constitutive relations of individual micro-protrusions. However, these models are only applicable to nominally flat and single-curvature surfaces, making them unsuitable for complex engineering component surfaces.

[0004] With the rapid development of high-tech fields such as microdevices, high-speed rail, and aerospace, the above methods cannot meet the requirements of balancing computational cost and computational accuracy. Summary of the Invention

[0005] The purpose of this invention is to address the shortcomings of existing technologies by providing a bidirectional coupling cross-scale calculation method between macroscopic structural deformation and microscopic interface contact. This method combines the point-to-surface contact algorithm in the finite element method with microscopic statistical contact theory, enabling rapid calculation of the contact pressure and actual contact area of ​​engineering component interfaces.

[0006] To achieve the above objectives, the technical solution of the present invention is as follows:

[0007] A multi-scale calculation method for contact pressure and contact area on rough surfaces includes:

[0008] Step 1: Calculate the morphological parameters of the equivalent rough surface based on the contour information of the rough surface obtained from the measurement.

[0009] Step 2: Establish a contact model for a single micro-protrusion based on morphological parameters, and divide it into constitutive relationships for the fully elastic stage, the first elastoplastic stage, the second elastoplastic stage, and the plastic flow stage according to the indentation depth. Then, establish a normal contact stiffness model for the local region based on statistical methods.

[0010] Step 3: Construct the global weak form equations for the cross-scale contact problem based on the variational principle;

[0011] Step 4: Establish a discretized calculation model that includes the contribution of roughness deformation based on the finite element method and contact discrimination criteria;

[0012] Step 5: Establish an iterative solution algorithm for contact nonlinearity to solve for the contact pressure distribution and the actual contact area.

[0013] Furthermore, in step one, the micro-morphological parameters of the rough surface include the root mean square σ, the average radius of curvature of the micro-protrusions κ, and the micro-protrusion density η. These micro-morphological parameters are obtained through three surface moments m0, m2, and m4.

[0014] ,

[0015] ,

[0016] Where z(x) is the surface profile information obtained from actual measurement; x represents the spatial coordinates in the direction of surface profile measurement; E[·] represents the mathematical expectation; the superscript "*" represents the equivalent moments of the two surfaces, i.e. .

[0017] Furthermore, in step two, a contact model for a single micro-protrusion is established based on the morphological parameters obtained in step one. Constitutive relationships are then established for the fully elastic stage, the first elastoplastic stage, the second elastoplastic stage, and the plastic flow stage, based on the indentation depth. Finally, a normal contact stiffness model for the local region is established using statistical methods, expressed as follows:

[0018] ,

[0019] Where, ηA n Indicates the nominal contact area A n The number of micro-protrusions on it; E1 and E2 are the equivalent Young's modulus, where υ1 and υ2 represent the Poisson's ratios of object 1 and object 2, respectively, and E1 and E2 are Young's moduli; d is the indentation depth; ω ek is the critical indentation depth. Z and H * These represent the hardness coefficient and hardness, respectively. Let be the probability density function of the height distribution of the micro-protrusions.

[0020] Furthermore, in step three, the variational principle is used to establish the global weak form equations for the cross-scale contact problem.

[0021] ,

[0022] Among them, Π tot Total global energy; Π int Π ext_Ω and Π ext_Г These represent the virtual work done by internal strain energy and the virtual work done by external body force and surface force, respectively; Π con This is virtual work done by contact force; It is the elastic tensor; For strain tensor; For displacement field; For displacement field; For traction displacement; Indicates the boundary area; The internal region is represented; the virtual work of contact force is decomposed into two parts: the contribution of smooth matrix deformation and the contribution of roughness deformation, expressed as follows:

[0023] ,

[0024] ,

[0025] In the above formula, the matrix deformation contribution is obtained by the penalty function method constraint treatment, and the roughness deformation contribution is obtained by the contact stiffness model in step two; As a penalty factor; g n is the normal gap; k is the upper limit factor.

[0026] Furthermore, in step four, the region and boundary are discretized to form a series of nodes and meshes, and contact discrimination is performed according to the following criteria.

[0027] No contact: ,

[0028] Only the slightly raised part is in contact: ,

[0029] Micro-protrusions make contact with the substrate: ,

[0030] Traversing the contact boundary nodes forms a series of discretized contact pairs.

[0031] ,

[0032] ,

[0033] in, The variation of the gap; Let n be the contact force at the j-th discrete node; c1 This represents the number of discrete points where the micro-protrusion partially contacts the substrate; n c2 This represents the number of discrete points where only the micro-protrusions are in contact; The shape function matrix of the contact element;

[0034] The global stiffness matrix is ​​obtained by assembling the contact stiffness matrix and the structural matrix from the above equation.

[0035] .

[0036] Furthermore, if the condition is determined to be non-contact, then discretization is unnecessary;

[0037] If the micro-protrusion is in contact with the substrate, the discretization expression is as follows:

[0038] ,

[0039] If it is determined that only the micro-protrusion is in contact, then the discretization expression is as follows:

[0040] .

[0041] Furthermore, in step five, the Newton-Raphson iterative algorithm is used to solve the obtained discretized equations. An improved Aikinken acceleration method is used to reduce the number of iterations. When the normalized L2 norm of the residual force array is less than a given small amount, the calculation is considered to be converged and the iteration is terminated. After the iteration is terminated, the contact area and the corresponding nodal contact force are obtained. The nodal contact force is divided by the nominal area of ​​the element to calculate the contact pressure and distribution characteristics.

[0042] Furthermore, in step five, the actual contact area is calculated based on the size of the interface gap. The expression for the actual contact area of ​​the rough surface is:

[0043] ,

[0044] A e A ep1 A ep2 and A p These represent the actual contact areas in the elastic stage, the first elastoplastic stage, the second elastoplastic stage, and the plastic flow stage, respectively.

[0045] The beneficial effects of the above-described technical solution of the present invention are as follows:

[0046] This invention provides a novel cross-scale calculation method for contact pressure and contact area on rough surfaces. The proposed method comprehensively considers the influence of surface micro-morphology, macroscopic configuration of the contacting object, and material parameters, simplifying the complex micro-morphology of the surface. The proposed method solves the problems of cumbersome mesh generation, huge computational cost, and difficult contact convergence in traditional finite element methods.

[0047] This invention can be applied to calculate the pressure distribution and actual contact area of ​​parts in contact, obtain the displacement and stress distribution of parts under contact conditions, and determine the location of danger.

[0048] The present invention also provides an apparatus comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the steps of the above-described method.

[0049] The present invention also provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of the above-described method. Attached Figure Description

[0050] Figure 1 This is a schematic flowchart of the method for calculating the contact pressure and contact area of ​​a rough surface according to the present invention;

[0051] Figure 2 This is a schematic diagram of the macroscopic configuration of the contacting object in an embodiment of the present invention;

[0052] Figure 3 This refers to the contour height of the rough surface in the embodiments of the present invention;

[0053] Figure 4 This is the relationship between contact stiffness and normalized normal load in the embodiments of the present invention;

[0054] Figure 5 This refers to the mesh division and boundary classification in this embodiment of the invention;

[0055] Figure 6 This invention provides a contact discrimination criterion applicable to rough surfaces.

[0056] Figure 7 This is the iterative process for solving cross-scale contact nonlinearity in this invention;

[0057] Figure 8 This is a contact pressure distribution diagram under different normal loads in an embodiment of the present invention;

[0058] Figure 9 These are normal stress cloud diagrams of objects under different normal loads in embodiments of the present invention. Detailed Implementation

[0059] To enable those skilled in the art to gain a deeper understanding of the technical solutions of the present invention, the specific embodiments of the present invention will be described in detail and clearly below with reference to the accompanying drawings. Figure 1 The process of the present invention is illustrated. It should be noted that the described embodiments represent only one representative embodiment of the present invention, and not all embodiments.

[0060] Step 1: Calculate the morphological parameters of the equivalent rough surface based on the contour information of the contact surface obtained by measurement;

[0061] This embodiment describes the contact between a semi-cylindrical indenter (object 1) and a rectangular base (object 2), both having rough surfaces, with a macroscopic configuration as follows: Figure 2 As shown, the isotropic rough surface consists of numerous micro-protrusions with the same curvature κ and different heights z, whose root-mean-square height and distribution density are σ and η, respectively. Two-dimensional profile height data z(x) are obtained by measuring the two contact surfaces using a laser scanning digital microscope, as shown... Figure 3 As shown, the number of measurement points is N. s The surface moment is calculated using the height data z(x):

[0062] ,

[0063] The superscripts 1 and 2 represent contact surface 1 and contact surface 2, respectively. The surface moment of the equivalent rough surface is calculated using the following formula.

[0064] ,

[0065] Therefore, the morphology parameters of the equivalent rough surface are:

[0066] ,

[0067] σ can be calculated using formula (3) as 2.02 × 10 -6 m, η = 6.51 × 10 9 1 / m 2 κ = 2.01×10 -6 m.

[0068] Step 2: Construct a contact stiffness model for a local region of the rough surface based on the theory of microscopic contact mechanics;

[0069] Based on the morphological parameters obtained in step one, a contact model for a single micro-protrusion is established, and constitutive relations are divided into a fully elastic stage, a first elastoplastic stage, a second elastoplastic stage, and a plastic flow stage according to the indentation depth. When in the fully elastic stage, according to Hertzian contact theory, the normal contact load for a single micro-protrusion in contact is...

[0070] ,

[0071] The normal contact stiffness of a single microprotrusion is the differential of the normal load with respect to the indentation depth d, and can be expressed as:

[0072] ;

[0073] In the first elastoplastic stage, the normal contact load when a single micro-protrusion is in contact is

[0074] ,

[0075] Its contact stiffness is

[0076] ;

[0077] In the second elastoplastic stage, the normal contact load when a single micro-protrusion is in contact is

[0078] ,

[0079] Its contact stiffness is

[0080] ;

[0081] During the plastic flow stage, the normal contact load at the contact of a single micro-protrusion is:

[0082] ,

[0083] Its contact stiffness is

[0084] ,

[0085] in, ω is the equivalent Young's modulus; d is the indentation depth; ω e k is the critical indentation depth. Z and H * These represent the hardness coefficient and hardness, respectively.

[0086] The probability density function of the height distribution of micro-protrusions on a rough surface is a Gaussian distribution function.

[0087] ,

[0088] Therefore, by integrating over the surface height, the normal contact stiffness of the rough surface in the four deformation stages is obtained as follows:

[0089] ,

[0090] ,

[0091] ,

[0092] ,

[0093] Among them, K e K ep1 K ep2 and K p These represent the contact stiffness generated by the micro-protrusions during the elastic stage, the first elastoplastic stage, the second elastoplastic stage, and the plastic flow stage, respectively; ηA n Indicates the nominal contact area A n The number of micro-protrusions on the surface. The complete contact stiffness of a rough surface is composed of four superimposed parts.

[0094] ,

[0095] Substituting formulas (13), (14), (15), and (16) into formula (17), we obtain the complete calculation expression for the normal contact stiffness of the rough surface as follows:

[0096] ,

[0097] Where, ηA n Indicates the nominal contact area A n The number of micro-protrusions on it; ω is the equivalent Young's modulus; d is the indentation depth; ω e k is the critical indentation depth. Z and H * These represent the hardness coefficient and hardness, respectively.

[0098] In this embodiment, the material parameters are: GPa, H * = 3340 MPa. Figure 4 The relationship between the contact stiffness calculated by formula (18) and the normalized normal load shows that the contact stiffness increases with the increase of the normal load.

[0099] Step 3: Construct the global weak form equations for the cross-scale contact problem based on the variational principle;

[0100] The global weak form of the cross-scale contact problem consists of macroscopic configuration deformation and microscopic interface micro-protrusion deformation:

[0101] ,

[0102] Among them, Π int Π ext_Ω and Π ext_Г These represent the virtual work of internal strain energy and the virtual work of external body force and surface force, respectively. The virtual work of contact forces on the configuration matrix; Let represent the virtual work generated by the deformation of all micro-protrusions in the contact area. From the variational principle, its expressions are as follows:

[0103] ,

[0104] ,

[0105] ,

[0106] ,

[0107] ,

[0108] In formula (20) Indicates the internal region of the configuration. It is a fourth-order elastic tensor. For strain tensor, For displacement field, For the variation of displacement; in formula (21), f is the body force inside the configuration; in formula (22) Indicates the boundary of the configuration. For Newman's border The traction force on the surface; in formula (23) and Representing the penalty factor and the contact boundary, respectively. The gaps between the interfaces on the screen is the variational form of the gap; k is the upper bound factor.

[0109] Substituting equations (20), (21), (22), (23), and (24) into equation (19) yields the global weak form equation for the cross-scale contact problem.

[0110]

[0111] Step 4: Establish a discretized calculation model that includes the contribution of roughness deformation based on the finite element method and contact discrimination criteria;

[0112] This embodiment deals with a plane stress problem. Without loss of generality, a three-node triangular mesh is used to discretize object 1 and object 2 into a series of nodes and meshes. Their boundaries are then divided into Dirichlet boundaries, Newman boundaries, and contact boundaries, as shown below. Figure 5 As shown. Considering the surface roughness and its height distribution probability density function, the contact discrimination criterion is as follows:

[0113] No contact: ,

[0114] Only the slightly raised part is in contact: ,

[0115] Micro-protrusions make contact with the substrate: ,

[0116] Traverse the contact boundary nodes ( Figure 6 The proposed contact discrimination criterion applicable to rough surfaces is shown, forming a series of discretized matrix forms of contact pairs. If object 1 and object 2 do not contact, no discretization is required.

[0117] If the micro-protrusions of objects 1 and 2 are in contact with the base, then the discretization is expressed as equation (29), with the specific expression being:

[0118] ,

[0119] If objects 1 and 2 only have slightly raised parts in contact, then the discretization is expressed as equation (30), with the specific expression being:

[0120] ,

[0121] The contact stiffness matrix is ​​formed from the above equation. With structure matrix Assemble the assembly to obtain the global stiffness matrix

[0122] ,

[0123] in, and All are element contact stiffness, structural matrix From the element stiffness matrix Assembly and assembly

[0124] , ,

[0125] In the above formula, and These are the strain matrix and the elasticity matrix, respectively.

[0126] A discretized calculation model incorporating roughness deformation contributions is constructed using formula (32), specifically as follows:

[0127] ,

[0128] Where U is the discretized displacement matrix, F ext F int and F c These are the discretized external load array, the discretized internal force array, and the discretized contact force array, respectively.

[0129] Step 5: Establish an iterative solution algorithm for contact nonlinearity to solve for the contact pressure distribution and the actual contact area;

[0130] Equation (33) is equivalent to a convex optimization problem: finding a discretized displacement matrix U that minimizes G(U) while ensuring the contact criterion holds. The standard Newton-Raphson iterative algorithm is used to solve this optimization problem, updating the displacement matrix, load matrix, and contact stiffness matrix in each iteration. An improved Aikinken acceleration method is used to correct the displacement matrix in the Newton-Raphson iterative process, specifically:

[0131]

[0132] During the iteration process, the normalized L2 norm error of the residual force is used for convergence judgment, specifically:

[0133] ,

[0134] when If convergence is achieved, the calculation can be determined. Otherwise, it continues to iterate until the preset maximum number of iterations is reached. If convergence is still not achieved after reaching the maximum number of iterations, it is determined to be divergent and the contact settings and boundary conditions need to be rechecked. Figure 7 The iterative process of solving nonlinear contact problems across scales is demonstrated.

[0135] After the iteration terminates, the contact area and the corresponding nodal contact force can be obtained. The contact pressure and distribution characteristics are calculated by dividing the nodal contact force by the nominal area of ​​the element. Figure 8 The figure shows the contact pressure distribution under normal loads of 5, 50, 500, and 5000 N. As the load increases, the contact pressure distribution on the rough surface gradually converges to that on the smooth surface, and the contribution of matrix deformation to the total energy gradually increases. The displacement field, stress field, and other mechanical parameters of the contacting object are directly calculated using the standard finite element method. Figure 9 The normal stress field contour plots are shown under different loads.

[0136] Based on statistical principles, the actual contact area is further calculated using the size of the interfacial gap. The contact areas of a single micro-protrusion in the four deformation stages (fully elastic stage, first elastoplastic stage, second elastoplastic stage, and plastic flow stage) are respectively...

[0137] ,

[0138] ,

[0139] ,

[0140] ,

[0141] Combining the probability density function of the height distribution of the micro-protrusions on the rough surface in formula (12) with the Gaussian distribution function, the actual contact area of ​​the four deformation stages is calculated as follows:

[0142] ,

[0143] ,

[0144] ,

[0145] ,

[0146] Among them, A e A ep1 A ep2 and A p These represent the actual contact areas in the elastic stage, the first elastoplastic stage, the second elastoplastic stage, and the plastic flow stage, respectively.

[0147] Formulas (40), (41), (42), and (44) are superimposed to form a complete expression for calculating the actual contact area of ​​the rough surface.

[0148] .

[0149] The actual contact area of ​​each contact pair is calculated using formula (44), and the summation yields the actual contact area of ​​the entire contact region.

[0150] The above description represents the preferred embodiments of the present invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A method of cross-scale calculation of rough surface contact pressure and contact area, characterized in that, The method comprises the following steps: Step 1, obtaining the profile information of the rough surface by measurement and calculating the topography parameters of the equivalent rough surface; Step 2, establishing the contact model of a single micro-protrusion according to the topography parameters, and dividing the contact relationship into the constitutive relationship of the complete elastic stage, the first elastic-plastic stage, the second elastic-plastic stage and the plastic flow stage according to the indentation depth, and then establishing the normal contact stiffness model of the local area according to the statistical method; Step 3, constructing the global weak form equation of the cross-scale contact problem according to the variational principle; Step 4, establishing the discretization calculation model containing the roughness deformation contribution based on the finite element method and the contact criterion; Step 5, establishing the iterative solution algorithm of the contact nonlinearity, and solving the contact pressure distribution and the real contact area.

2. The method of claim 1, wherein, In step 1, the rough surface micro-topography parameters include the root mean square σ, the average micro-protrusion curvature radius κ and the micro-protrusion density η, and the rough surface micro-topography parameters are obtained through three surface moments m0, m2 and m4 , , where z(x) is the surface profile information obtained from actual measurement; x represents the spatial coordinate in the surface profile measurement direction; E[•] represents the mathematical expectation; the superscript "*" represents the equivalent moments of the two surfaces, i.e. .

3. The method of claim 2, wherein, In step 2, the contact model of a single micro-protrusion is established according to the topography parameters obtained in step 1, and the constitutive relationship is divided into the constitutive relationship of the complete elastic stage, the first elastic-plastic stage, the second elastic-plastic stage and the plastic flow stage according to the indentation depth, and then the normal contact stiffness model of the local area is established according to the statistical method, and the expression is ， where ηA n represents the number of microprotrusions on the nominal contact area A n ; is the equivalent Young's modulus, where υ1 and υ2 represent the Poisson's ratios of object 1 and object 2, respectively, and E1 and E2 are the Young's moduli; d is the indentation depth; ω e is the critical indentation depth; k Z and H * represent the hardness coefficient and hardness, respectively; is the probability density function of the microprotrusion height distribution.

4. The method of cross-scale calculation of rough surface contact pressure and contact area according to claim 3, characterized in that, In step 3, the global weak form equation of the cross-scale contact problem is established by using the variational principle , where Π tot is the global total energy; Π int , Π ext_Ω and Π ext_Г represent the internal strain energy virtual work, the external body force and surface force virtual work, respectively; Π con is the contact force virtual work; is the elastic tensor; is the strain tensor; is the displacement field; is the displacement field; is the traction displacement; denotes the boundary region; denotes the internal region; the contact force virtual work is decomposed into two parts, the smooth base deformation contribution and the roughness deformation contribution, and the expressions are ， ， In the above formula, the base deformation contribution is obtained by the penalty function method constraint processing, and the roughness deformation contribution is obtained by the contact stiffness model in step two. is a penalty factor; g n is the normal gap; k is the upper limit factor.

5. The method of claim 4, wherein, In step 4, the region and the boundary are discretized to form a series of nodes and grids, and the following contact criterion is used No contact: , Only micro-protrusion parts contact: , The microprotrusions are in contact with the base portion: , The contact boundary nodes are traversed to form a series of discrete forms of contact pairs ， , wherein is the variation of the gap; is the contact force on the jth discrete node; n c1 represents the number of discrete points where the microprotrusion is in contact with the base portion; n c2 represents the number of discrete points where only the microprotrusion is in contact; is the contact element shape function matrix; The global stiffness matrix is obtained by assembling the contact stiffness matrix and the structure matrix 。 6. The method of cross-scale calculation of rough surface contact pressure and contact area according to claim 5, characterized in that, If it is determined that there is no contact, no discretization processing is required; If it is determined that the micro-protrusion and the substrate are partially in contact, the discretization specific expression is ， If it is determined that only the micro-protrusion is in contact, the discretization specific expression is 。 7. The method of cross-scale calculation of rough surface contact pressure and contact area according to claim 6, characterized in that, In step 5, the Newton-Raphson iterative algorithm is used to solve the discretization equation obtained, wherein the improved Aitken acceleration method is used to reduce the number of iteration steps, and when the normalized L2 norm of the residual force array is less than a given small amount, it is determined that the calculation converges, and the iteration is terminated; after the iteration is terminated, the contact area and the corresponding node contact force are obtained, and the node contact force is divided by the nominal area of the element to calculate the contact pressure and the distribution characteristics.

8. The method of claim 7, wherein, In step 5, the real contact area is calculated through the size of the interface gap, and the expression of the real contact area of the rough surface is: ， A e , A ep1 , A ep2 and A p represent the real contact area in the elastic phase, the first elastoplastic phase, the second elastoplastic phase and the plastic flow phase, respectively.

9. An apparatus comprising a memory, a processor, and a computer program stored on the memory and executable on the processor, wherein, The processor implements the steps of the method of any one of claims 1 to 8 when executing the computer program.

10. A computer-readable storage medium having stored thereon a computer program, characterized in that, The computer program is executed by the processor to implement the steps of the method of claims 1 to 8.

Images