ISIJ International
Online ISSN : 1347-5460
Print ISSN : 0915-1559
ISSN-L : 0915-1559
Forming Processing and Thermomechanical Treatment
Multi-step Simulation of Vacuum-carburizing Reactor Based on Kinetics Approach
Soichiro Makino Masahide InagakiHideaki IkehataKouji TanakaHiroyuki InoueKoji Inagaki
Author information
JOURNAL OPEN ACCESS FULL-TEXT HTML

2020 Volume 60 Issue 9 Pages 2000-2006

Details
Abstract

A novel and efficient simulation technique for the purpose of optimization of vacuum-carburizing process was proposed. This method consists of three steps: calculation of gas convection and diffusion, calculation of only gas diffusion, and calculation of carbon diffusion in steel. The first step provides the gas convection velocity that is employed in the second step. Adsorption rate of carbon on the steel surface is obtained in the second step, and carbon concentration in the steel is calculated in the third step based on the adsorption rate of carbon.

Experiments were conducted to verify the proposed method in both laboratory- and industrial-scale reactors. Comparison of the computational predictions to the experimental data revealed that the proposed simulation technique enabled accurate prediction of the adsorption rate of carbon on the steel surface at various temperature conditions, the amount of carburized carbon at each operating time, and the profile of carbon concentration in the steel that is, in other words, the carburized depth. In addition, the calculation of the industrial-scale reactor, whose simulation model consisted of approximately seven million computational meshes, was completed within about two days. Therefore, the proposed simulation technique could be used to control and optimize the process in industrial vacuum-carburizing reactors.

1. Introduction

Vacuum carburizing, which involves direct carburization from a source gas under reduced pressure, can increase the carbon concentration on a steel surface compared with that of traditional gas carburization based on the Boudouard equilibrium reaction, and can therefore reduce the processing cost and energy requirements.1,2,3) For the case of vacuum carburizing, the carbon adsorption rate from the source gas onto the steel surface is so high that the distributions of the source gas and by-product gas concentrations in the reactor, which can affect the carburization, cannot be neglected. It is therefore difficult to control carburizing process by using conventional analysis based on the equilibrium partial pressure for gas carburization. Thus, developing a new method of prediction and analysis that considers the distribution of gas concentration in the reactor is essential.

Khan et al.4) investigated the thermal decomposition behaviour of propane gas in vacuum carburization based on experiment and reaction kinetics theory, but in their study they considered only the thermal decomposition of source gas and they considered neither the carbon adsorption on a steel surface from the source gas nor the carbon diffusion in the steel, which are important factors in vacuum carburization. Yada et al.5) proposed a pioneering simulation method that calculates the diffusion of carbon into steel while calculating the gas flow and gas diffusion in vacuum reactors. With that method, the carburization can be predicted while considering the source gas concentration distribution in the reactor, and it therefore appears promising as a new method of carburizing simulation. In their report, however, they addressed only carburization of one thin plate as the domain of the analysis. Hence, if an application to industrial vacuum-carburizing reactors, which are actual large-scale multiple-workpiece mass production processes, is considered, the following problems will be occurred from the perspectives of simulation scale and costs.

In many processes, the distribution of carbon concentration in the steel forms mainly on the order of 1 mm from the surface. In particular, it involves a sharp concentration gradient on the order of several hundred μm from the surface, and therefore in the surface vicinity a minimum computational grid spacing of 10 μm or less is desirable (in the previous study,5) the minimum spacing was set to approximately 10 μm or less). This leads that a minimum grid spacing for the calculation of gas phase has to be small as well, and also the computational time steps have to be extremely small values from the viewpoint of simulation stability in terms of the Courant number and the diffusion number.6) Generally, the number of computational meshes for a carburizing reactor is huge because the reactor scale is generally several-score centimeters to several meters. Moreover, if workpieces have three-dimensional (3D) complex shape, the number of meshes becomes huger. Accordingly, for the case of workpieces with a complex shape in an actual 3D carburization reactor, it seems to be difficult to conduct the simulation of the carburization processes by using the method of the previous study in the perspectives of both computational scale and time.

In the present study, we therefore propose a simulation method that enables prediction of gas flow, gas concentration distributions, and carbon concentration distribution in steel in a practical computational time. We also verify the proposed method with both a single-workpiece test reactor and a multi-workpiece industrial reactor.

2. Proposal of Multistep Computational Method

As described in the previous section, because there are several length and time scales in a vacuum-carburizing reactor, such as gas convection, gas diffusion and carbon diffusion in steel, it can be assumed that strong-coupling calculation of these phenomena will result in extremely poor computational efficiency. In order to solve this problem, as shown in Fig. 1, we focus on the differences in time scale in the carburizing reactor and propose a multistep method that separately computes gas convection and gas diffusion in Step 1, gas diffusion in Step 2, and carbon diffusion in steel in Step 3.

Fig. 1.

Schematic overview of the multi-step simulation for a vacuum-carburizing reactor.

Let us first consider Steps 1 and 2. In a vacuum-carburizing reactor, the drastic gas expansion occurs and the gas velocity Uin becomes very large at a gas inlet because the source gas is supplied to the vacuum reactor from atmospheric outside. With the inside of the reactor in a vacuum state, the transport of the source gas occurs mainly as gas diffusion rather than as gas convection. For the gas convection time scale τ c g can be estimated as Lg/Uin, where Lg is the mean distance between the inlet section and the workpiece. The gas diffusion time scale τ d g can be estimated as Lg2/(Dgπ2), where Dg is the gas diffusion coefficient. Thus, the time scales ratio Γ is then   

Γ= τ d g τ c g = L g2 D g π 2 L g U in = L g U in D g π 2 , (1)
Accordingly, if Γ >> 1, the time for the gas convection to reach a steady state becomes sufficiently small compared with the time for gas diffusion to reach a steady state. Hence, one can calculate only the gas diffusion after obtaining the steady state of gas convection under Γ >> 1. For the case of a carburizing reactor that inputs source gas with a small nozzle, Γ can be estimated as Γ > 101 because the inlet velocity Uin becomes a large value. On the other hand, for the case of a reactor incorporating baffle plates, Γ can be estimated as Γ = 10−1 − 100 because Uin becomes comparatively lower. Accordingly, it is necessary to ascertain the value of Γ in advance for the isolation of Steps 1 and 2.

Next, in regard to Steps 2 and 3, we obtain the ratio between the gas diffusion time scale τ d g and the in-steel carbon diffusion time scale τ d s as   

Λ= τ d s τ d g = L s2 D s π 2 L g2 D g π 2 = D g D s ( L s L g ) 2 , (2)
where Dg is the coefficient of carbon diffusion in steel and Ls is the distance of carbon diffusion in steel (for example, the carburization depth). If Λ >> 1, the time until gas diffusion reaches a steady state becomes sufficiently small compared with the time until carbon diffusion in steel reaches a steady state so that after the solution for gas diffusion is obtained the in-steel carbon diffusion can then be calculated using that result. Here, for example, assuming the general condition for vacuum carburizing, where the temperature is 1100–1300 K and the pressure is 3000 Pa or less, we then have Λ = 102 − 104. There is thus a large difference in time scale between Step 2 and Step 3, and then approximately we can calculate Step 2 and Step 3 separately.

3. Simulation Method

As discussed above, physical phenomena in a vacuum-carburizing reactor can be simulated by the multistep method that separately computes Step 1, Step 2, and Step 3, based on the estimated values of Γ and Λ. Here, numerical procedures in each step will be introduced.

3.1. Step 1

3.1.1. Governing Equations

The governing equations for gas flow consist of mass and momentum conservation equations for the gas mixture, an energy conservation equation, and mass conservation equations for the chemical species as follows:   

ρ t + (ρ u j ) x j = k ( w s ) k , (3)
  
(ρ u i ) t + (ρ u i u j ) x j =- p x i + x j ( η( u i x j + u j x i ) ) , (4)
  
(ρ Y m ) t + (ρ u j Y m ) x j = x j ( ρ D m Y m x j ) + ( w s ) m , (5)
  
P 0 = ρRT M , (6)
where ρ is the density of the gas mixture, t is the time, xj is the coordinate direction (x1 = x, x2 = y, x3 = z), uj is the velocity vector, ws is the source term related to the surface reaction rate which will be discussed later, p is the dynamic component of pressure, η is the viscosity of the gas mixture, Ym is the mass fraction of each species; the subscript m indicates the material number of each species, Dm is the molecular diffusion coefficient of each species, P0 is the thermodynamic component of pressure, R is the ideal gas constant, and M is the molecular weight of the gas mixture. The summation convention is applied only to the direction, i, j, k, but not to the material number, m. In the present study, an C2H2−H2 gas mixture was examined. For the wall surface conditions, non-slip condition and the Neumann condition are adopted for gas flow and mass fraction, respectively. Note here that the mass flux of each species based on the surface reaction rate ws is calculated on the surface of workpieces.

3.1.2. Thermophysical Properties

The temperature dependences of the viscosity and the molecular diffusion coefficient are calculated based on the semi-theoretical assumption by Chapman and Enskog.7) The single component viscosities are given by   

η m =2.6693× 10 -6 M m T σ m 2 Ω m (2,2) , (7)
where σm is the Lennard-Jones collision diameter, and Ω m (2,2) is the collision integral, which is dependent on the reduced temperature defined by   
T * = k B ε m T. (8)
In the above expression, kB is the Boltzmann constant and εm is the parameter which indicates the Lennard-Jones potential well depth. A simple empirical formula is used to determine the gas mixture viscosity,   
η= X m η m , (9)
where Xm is the molar fraction.

The binary molecular diffusion coefficients are given by   

D mn =1.8583× 10 7 T 3/2 1/ M m +1/ M n p σ mn 2 Ω mn (1,1) . (10)
Here, σmn is the mean collision radius between specie m and specie n given by   
σ mn = 1 2 ( σ m + σ n ). (11)
The collision integral Ω i (1,1) is also dependent on the reduced temperature given in Eq. (12).   
T ** = k B ε mn T. (12)
The parameter εmn which indicates the Lennard-Jones well depth is calculated using the well depth parameters of specie m and specie n as below:   
ε mn = ε m ε n . (13)
For the effective diffusion coefficient, the following semi-empirical formula by Wilke8) and Bird9) is also used:   
D m = 1- X m mn ( X m / D mn ) . (14)

3.1.3. Surface Reaction Model between Source Gas and Steel Surface

Regarding the transport of carbon between the source gas and the steel surface, Morita et al.10,11) reported that the source gas is decomposed close to the steel surface, and consequently graphite is generated on the surface and it is absorbed into the steel immediately. In addition, Yada et al.5) conducted a numerical simulation of vacuum-carburizing reactor and proposed a chemical reaction model which consists of C2H2 adsorption on the steel. In accordance with these previous studies, here, the C2H2 adsorption is described as a first-order reaction with respect to the C2H2 concentration, as follows:   

C 2 H 2 2C * +H 2 , (15)
where C* indicates the adsorbed carbon.

In this model, the rate of C2H2 adsorption, ws, is given by:   

w s = k s C C2H2 , (16)
where ks is the reaction rate constant, CC2H2 is the molar concentration of C2H2. The reaction rate constant ks is given by,   
k s =Aexp(-E/ RT) , (17)
where A is a pre-exponential factor, and E denotes the activation energy.

3.2. Step 2

In the calculation of Step 2, only the gas diffusion is considered.

Using the steady-state solution of gas convection velocity, U j ¯ , obtained in Step 1, mass conservation equations for each species as follows are calculated,   

(ρ Y m ) t + (ρ U j ¯ Y m ) x j = x j ( ρ D m Y m, x j ) + ( w s ) m , (18)
  
P 0 = ρRT M . (19)
The surface reaction rate ws can be calculated by Eq. (16). The diffusion coefficient can also be obtained by Eq. (14), and the boundary conditions are the same as those of Step 1.

3.3. Step 3

In Step 3, the governing equation is a mass conservation equation of carbon in the steel (austenite phase) only in the depth direction (ξ direction) as follows:   

( ρ A Y C ) t = ξ ( ρ A D C Y C ξ ) . (20)
Here, ρA is the density of austenite, YC and DC are the mass fraction of carbon in austenite and the diffusion coefficient of carbon in austenite, respectively. DC can be calculated by the Wells’ equation:12)   
D C = D 0 exp( - Q R 0 T ) , (21)
where,   
D 0 =5.0398× 10 -5 exp(-1.4589 W C ), (22)
  
Q=-2   908.1 W C 2 -21   511 W C +154   490, (23)
  
W C = Y C ×100. (24)
Here, by using the solubility limit of carbon in steel YCsat and the carbon fraction of base material YC0, the boundary condition of YC is determined as follows:

at ξ = 0 (on the steel surface);   

if    Y C < Y Csat    ;    ρ A D C Y C ξ = w s , (25)
  
if    Y C Y Csat    ;    Y C = Y Csat . (26)
at ξ = ∞ (in the base material);   
Y C = Y C0 . (27)
In this study, YCsat and YC0 were set to 0.0135 (1.35 mass%) and 0.002 (0.2 mass%), respectively, in accordance with the previous studies.13,14)

3.4. Simulation Procedure

Figure 2 shows a summary of the simulation procedure. First, in Step 1, Eqs. (3), (4), (5), (6) are calculated using Eqs. (9), (14), and (16). Then, if Γ >> 1, after the spatial average value of gas velocity reaches a steady state, the value is transferred to Step 2 and Eqs. (18) and (19) are computed with Eqs. (14) and (16). This computation is performed until the mean value of the carbon adsorption rate on the surface of workpieces reaches a steady state. If Γ ≤ 1, then because gas convection and gas diffusion cannot be isolated, Step 1 computing must be performed until the mean carbon adsorption rate on the surface of workpieces reaches a steady state, without performing Step 2. The steady-state carbon adsorption rate thus obtained is transferred to Step 3, Eqs. (20) and (21) are solved with the boundary condition of Eq. (25) or (26), and the carbon concentration distribution in the direction of depth at a selected position on the workpiece surface and its time change can then be computed. Furthermore, since Λ >> 1 can be estimated in nearly all vacuum carburizing reactors as described above, it is assumed that we can calculate Steps 1 and 3 separately, or Steps 2 and 3 separately as well in this study. Hence, the condition where Λ >> 1 does not work is not applicable to this study.

Fig. 2.

Computational flow-chart of the multi-step simulation.

4. Experimental Method

Two experiments were performed to verify the simulation results; one using a single-workpiece test reactor (the “test reactor”) and the other using a multi-workpiece industrial reactor (the “industrial reactor”). Figure 3 shows the configuration of the test reactor. As shown, a cylindrical workpiece of 18 mm diameter and 50 mm height (case-hardening steel, SCr420) was suspended in a cylindrical pipe reactor of 55 mm inner diameter, and 100% C2H2 gas was input from below at 62 mL/min to perform the carburization. The temperature in the reactor Ttr was changed to 1223 K, 1273 K, 1323 K, and 1373 K, and the pressure in the reactor was 300 Pa. After the carburizing process, the C2H2 gas supply was stopped and the test piece was immediately dropped for rapid cooling into a nitrogen-gas cooling bath located in the lower part of the carburizing reactor. The four temperature conditions where the carburizing time tc was set to 5 s were investigated. In addition, the 11 conditions of carburizing time where Ttr is constant at 1373 K were investigated. In each investigation, the carburizing amount was assessed from the weight of the test piece before and after the carburization processing.

Fig. 3.

Illustration of the test reactor for vacuum carburizing. (Online version in color.)

The experiment with the industrial reactor, as shown in Fig. 4, was performed in a box-type vacuum-carburizing reactor with 28 gear-type workpieces (SCr420) each around 130 mm in outer diameter in a staggered array on a jig (the reactor in Fig. 4 is symmetrical about the xy section plane in the interior so that only half of one side is visible, and the jig is not shown). To analyse the carbon concentration distribution in steel, at the positions near the inlet section (indicated as A in the figure) and away from the inlet section (indicated as B in the figure), a cylindrical workpiece (SCr420) of 18 mm in diameter and 50 mm in height was tied on the jig by steel wire. C2H2 gas of 100% was input from the four inlet sections (only two sections are shown in Fig. 4 because of the plane symmetry) at the side wall. In one inlet section, there were nozzle orifices of 2 mm diameter in 5 places, and C2H2 gas input was performed from the respective nozzle orifices in the −x and x directions, the −y and y directions, and the z direction. The flow rates were set to be 3 L/min per inlet section. Hence, the input from each nozzle orifice was thus 0.6 L/min. Drainage was performed from one position in a different plane from the inlet section. The internal temperature of the reactor was 1223 K, the internal pressure was 1200 Pa, and the carburizing time was 240 s. After the processing, the C2H2 gas supply was stopped immediately followed by dropping into an oil cooling bath outside the reactor for rapid cooling. Then, following test piece cutting in half in the height direction and resin embedding, the alumina-polished surfaces were analysed for the carbon concentration distribution on an electron probe micro analyser (EPMA) (JXA-8200, JEOL, Ltd.). The analysis line of EPMA in the cross section of the test piece is indicated as an allow in Fig. 5. This section was in the position nearest each workpiece, and was presumed to approximately reflect the carburized state of the workpiece.

Fig. 4.

Illustration of the industrial reactor for vacuum carburizing. (Online version in color.)

Fig. 5.

Illustration of the location of test pieces for EPMA. The circles which are indicated by A and B are the test pieces and the arrow in the test piece indicates EPMA line. (Online version in color.)

5. Simulation Conditions

Simulation was performed for the above-described test reactor and industrial reactor. In the computing for the test reactor, the minimum mesh spacing for calculations close to the workpieces was 0.5 mm and a 3D model with a total mesh number of approximately 110000 was constructed. The reactor internal temperature, pressure, flow rates, and other conditions were set to match those of the experiment. In the test reactor, because Γ < 101, we calculated the carbon adsorption rate on the workpiece surfaces through only Step 1 without performing Step 2. In comparing the experimental results, the mean value of the carbon adsorption on the workpiece surfaces was transferred to Step 3, and the total mean carburization amount of the workpiece was computed.

In the computing for the industrial reactor, as shown in Fig. 6, we constructed a 3D model with a total mesh number of approximately 6800 000 where a minimum mesh spacing closed to the workpiece was 1 mm (in Fig. 6, the xy section in the foreground is symmetrical and only half of one side is shown. Note that the view direction is opposite to that in Fig. 4 in terms of z direction). For simplicity, the jig positioning the workpiece and the workpiece gear teeth are not shown. The reactor internal temperature, pressure, flow rate, and other conditions were set to match those of the experiment. For the industrial reactor, because Γ >> 1 was established, Steps 1, 2, and 3 were separately computed. In the simulation, the cylindrical workpiece used for the analysis in the experiment was not considered, and in the comparison with the experimental results, the local carbon adsorption rate of the workpiece surfaces closest to the position of the workpiece used for the analysis in the experiment was computed, and that value was transferred to Step 3 to compute the local carburization amount and the distribution of carbon concentration in steel.

Fig. 6.

Computational domain and grid of the industrial reactor. Note that the view direction is opposite to that in Fig. 4 in terms of z direction. (Online version in color.)

All simulations were performed using the in-house program. Steps 1 and 2 were performed using a thermal-fluid analysis program coupled with chemical reaction based on the finite-volume method15) and Step 3 was computed using the finite-difference method. In Step 3, the calculated mesh spacing used an equal-interval mesh of 7.8 μm.

6. Results and Discussion

6.1. Proposal of Surface Reaction Model

Based on the experimental results of four temperature conditions under tc = 5 s in the test reactor, the pre-exponential factor A and the activation energy E were determined as follows:   

A=7.80× 10 10    [m/s], (28)
  
E=2.65× 10 5    [J/mol/K]. (29)
The comparison of the computational prediction using the propose reaction model to the experimental data for the adsorption rate of carbon as a function of temperature is shown in Fig. 7. In the experiment, the adsorption rate of carbon was given by the weight difference of the test piece before and after the carburizing divided by the carburizing time (5 s) and the surface area of the test piece. Here, it is estimated that the carburizing process is rate-limited by the adsorption rate of carbon onto the steel surface in the first 5 s under these conditions of temperature and pressure (the detail will be discussed later). From Fig. 7, it can be stated that the proposed reaction model enables accurate prediction of the adsorption rate of carbon from C2H2 gas onto the steel surface.
Fig. 7.

Comparison of the computational prediction to the experimental data for adsorption rate of carbon as a function of temperature.

6.2. Variation of Carburized Carbon with Operating Time

The comparison of the computational prediction to the experimental data at Ttr = 1373 K for the amount of carburized carbon as a function of operating time is shown in Fig. 8. In the early stage of carburizing, the amount of carburized carbon increases linearly with the operating time, while the increase rate decreases in the period of tc > 40 s. This early stage corresponds to the period where the carbon concentration does not reach the solubility limit in the steel, and the carburizing is rate-limited by the adsorption rate of carbon onto the steel surface in this period. Whereas, the stage where the increase rate decreases mentioned above is the period after the solubility limit, and the carburizing is rate-limited by the diffusion rate of carbon in the steel in this period. Therefore, the present simulation is incapable of accurately predicting the change of rate-limit and the amount of carburized carbon in the steel.

Fig. 8.

Comparison of the computational prediction to the experimental data for amount of carburized carbon as a function of operating time.

6.3. Distribution of Carbon Concentration in the Industrial Reactor

Figure 9 shows the results of simulation visualizing the xz section from the upward direction (positive y direction) of the industrial reactor. It can be observed that C2H2 gas is diffused toward the centre of the reactor from the nozzles in the inlet section. As a result, a region of high C2H2 gas concentration near the nozzle and a region of low concentration in the centre of the reactor can be observed, and it is thus clear that a significant C2H2 gas concentration distribution forms in the reactor. It can also be confirmed that the carbon adsorption rate on the workpieces varies with the C2H2 gas concentration distribution. Here, to compare the amount of carburizing of the workpieces in position A near the inlet section and those in position B away from the inlet section, Fig. 10 shows the comparison of the carbon concentration distribution in the direction of workpiece depth obtained by simulation and that obtained by EPMA in the experiment. As shown, the carburizing progresses more deeply in position A than in position B. In addition, the simulation accurately predicts these experimental results. Thus, this result indicates that C2H2 gas concentration at position A is higher than at position B and the carbon adsorption rate to the steel surface is therefore higher because position A is near the inlet section as described above. Figure 11 shows the time change in carburization amount per unit area obtained at each position, A and B in Fig. 9, in the simulation. The arrows in the figure indicate the times when the carbon concentration at the steel surface reached the solubility limit. As may be seen, at position A, where the carbon adsorption rate is higher, the surface carbon concentration reached the solubility limit in approximately 10 s and thereafter the carburization rate slows. While at position B, where the carbon adsorption rate is relatively low, the solubility limit is not reached until after approximately 235 s. Therefore, by this simulation, both accurate prediction of the carbon concentration distribution in steel and evaluation of the time change in carburization can be obtained for each workpiece. This promises a capability of application for improvement of the vacuum-carburizing reactor by gas nozzle positioning, optimum workpiece positioning, and by optimizing the process by adjusting the amount of source-gas supply, its concentration, and the carburizing time.

Fig. 9.

Distributions of mass fraction of C2H2 (colours) and adsorption rate of carbon on the surface of workpieces (mono colour). Additionally, some velocity vectors of the gas are plotted.

Fig. 10.

Comparison of the computational predictions to the experimental data for carbon concentration profiles in the steel.

Fig. 11.

Amounts of carburized carbon at positions A and B as a function of operating time.

6.4. Simulation Time

The calculation of the industrial reactor discussed above was conducted on a multi-core processor computer (CPU: Xeon X7560 (Nehalem/8 core/2.27 GHz) × 4, Memory: 128 GB), and the simulation times were around 24 hours for Step 1, 24 hours for Step 2, several minutes for Step 3, and approximately 2 days in total. Therefore, the simulation using the propose technique will complete in several days even in the case of the larger-scale reactor.

7. Conclusion

With a focus on gas convection, gas diffusion, and in-steel carbon diffusion time-scale differences in a vacuum-carburizing reactor, we have proposed a simulation technique for accurate, efficient prediction of gas flow, source gas concentration distribution, and carbon concentration distribution in steel in the reactor.

In regard to the carbon adsorption rate from C2H2 gas onto the steel surface, which is one of the key factors, we used the experimental results in the test reactor and proposed the adsorption reaction model of carbon based on reaction kinetics. In addition, it was shown that with the proposed simulation technique the time change in carburization amount and carbon concentration distribution in selected workpieces in the industrial reactor can be accurately predicted in approximately 2 days.

Gas flow in a reactor can be analysed coupled with carburization depth in steel in a reasonable amount of time by this simulation technique, which therefore promises a capability for application in improvement of a vacuum-carburizing reactor including gas nozzle positioning, optimization of workpiece position, as well as in optimization of process conditions including adjustment of source-gas supply amount, concentration, and carburizing time.

References
 
© 2020 by The Iron and Steel Institute of Japan
feedback
Top