To provide a calculation method for contact surface pressure allowing the high-accuracy calculation of the contact surface pressure as compared to a uniform surface pressure element, allowing easy mesh division as compared to a triangular element, having a short calculation time as compared to a finite element method or a boundary element method, and capable of being easily used.

A square area R slightly larger than a contact area is set, and the area is cut into a plurality of square elements e in a mesh form arranged in two side directions orthogonal to each other. The arbitrary square element e is subdivided. It is assumed that a linear contact surface pressure distribution according to a shape function is generated in each the square element e. An equality wherein a value obtained by superimposing an influence exerted on displacement in the contact surface normal direction of each node by the square element e about all the square elements is equal to the displacement in the contact surface normal direction of each the node is obtained. An equality wherein the sum of loads supported by the square elements e is equal to an external force applied between contact objects is obtained. The contact surface pressure of each square element satisfying both the equalities is obtained.