JP2009098030A - Simulation method and program for heat treatment - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To speedily and precisely compute a temperature distribution of a work piece in a cooling operation after a heating treatment, without the need for a complex boiling analysis of a cooling liquid. <P>SOLUTION: In a heat transfer simulation section 2, a heat transfer at an interface of the work piece (solid) and a fluid (cooling oil) is obtained by an interface heat transfer computing section 5 from a simplified computation, using "equivalent heat transfer coefficient", and a heat transfer in the work piece and a heat transfer in the fluid, taking place along with the heat transfer at the interface are computed, respectively, by a solid heat transfer computing section 4 and a fluid heat transfer computing section 6, whereby temperature distributions of all gears are obtained as a function of time. Next, in a structure analysis simulation section 3, each gear is set as an object, and its stress and deformations are computed by using a temperature expressed as a function of time which is obtained by the heat transfer simulation section 2. <P>COPYRIGHT: (C)2009,JPO&INPIT

Description

  The present invention relates to a heat treatment simulation method and a simulation program for simulating the temperature distribution of a workpiece during cooling after heat treatment.

  A heat treatment such as quenching of steel is widely applied to parts of various products because the strength can be increased without increasing the weight of the material. For example, in vehicles such as automobiles, weight reduction of component parts is promoted for the purpose of environmental improvement and energy saving, and parts requiring strength such as gears used in transmissions are quenched by heat treatment. Thus, it is possible to reduce the thickness and weight while maintaining the strength.

  However, quenching improves strength, but exists as a phenomenon in which distortion (deformation or deformation) cannot be avoided. In vehicles such as the automobiles mentioned above, the quietness of the car is required as part of improving comfort when riding, and high precision of the parts dimensions is required, so quenching distortion is a concern, There is a need for optimization of the heat treatment process.

  Conventionally, trial and error methods based on past experience have been used for heat treatment optimization, but in recent years, techniques for simulating and optimizing heat treatment have been attempted. Various proposals have been made.

For example, Patent Document 1 discloses a technique capable of easily and accurately estimating a quenching range by simulation using a finite element method, and Patent Document 2 discloses a heat treatment using a finite element model. Thus, a technique capable of analyzing distortion with high accuracy is disclosed.
JP 2001-49333 A JP 2003-194754 A

  However, the techniques disclosed in Patent Documents 1 and 2 do not consider boiling generated at the interface between the workpiece (solid) and the coolant (fluid) with respect to the heat transfer when the workpiece is cooled. It is difficult to accurately grasp the temperature distribution.

  Such heat transfer accompanied by boiling between a solid and a fluid can be simulated by numerically solving an unsteady heat conduction equation, but requires extremely complicated analysis. For this reason, a large-scale and high-speed computer device is required, which not only increases the cost required for the simulation, but is difficult to apply when heating and cooling a large number of workpieces at a time in an actual production process.

  The present invention has been made in view of the above circumstances, and a heat treatment simulation method capable of calculating the temperature distribution of a workpiece during cooling after heat treatment at high speed and with high accuracy without requiring complicated boiling analysis of the coolant. And it aims at providing a simulation program.

  In order to achieve the above object, a heat treatment simulation method according to the present invention is a heat treatment simulation method for simulating a temperature change in a coolant of a heat-treated workpiece, and an interface involving boiling between the workpiece and the coolant. The equivalent heat transfer coefficient for each temperature calculated in advance based on the cooling characteristics of the coolant is set, and the heat transfer from the workpiece to the coolant is calculated for each analysis element based on the equivalent heat transfer coefficient. And a step of calculating a heat distribution for each analysis element in a predetermined minute time step and calculating a temperature distribution of the work surface on a time unit basis.

The heat treatment simulation program according to the present invention is a heat treatment simulation program that can be executed by a computer for simulating a temperature change in a coolant of a heat-treated workpiece, and at the interface with boiling between the workpiece and the coolant, Setting an equivalent heat transfer coefficient for each temperature calculated in advance based on the cooling characteristics of the coolant, and calculating heat transfer from the workpiece to the coolant for each analysis element based on the equivalent heat transfer coefficient; ,
Calculating the heat transfer for each analysis element in predetermined minute time increments and calculating the temperature distribution of the work surface on a time unit basis.

  According to the present invention, the temperature distribution of the workpiece during cooling after the heat treatment can be calculated at high speed and with high accuracy without requiring complicated boiling analysis of the coolant, and a large number of workpieces can be obtained in an actual production process. Applicable when heating and cooling at a time, it is possible to predict the distortion generated in the workpiece with high accuracy.

  Embodiments of the present invention will be described below with reference to the drawings. 1 to 10 relate to an embodiment of the present invention, FIG. 1 is an explanatory view showing a cooling tank model, FIG. 2 is an explanatory view showing a temperature change of a workpiece due to thermal boiling of a fluid, and FIG. 3 is a heat transfer model. FIG. 4 is a functional block diagram of the simulation apparatus, FIG. 5 is an explanatory diagram showing the concept of heat transfer coefficient measurement, FIG. 6 is an explanatory diagram showing the temperature change of the metal rod in the cooling oil, and FIG. FIG. 8 is a flowchart of a temperature distribution simulation program, FIG. 9 is a flowchart of an interfacial heat transfer calculation program, and FIG. 10 is an explanatory diagram showing a gear distortion analysis result.

  The present invention simulates the temperature distribution of the workpiece surface in the cooling process after the heating process, taking into account the boiling of the cooling oil, with high accuracy for the distortion caused by quenching the workpiece such as gears and bearing races. Is predictable.

  Quenching is a heat treatment in which the steel is heated to the austenite region and then carburized and then quenched in a cooling medium such as oil to harden it as a martensite structure. The strain generated during quenching is the sum of transformation strain and thermal strain. Transformation strain is a change caused by volume expansion accompanying martensitic transformation, and thermal strain is caused by restraint between materials inside the object due to thermal expansion and contraction caused by temperature difference in the cooling process. Due to the restraint and the generation of thermal stress, deformation due to this thermal stress such as shape change or volume change of the object.

  In general, the thickness of the martensite layer produced by quenching is extremely thin (for example, about 0.3 mm), so that the size change due to expansion is small, and the main cause of quenching distortion is uneven cooling, that is, This is thought to be due to uneven cooling. The deformation at the time of cooling is thought to be caused by the difference in cooling rate between the thick part and thin part, the difference between the cooling rate between the surface and the inside, and the cooling rate difference depending on the position. It is necessary to accurately grasp the temperature distribution of the workpiece surface at the time.

  Taking the quenching of gears as an example, in the actual quenching process, dozens of gears are heat-treated at the same time, and the heat-treated gears are immersed in cooling oil and rapidly cooled, so that the surface texture of the material is changed. Harden. FIG. 1 shows an example in which a cooling tank is modeled. In a cooling oil 51 in the cooling tank 50, a substantially U-shaped duct-shaped cooling device 52 is provided.

  The plurality of heat-treated gears 100 are immersed in the cooling oil 51 in a state where they are put in the flange 101, and are arranged at the opening end of the main body duct 52a of the cooling device 52. An agitator 53 for agitating the cooling oil 51 is attached to the other opening end side of the main body duct 52a. The agitator 53 forces the cooling oil 51 from the inside of the main body duct 52a toward the gear 100. Sent. Further, a rectifying plate 54 for adjusting the flow velocity at the outlet of the cooling oil 51 to make the cooling portion uniform is disposed in the main body duct 52a.

  FIG. 2 shows an example in which the gear 100 and the cooling oil 51 are modeled by analysis elements (three-dimensional mesh), and the temperature change of the gear 100 in the cooling oil 51 is simulated. In this simulation, as shown in FIG. 2A, a workpiece (gear) 100 that is initially heat-treated at a predetermined temperature and has a uniform temperature state on the surface and inside of the member as shown in FIG. Then, it is immersed in the cooling oil 51, and immediately after the immersion is set to time t = 0. The temperatures of the gear 100 and the surrounding cooling oil 51 when the times t1, t2, t3 (t1 <t2 <t3) have elapsed are calculated using the heat transfer coefficient between the gear 100 and the cooling oil 51. The temperature distribution is as shown in FIGS. In FIG. 2, the temperature is indicated by the concentration expression of the gear 100, and the concentration decreases as the temperature decreases.

  At this time, the temperature of the gear 100 rapidly decreases due to heat transfer to the cooling oil 51, but boiling occurs at the interface between the gear 100 and the cooling oil 51, so that complicated analysis is required. Therefore, if the heat transfer coefficient at the interface between the workpiece (metal) and the cooling oil (liquid) is calculated without considering boiling, the actual temperature distribution will be different, and the calculated workpiece temperature will be higher than the actual workpiece temperature. Get higher.

  This is because the cooling oil boils in the middle and a boundary layer is generated between the solid and the liquid, and there is a temperature region that is cooled by latent heat. As shown in FIG. 3, the heat transfer phenomenon in the cooling process is divided into three zones: a solid layer (gear) Rs, a boundary layer Rc, and a fluid layer (cooling oil) Rf. The cooling at the time is divided into three stages: immediately after quenching, boiling stage and convection stage. Immediately after quenching, the cooling rate is the slowest in the vapor film stage, and the cooling rate is the fastest in the boiling stage. When the surface temperature of the workpiece reaches a predetermined temperature (for example, about 400 ° C.), the convection stage is entered and the cooling rate becomes slow again.

  The boiling stage in the boundary zone between the solid layer (gear) and the fluid layer is a nucleate boiling stage where the generated bubbles do not disappear even if they leave the heat transfer surface, and the bubbles cover the heat transfer surface and form a film In the nucleate boiling stage, the heat transfer is faster due to evaporation, so the thermal conductivity is large, and in the film boiling stage, the heat transfer is slowed by the adiabatic action of the air film. The degree becomes smaller. In heat transfer involving boiling, the transition between nucleate boiling and film boiling becomes a problem.

  Such a problem can be simulated by thermal fluid analysis by computational fluid dynamics (CFD). In general, the behavior of a thermal fluid is governed by the conservation laws of mass, momentum, and energy, and the Navier-Stokes formula, mass conservation formula, and energy conservation formula are basic equations. In CFD, the Navier-Stokes equation is solved progressively using the Finite Volume Method (FVM), etc., and the force and moment obtained by pressure integration are added to the equation of motion. This is reflected in the handling of the analysis element. Thereby, the flow velocity and temperature of the object can be understood, and the behavior (temperature and flow velocity) of the cooling oil and the object in the quenching process can be reproduced.

  However, in order to handle the boiling phenomenon of fluid, a calculation in consideration of bubbles generated by boiling is necessary. For this reason, simply determining the solution of the unsteady heat conduction equation using CFD not only requires a large amount of calculation and calculation time, but also increases the cost required for the simulation. I have to say otherwise.

  For this reason, as shown in FIG. 4, the simulation apparatus 1 of this embodiment includes a heat transfer simulation unit 2 that simply approximates the boiling of the boundary layer and calculates heat transfer between the workpiece and the cooling oil, A structural analysis simulation unit 3 that performs a structural analysis using the result of the thermal fluid analysis in the movement simulation unit 2 as a boundary condition, and the heat transfer simulation unit 2 and the structural analysis simulation unit 3 are sequentially coupled to deform the gear 100. I try to predict.

  Specifically, the simulation apparatus 1 is configured using a single computer such as a microcomputer or a personal computer, or a plurality of computers connected to each other via a network. Below, the example which comprises the simulation apparatus 1 with a single computer for convenience is demonstrated.

  The simulation apparatus 1 includes an arithmetic device 10 that forms a heat transfer simulation unit 2 and a structural analysis simulation unit 3, an input device 11 such as a keyboard and a mouse, a display device 12 such as a CRT and a liquid crystal display, and an external storage such as a magnetic disk and an optical disk. The apparatus 13 etc. are provided.

  The arithmetic device 10 includes an internal memory such as a CPU, a ROM and a RAM, an input / output interface, and the like, and includes a simulation program stored in an internal ROM, an external storage device 13 and an external storage medium, or a network (not shown) The temperature distribution of the workpiece surface when the simulation program loaded from the outside via the communication device is executed by the CPU and the workpiece (object) to be analyzed instructed via the input device 11 is immersed in the cooling bath 50 And the stress and deformation of each part of the workpiece are calculated based on this temperature distribution.

  That is, the arithmetic unit 10 first calculates the flow rate of the cooling oil 51 and the accompanying heat transfer between the gear 100 and the cooling oil 51 in the heat transfer simulation unit 2, and sets all the gears 100 as a function of time. Obtain the temperature distribution of. Next, the structural analysis simulation unit 3 calculates the stress and the deformation amount for each gear 100 using the temperature as a function of the time obtained by the heat transfer simulation unit 2.

  Specifically, the heat transfer simulation unit 2 includes a solid heat transfer calculation unit 4, an interface heat transfer calculation unit 5, and a fluid heat transfer calculation unit 6, and provides an interface between a workpiece (solid) and a fluid (cooling oil). The heat transfer is simply calculated by the interfacial heat transfer calculation unit 5, and the heat transfer in the work and the heat transfer in the fluid accompanying the heat transfer at the interface are respectively calculated by the solid heat transfer calculation unit 4 and the fluid heat transfer calculation unit 6. calculate.

  The calculation in the solid heat transfer calculation unit 4 and the fluid heat transfer calculation unit 6 is a calculation based on a general heat conduction equation, while the calculation in the interface heat transfer calculation unit 5 is performed between the solid layer and the fluid layer. Introducing the concept of "equivalent heat transfer coefficient" to the boiling analysis in the boundary layer of, reducing the amount of calculation and speeding up.

The equivalent heat transfer coefficient is a coefficient for simply analyzing the heat transfer that changes due to the boiling phenomenon of the boundary layer, and is experimentally obtained using the cooling oil that is actually used. That is, boiling heat transfer is a heat transfer accompanied by a phase change, but Newton's cooling law can be applied under conditions where the temperature difference is small. The amount of heat (heat flux) Q at this time is expressed by the following equation (1), where As is the surface area of the solid (wall surface of the boundary layer; workpiece), Ts is the solid temperature, and Tf is the temperature of the cooling medium (cooling oil). It can be described by.
Q = −α · As · (Ts−Tf) (1)
Where α is the thermal conductivity of the solid surface

On the other hand, the amount of heat Q transferred from the solid can be expressed by the heat capacity and cooling rate of the solid. If the specific heat of the solid in the cooling medium is Cs, the density is ρs, and the volume is V, the amount of heat Q is ( 2) It can be expressed by an equation.
Q = −Cs · Ts · ρs · Vs · dTs / dt (2)

The thermal conductivity α in equation (1) is a state quantity that varies depending on the shape of the wall surface, the type of fluid, the state of flow, etc. In this embodiment, the thermal conductivity α in equation (1) is usually Equation (2) is applied as the heat transfer coefficient (equivalent heat transfer coefficient) hm for each temperature in the fluid layer. As a result, the heat transfer coefficient hm can be expressed by the following equation (3).
hm = −Q / (As · (Ts−Tf))
= (Cs · Ts · ρs · Vs / (As · (Ts−Tf))) · dTs / dt (3)

Here, since the specific heat Cs, the density ρs, the volume V, and the surface area As of the solid to be analyzed are known, when these are collectively replaced with a constant K, the expression (3) is expressed by the following expression (4): Can be expressed as
hm = (K · Ts / (Ts−Tf)) · dTs / dt (4)

  Accordingly, the heat transfer coefficient hm expressed by the equation (4) is calculated for each temperature based on the data measured in advance, and the solid (work) and the cooling medium (cooling oil) through the boundary layer with boiling are used. When calculating the heat transfer, the temperature distribution can be obtained accurately by analyzing the interface with the solid by varying the heat transfer coefficient for each temperature.

  The heat transfer coefficient hm is calculated based on the concept of data measurement schematically shown in FIG. That is, the metal bar 200 is used as a member corresponding to the workpiece to be analyzed, and the metal bar 200 is heated to the heating temperature in the actual quenching process, and then inserted into the cooling oil 201 that is actually used. The measuring device 202 measures the temperature of the metal bar 200 for each hour. And the cooling performance of the cooling oil 201 is investigated by this measurement, and the heat transfer coefficient hm is obtained.

  As the metal rod 200, it is desirable to use a rod having a diameter as small as possible, and the material is the same as that of the object to be analyzed, for example, when the material of the gear to be analyzed is chrome molybdenum steel (SCM). However, the heat transfer of the cooling oil can be efficiently measured by using a silver (Ag) rod having a high thermal conductivity. In addition, the calculation accuracy can be improved by making the surface roughness of the metal bar 200 the same as that of the analysis target as much as possible.

  For example, when the quenching temperature of the gear to be analyzed is 850 ° C., the metal rod 200 heated to 850 ° C. is immersed in the cooling oil 201 held at the cooling oil temperature 150 ° C. in the actual cooling process. The temperature of the metal rod 200 is measured in predetermined minute time increments (the time increment is preferably as short as possible), with the time immediately after immersion as 0.

  By this measurement, as shown in FIG. 6, the relationship between the temperature of the metal bar 200 and time can be obtained. This relationship is obtained by discrete measurement of the solution curve for the temperature Ts in the equation (4). The heat transfer coefficient hm is obtained from the temperature difference for each minute time step, and the minimum value of the heat transfer coefficient hm is obtained. As shown in FIG. 7, the amount of heat (boiling heat) Qf due to boiling of the interface can be obtained for each temperature.

  Thus, by adding the boiling heat Qf to the heat transfer amount Qk of the interface calculated by the heat transfer between the normal solid and the fluid, the total heat of the interface at the time of boiling can be obtained without requiring a complicated boiling analysis. The movement amount Qt can be obtained quickly and accurately.

  Next, simulation processing for obtaining the temperature distribution of the work immersed in the cooling oil will be described with reference to the flowchart of FIG.

  In this process, first, in the first step S1, the temperature of the solid part is defined. The temperature of the solid part is defined as the temperature of the solid layer composed of analysis elements that model the workpiece (gear) to be analyzed. The analysis element is composed of, for example, a three-dimensional mesh by a finite volume method, and the temperature of all the meshes of the solid layer is initially set to the heating temperature of quenching (for example, 850 ° C.).

  Next, in step S2, the temperature of the coolant portion is defined. The temperature of the coolant section is defined as the temperature of the fluid layer composed of the analysis element (three-dimensional mesh) that models the cooling oil into which the workpiece is immersed. Is initially set to a temperature (for example, 150 ° C.).

  In the subsequent step S3, meshes that form a solid surface that comes into contact with the liquid are grouped into a group G. Then, in step S4, the simulation time TM is reset to 0, and in step S5, the simulation time TM is incremented by a minute time step Δt. In step S6_SOL, the heat conduction inside the solid (work) is calculated. The process of step S6_CL for calculating the heat transfer at the interface (interface between the workpiece and the cooling oil) and the process of step S6_LIQ for calculating the heat conduction in the liquid (cooling oil) are switched and executed.

  The interface heat transfer calculation in step S6_CL is calculated for each mesh using the boiling heat Qf calculated from the heat transfer coefficient hm described above. Specifically, it is calculated by an interfacial heat transfer calculation program shown in the flowchart of FIG. The processing by this interface heat transfer calculation program will be described later.

  On the other hand, the heat transfer calculation inside the solid in step S6_SOL is based on the heat transfer amount of the interface calculated in step S6_CL, the heat amount for each mesh from the solid to the surface, and the boundary condition due to the heat transfer between adjacent meshes is unsteady. Apply to the heat equation and calculate. Similarly, the heat transfer calculation inside the liquid in step S6_LIQ takes into account the influence of the flow simulated for each mesh that diffuses to the surroundings using the cooling tank model (see FIG. 1) based on the heat transfer amount at the interface. To calculate unsteady.

  After performing the heat transfer calculation inside the solid, interface, and liquid as described above, the process proceeds to step S7, where the temperature of the solid part and the temperature of the liquid part are calculated from the amount of heat transferred for each mesh, and the temperature before the heat transfer calculation. Is updated for each mesh. In step S8, it is determined whether or not the simulation time TM has reached a preset end time. If the simulation time TM has not reached the end time, the process returns to step S5 and the simulation time TM is incremented by a minute time step Δt. The above processing is repeated until the time TM reaches the end time.

  Next, the interface heat transfer calculation in step S6_CL will be described with reference to the flowchart of FIG.

  In the interfacial heat transfer calculation program of FIG. 9, first, the mesh number i of the mesh G_Noi belonging to the group G is reset to 0 in step S11. In step S12, the mesh number i is incremented, and in step S13, heat is transferred from the solid (work) to the liquid (cooling oil) based on the temperature difference between the mesh G_Noi and the mesh on the liquid side adjacent to the mesh G_Noi. The amount of heat transferred Qk is calculated.

  Next, it progresses to step S14 and calculates temperature Ti of the mesh G_Noi after moving calorie | heat amount Qk. The temperature Ti in step S14 is based on the updated heat quantity Qt described below, and the first time is the same as the temperature in step S13.

  Then, in step S15, the boiling heat Qf is obtained based on the temperature Ti of the mesh G_Noi, and in step S16, the boiling heat Qf based on the equivalent heat transfer coefficient hm is added to the previously calculated moving heat amount Qk of the mesh G_Noi. Then, it is determined and updated as the amount of moving heat taking into account the boiling of the mesh G_Noi (Qt = Qk + Qf).

  After determining the amount of heat Qt of movement of the mesh G_Noi, the process proceeds to step S17 to check whether or not the calculation for all the meshes in the group G has been completed. If calculation for all meshes in group G has not been completed, the process returns from step S17 to step S12 to increment mesh number i, and the above processing is repeated to complete calculation for all meshes in group G. When you do, exit this program.

  When the temperature distribution on the workpiece surface is obtained for each time by the above heat transfer calculation, the structural analysis simulation unit 3 predicts the strain generated in the workpiece from this temperature distribution. For example, using a temperature distribution for each hour when the gear 100 is immersed in the cooling oil 51, a load is set by thermal stress deformation, and a tooth taper / taper amount of the tooth trace is calculated. In this case, by comparing the time-dependent temperature distribution of the target gear material with the continuous cooling transformation curve (CCT diagram), the martensite distribution is estimated, and the transformation stress due to the volume change of this martensite portion is taken into account, Calculation accuracy can be improved.

  The prediction result of this distortion is visualized by post processing and displayed on the display device 12. For example, as shown in FIG. 10, the gear 100 is three-dimensionally displayed on the display screen 250, and the degree of distortion of each part is identified and displayed by hue or the like. The relationship between the hue and the degree of distortion is indicated by a bar-shaped graph 251 at the bottom of the gear 100.

  As described above, in the present embodiment, the temperature change in the coolant accompanying the boiling of the heat-treated workpiece can be grasped quickly and accurately without requiring a complicated boiling analysis. This not only reduces the cost required for simulation by eliminating the need for a large-scale and high-speed computer device, but also applies distortion when heating and cooling a large number of workpieces at the same time in an actual production process. Can be predicted with high accuracy.

Explanatory drawing showing the cooling tank model Explanatory drawing showing temperature change of workpiece due to thermal boiling of fluid Illustration of heat transfer model Functional block diagram of the simulation device Explanatory diagram showing the concept of heat transfer coefficient measurement Explanatory drawing showing temperature change of metal rod in cooling oil Explanatory diagram showing the relationship between the temperature at the interface and the amount of heat transferred Flow chart of temperature distribution simulation program Flow chart of interface heat transfer calculation program Explanatory drawing showing the result of gear distortion analysis

Explanation of symbols

DESCRIPTION OF SYMBOLS 1 Simulation apparatus 2 Heat transfer simulation part 3 Structural analysis simulation part 4 Solid heat transfer calculation part 5 Interface heat transfer calculation part 6 Fluid heat transfer calculation part Q Calorie | heat amount Qf Boiling heat hm Equivalent heat transfer coefficient

Claims (4)

A heat treatment simulation method for simulating a temperature change in a coolant of a heat-treated workpiece,
An equivalent heat transfer coefficient for each temperature calculated in advance based on the cooling characteristics of the coolant is set at the interface with boiling between the work and the coolant, and from the work based on the equivalent heat transfer coefficient Calculating heat transfer to the coolant for each analytical element;
Calculating the heat transfer for each analysis element in predetermined minute time increments and calculating the temperature distribution of the workpiece surface in time units.   The heat treatment simulation method according to claim 1, wherein the equivalent heat transfer coefficient is calculated in advance by applying Newton's cooling law for each temperature of the workpiece surface.   3. The heat treatment simulation according to claim 2, wherein a temperature change of the member corresponding to the workpiece in the coolant is measured every predetermined time, and the equivalent heat transfer coefficient is calculated from a measured value of the temperature change. Method. A computer-executable heat treatment simulation program for simulating a temperature change in a coolant of a heat-treated workpiece,
An equivalent heat transfer coefficient for each temperature calculated in advance based on the cooling characteristics of the coolant is set at the interface with boiling between the work and the coolant, and from the work based on the equivalent heat transfer coefficient Calculating heat transfer to the coolant for each analytical element;
Calculating a heat transfer for each analysis element at predetermined minute time intervals, and calculating a temperature distribution of the work surface on a time unit basis.