A computer-implemented contact force evaluation method for evaluating a contact force against a supporting member of a tube bundle, required for preventing the tube bundle from causing a self-excited vibration due to an excitation force (F ex) of a fluid, the tube bundle being a heat-transfer tube bundle of a steam generator which is positioned in the fluid and supported by the supporting member, the method comprising: a contact force setting step of setting the contact force, based on a value received from an input unit; a probability density function calculation step of calculating a probability density function (Pdf(F rv)) of a reaction force (F rv) received by the supporting member from the tube bundle in response to the excitation force (F ex) which is input to the tube bundle, using a vibration analysis model of the tube bundle and the supporting member; a probability calculation step of calculating a probability that a reaction force (F rv) equal to or higher than the set contact force occurs, based on the calculated probability density function (Pdf(F rv)); an evaluation step of evaluating the set contact force, based on the calculated probability, and based on the results of the evaluation step an outputting step of outputting the contact force to be applied between the tube bundle and the supporting member; wherein the evaluation step includes: estimating a wear amount which occurs between the tube bundle and the supporting member during a predetermined period, based on the calculated probability; and evaluating the set contact force, based on the estimated wear amount, and wherein the probability density function (Pdf(F rv)) is calculated by synthesizing a first probability density function (Pdf(F slide)) and a second probability density function (Pdf(F lift)), the first probability density function (Pdf(F slide)) corresponding to a first reaction force (F slide) received by the supporting member in a tangential direction from a contact surface between the supporting member and the tube bundle, the second probability density function (Pdf(F lift)) corresponding to a second reaction force (F lift) received by the supporting member in a normal direction from the contact surface. a contact force setting step of setting the contact force, based on a value received from an input unit; a probability density function calculation step of calculating a probability density function (Pdf(F rv)) of a reaction force (F rv) received by the supporting member from the tube bundle in response to the excitation force (F ex) which is input to the tube bundle, using a vibration analysis model of the tube bundle and the supporting member; a probability calculation step of calculating a probability that a reaction force (F rv) equal to or higher than the set contact force occurs, based on the calculated probability density function (Pdf(F rv)); an evaluation step of evaluating the set contact force, based on the calculated probability, and based on the results of the evaluation step an outputting step of outputting the contact force to be applied between the tube bundle and the supporting member; wherein the evaluation step includes: estimating a wear amount which occurs between the tube bundle and the supporting member during a predetermined period, based on the calculated probability; and evaluating the set contact force, based on the estimated wear amount, and wherein the probability density function (Pdf(F rv)) is calculated by synthesizing a first probability density function (Pdf(F slide)) and a second probability density function (Pdf(F lift)), the first probability density function (Pdf(F slide)) corresponding to a first reaction force (F slide) received by the supporting member in a tangential direction from a contact surface between the supporting member and the tube bundle, the second probability density function (Pdf(F lift)) corresponding to a second reaction force (F lift) received by the supporting member in a normal direction from the contact surface.