- Email: [email protected]
A unified spray forming model for the prediction of billet shape geometry
Acta Materialia 52 (2004) 5275–5288 www.actamat-journals.com
A unified spray forming model for the prediction of billet shape geometry J.H. Hattel a
a,*
, N.H. Pryds
b
Process Modelling Group, Department of Manufacturing Engineering and Management, Technical University of Denmark, Building 425, DK-2800 Lyngby, Denmark b Materials Research Department, Risø National Laboratory, DK-4000 Roskilde, Denmark Received 15 April 2004; received in revised form 7 July 2004; accepted 10 July 2004
Abstract In the present work a unified model for simulating the spray forming process has been developed. Models for the atomization and the deposition processes have been coupled together in order to obtain a new unified description of the spray forming process. The model is able to predict the shape and the temperatures of a spray-formed billet and takes into account the thermal coupling between the gas and the droplets, the change in droplet size distribution along the r-axis in the spray cone and the shading effect. The deposition describes the evolution of the preform with time. For this stage a novel 3D model, which allows the atomizer to be placed asymmetrically over the substrate and also includes the withdrawal of the deposit, was developed. This makes it possible to model not only the growth of a Gaussian shaped preform in which case the spray axis and the rotation axis coincide, but also the surface evolution during billet growth. For this purpose, shading must be taken into account as a core part of the surface evolution algorithm. The unified model involves coupling of three sub models for the atomization, the deposition and the shape of the billet. This coupling, which is a central part of the present work, is also described. Results from the integrated model are presented and the potential for better process understanding as well as process optimization is evident. 2004 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. Keywords: Spray forming; Atomization; Modelling; Rapid solidification; Integrated model
1. Introduction The basic principle of the spray forming process is that molten metal is atomized in an inert or reactive atmosphere to give a spray of liquid particles. The particles are then collected onto a substrate situated below the point of atomization. In the spray forming process, it is possible to spray form in different geometries such as tubes, billets or strips. The most common geometry is a cylindrical billet. This is produced either horizontally or vertically. In or-
*
Corresponding author. Tel.: +45 4525 4710; fax: +45 4593 4570. E-mail address: [email protected] (J.H. Hattel).
der to account for the growth of the billet, the substrate is moved continuously to keep a constant distance between the atomizer and the substrate. Despite continuous contributions in literature towards fundamental understanding of the spray forming process, many important features of the process remain unexplored. An enhanced understanding of the process should permit development of predictive models elucidating the interrelationships between the process/structure/properties and performance. In the literature, different models for the prediction of the spray deposition have been described. Among those are the geometrical model formulated by Frigaard [1,2]. The model describes the dynamics of a spray-formed billet using a geometrical approach and predicts the shape
1359-6454/$30.00 2004 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.actamat.2004.07.016
5276
J.H. Hattel, N.H. Pryds / Acta Materialia 52 (2004) 5275–5288
of the billet produced by both single and multi atomizers. Seok et al. [3,4] have modelled the shape of the developing billet. Their approach is also geometrical and they use both a 2D and a 3D model. Both the models suggested by Frigaard [1,2] and Seok et al. [3,4] are geometrical models in the sense, that they do not take the thermal condition of the spray and other thermal factors into account when predicting the shape of the billet. Models that take both thermal and some geometrical factors into account can be found in literature. A classical model is the one formulated by Mathur et al. [5]. This model takes factors such as both geometrical and thermal sticking efficiency into consideration; however, no attention is paid to the effect of shading. A review of these models can also be found in the references [6,7]. Recently, a 2D model that couples both the thermal condition of the spray, sticking efficiency and mass distribution for a Gaussian shaped preform has been suggested by Hattel et al. [8–11]. In this model the atomizer is allowed to be placed symmetrically over the preform only, making the spray and rotation axes coincide. The objective of the present paper is to describe the latest developed spray forming model. This is a fully coupled, integrated 3D model describing the entire spray forming process of billets. The atomization model takes into account the interaction between the atomization gas and the different droplet sizes. The deposition model, which is coupled to the atomisation model, is able to predict both shape and temperature of the spray-formed billet, based on the coupling of the two models. In order to achieve this, a model describing the effect of shadowing as well as the thermal conditions of the billet has been developed.
2. Atomization model In order to solve correctly the two-phase flow problem, i.e. the gas phase and the liquid metal phase, during the atomisation process, a fully coupled momentum and energy formulation should be applied in general. In the present model, the momentum and energy conservation of the droplets during flight are coupled together via Eq. (1) and the momentum equation, which describes the motion of a spherical droplet in a 1D gas flow. The heat balance between the droplets and the surrounding gas is set up by assuming a 1D Eulerian frame, i.e. fixed finite control volumes along the centreline of the spray cone, assuming that the injected gas is only slightly expanded along the radial direction [12]. This assumption enables to simulate the process in one dimension, see Fig. 1. The total heat balance between the droplet group di and the surrounding gas is used
Fig. 1. Spray forming process modelled by a 1D atomization model and a 3D cylindrical deposition model, including the effect of shading.
to calculate the change in droplet temperature of group Þ, according to the following equation [12]: di, dðT droplet i;j _ global ðd i ÞDðln d i ÞC p dT MP ¼ hAðT d T g Þndroplets AerðT 4d T 4w Þndroplets dfs ndroplets ; þ mdi H f dt
ð1Þ
where ndroplets is the number of droplets having a diam_ the overeter in the interval [di (Dd/2); di + (Dd/2)], M all mass distribution of the droplets, P(di) is assumed to be a log normal distribution of the particles, and the heat transfer coefficient, h, between the droplets and the surrounding gas is given by the Ranz–Marshall equation [13] 2 3 !1=2 cg gg 1=3 kg 4 vr dq 5; h¼ 2 þ 0:6 ð2Þ gg d kg where vr is the relative velocity of the droplet in the gas flow and subscript ÔgÕ denotes gas. The gas velocity is given by the expression [15] vg ¼
Jg expðz=kv Þ; qg Ag
ð3Þ
where the decay constant varies approximately with the initial gas velocity as [15] kv ¼ avn0 ;
ð4Þ
and the two constants, a and n, have been experimentally determined as 3 · 104 and 1.24–1.25, respectively [15].
J.H. Hattel, N.H. Pryds / Acta Materialia 52 (2004) 5275–5288
The mean diameter of the droplet size distribution is directly affected by the gas to metal ratio (GMR) and several suggestions for this relationship have been given in literature. In the present model, the following empirical equation representing the mass median diameter has been used [13]: " #0:5 _ M Dm gl 1þ ¼ Ka ; ð5Þ D0 gg We G_ where Dm is the mass median diameter of the droplets, D0 is the liquid stream diameter, Ka is a constant _ and depending on the experimental conditions and M _G are mass flow rates of the liquid and gas, respectively, and We is the Weber number. From this expression it is readily seen that increasing the GMR decreases the mass median diameter and thereby the droplet size. The temperature of the arriving particles at the deposit surface depend both on the particle size, i.e. smaller droplets cool faster than large ones, as well as on the position of the deposit substrate along the axis of symmetry. This leads to an enthalpy of the spray cone, which is position dependent and influences the surface temperature of the billet. During the solidification of the droplets, the temperature-time/distance curve for each droplet size is characterized by five stages: (1) liquid cooling stage, (2) nucleation, (3) recalescence, (4) continuation of the solidification and (5) solid cooling. Eq. (1) is used in each of these stages (with proper modification according to the stage) to describe the energy balance of the droplets. A more thorough description of the numerical formulations of these stages is given in [4]. The temperature of the surrounding gas, Tg, in Eq. (1) is not prescribed empirically, but is calculated by setting up a heat balance between the heat that is released from different droplet sizes and the gas the following way: n X _ g DT g;j ¼ Qtotal ¼ GC Qgroup i : ð6Þ d i¼1
This enables the droplets to interact with each other via the gas. Recently, Crowe [14,15] suggested a parameter, which helps to evaluate whether the thermal coupling between the surrounding gas and the droplets is necessary. This parameter is based on the conservation of the thermal energy between the gas and the droplets and can be expressed as the change in gas temperature divided by the difference between the gas and the droplet temperature, i.e. _ gÞ _ d =GC ðMC DT g p p T g T d ðNu0 =NuÞðsT uc =LÞ þ 1 _ gÞ _ d =GC ðMC p p ¼ ; 0:5 StT =½ð1 þ 0:3Re Pr0:3 Þ þ 1
ð7Þ
5277
where the thermal Stokes number is defined as StT = sTuc/L, sT is the thermal response time and uc is _ g L2 Þ. the gas velocity in the spray ð G=q Based on this analysis, the choice of applying a thermal coupling between the droplets and the surrounding gas can be made; when the Stokes number is small, the droplets have sufficient time to achieve thermal equilibrium with the gas and in this case the thermal coupling should be included. On the other hand, if the Stokes number is large there is insufficient time for the droplets to transfer heat into the gas and they will exhibit little change in the temperature meaning that the thermal coupling can be neglected. Moreover, the relation between the heat capacity of the gas and the droplets controls the change of gas temperature related to the total temperature difference. As demonstrated [9], in order to estimate the solid fraction of the spray accurately, the whole size distribution must be taken into consideration in a dynamic way so that the different droplet size groups interact with each other via the dynamic calculation of the gas temperature. This thermal coupling between different droplets is often simulated in literature by the weighted average of the solid fraction at a specific flight distance and it is calculated based on a ‘‘summation’’ of discrete events. The discrete event corresponds to individual droplets not interacting with each other and travelling in a gas with a temperature field known a priori to the calculation. A more comprehensive description of the atomization model is given in [9,12,16] as well as an evaluation of the thermal coupling in [9]. For the present calculations, thermal coupling was included.
3. Deposition model In a previous study, the authors proposed a deposition model in which rotational symmetry and that the atomizer is placed perpendicular over the centre of the rotating substrate [8,10] was assumed. This model results in a Gaussian shaped geometry, and it takes into account the thermal component of the sticking efficiency as proposed by Mathur [17]. In the present study, the model was extended to 3D as well as being coupled to a shading algorithm, which allows the atomizer to be placed asymmetrically over the substrate. In order to predict the preform shape when modelling a billet, a geometrical approach is used. This geometric model takes into account geometrical process parameters such as: mass flow and mass distribution, position of the atomizer, distance between atomizer and the preform, preform withdrawal speed and rotation speed of the preform [3,4,18]. Finally, the results from the geometrical model are used in the 3D heat conduction analysis made for the preform.
5278
J.H. Hattel, N.H. Pryds / Acta Materialia 52 (2004) 5275–5288
3.1. Mass distribution The mass distribution used in the present model is an axis symmetrical Gaussian function proposed by Mathur [17] _ M ða; b; rÞ ¼ a expðbr2 Þ; ð8Þ A where r is the shortest distance from a point to the symmetry axis line (spray axis) and a (kg s1 m2) and b (m2) are spray distribution parameters which depend on the total mass flow, the design of the atomizer and the distance from the nearest point on the spray axis to the atomizer (the spray distance). The mass distribution as given by Eq. (8) is shown as function of the spray distances ds and dl in Fig. 2. Due to mass conservation, the mass fluxes through the areas As and Al will be equal. The mass fluxes are obtained by integration of the mass distribution over the proper area Z rs Z 2p _s¼ M as expðbs r2 Þr dh dr 0 0 as 2 ¼ 2p ðebs rs 1Þ; 2bs Z 2p rl _ M l ¼ nt0 al expðbl r2 Þr dh dr 0 al 2 ¼ 2p ðebl rl 1Þ; 2bl
ð9Þ
where rs and rl are the radii of the areas As and Al, respectively. The geometrical relations of ds/dl = rs/rl _s¼M _ l , lead to the foland the conservation of mass, M lowing relation: 2 as b s dl ¼ ¼ : ð10Þ al bl ds
mass distribution at all distances (d1, d2, etc.) by combining Eqs. (10) and (8). The mass flow is then used for calculation of the deposition rate. 3.2. Billet forming mechanism Fig. 3 shows a schematic of the spray forming process and the parameters used in the analyses. The spray axis is inclined an angle, /, from the rotation axis of the substrate. Normally, in the spray forming process the substrate is rotated with an angular velocity, x, and withdrawn with the withdrawal velocity, V, while the atomization point is fixed. In order to simplify the description of the process, the substrate is assumed to be fixed and the atomizer is assumed to move upwards with the velocity, V, and rotate counter clockwise with the angular velocity, x, for convenience. The origin of the reference Cartesian coordinate system, O, is the centre of the substrate; O 0 is the intersection of the spray axis and the rotation axis at time t; At is the position of the atomizer at time t; A0 is the position at the start of the spray forming process; S0 is the intersection point of the spray axis and the substrate at the start of the spray forming process, le is the initial eccentric distance from the point O to S0; ra is the distance from the atomizer to the rotation axis of the substrate; h0 is the initial vertical distance from the atomizer to the substrate; d0 is the axis distance from the atomizer, At, to the intersection, O 0 ; d is the spray distance from the atomizer, At, to point I, which is the nearest point on the spray axis from point P; and r is the distance from point P to point I, which is used to evaluate the droplet flow rate at point P.
Thus, if the mass distribution is known at a certain distance ds (the start distance), it is possible to calculate the
Fig. 2. Schematic diagram of the mass distribution function at different spray distances of ds and dl.
Fig. 3. Schematic diagram of the parameters used when the ! analysing! ! ! spray forming mechanism. d ¼ IAt ; r ¼ IP ; d 0 ¼ O0 At and le ¼ OS 0 .
J.H. Hattel, N.H. Pryds / Acta Materialia 52 (2004) 5275–5288
3.3. Surface evolution
en ¼
The surface at time t = 0 is initially meshed by assigning a number of grid points on the substrate surface. The position of these grid points is then advanced with time as seen in Fig. 5(b). In each time step of the calculation, the new position of all the points on the surface is calculated using the mass flow distribution and geometrical considerations. The new position is found as [4]
ð11Þ
znew ¼ z þ Dh en ez ; where the growth distance, Dh, at point P, during a time interval, Dt, along the surface normal direction en, can be given as follows: Z tþDt _ Dh ¼ nðx; y; z; tÞMða; b; rÞen ef dt; ð12Þ t
where ef is the unit vector along the line between the point P and the atomizer, see Fig. 4. ex, ey and ez are the unit basis vectors of the Cartesians coordinates. n(x, y, x, t) is a function which accounts for the effect of shading (zero in the case of shading, one otherwise). In order to calculate the surface evolution, i.e. Dh, the _ and n must be obtained. The vector four terms, ef, en, M ef is calculated from the following relationship [4]: ef ¼
ap ; ja pj
ð13Þ
where a is the position vector for the atomizer (Fig. 4) and p is the position vector for the point under consideration, P. The surface normal vector, en, can be expressed as
ð14Þ
where n1 and n2 are tangential vectors to the preform surface at point P. Numerically, the surface normal is found by firstly finding the surface composed by the grid point itself and its four neighbouring grid points. The surface is described by an equation of the following form: z ¼ ax2 þ by 2 þ cx þ dy þ e:
xnew ¼ x þ Dh en ex ; y new ¼ y þ Dh en ey ;
n1 n2 ; j n1 n2 j
5279
The normal vector in a point P is then given by oz oz ox ðP Þ; oy ðP Þ; 1 en ðP Þ ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : 2 oz 2 oz ðP Þ þ oy ðP Þ þ 1 ox
ð15Þ
ð16Þ
To ensure continuity between the surfaces, the normals in the neighbouring points are also used when finding the value of the normal in the numerical calculations. As weighting function the following has been found to yield good results: 1 1 en ði; jÞ ¼ en ði; jÞ þ ðen ði þ 1; jÞ þ en ði 1; jÞ 2 8 þ en ði; j þ 1Þ þ en ði; j 1ÞÞ:
ð17Þ
_ as given by Eq. (8) is dependThe mass distribution, M, ent on the radial distance to the spray axis. This distance r, the shortest distance from a point P on the surface to the spray axis can be obtained by vector subtraction as follows: r ¼j i p j;
ð18Þ
where i is the position vector of point I, which is the nearest point on the spray axis from point P. p is the position vector of point P, see Fig. 4. The position vector i can be written as i ¼ a þ des
ð19Þ
where a is the position vector of the atomizer, es is the unit vector for the spray axis and d is the spray distance for point I to the atomizer. Using the orthogonality between the vector es and (i p), ði pÞ es ¼ ða þ des pÞ es ¼ 0
ð20Þ
the spray distance, d, can be obtained as d ¼ es ðp aÞ:
Fig. 4. The geometrical configuration used to calculate surface evolution.
ð21Þ
The distance, d, is then used to find a and b in Eq. (8) using Eq. (10). The position of the atomizer, at, moves upward with the velocity V and it rotates counter clockwise with the angular velocity, x. At time t, the unit vector of the spray axis, es, and the position vector of the atomizer, a, can be calculated as follows [4]:
5280
J.H. Hattel, N.H. Pryds / Acta Materialia 52 (2004) 5275–5288
2
cos xt
6 es ¼ ½ sin /; 0; cos / 4 sin xt 0
sin xt
0
cos xt 0
3
7 0 5; 1 ð22Þ
2
cos xt
6 a ¼ ½ra ; 0; h0 þ V t 4 sin xt 0
sin xt cos xt 0
0
3
7 0 5: 1
ð23Þ
Both es and a are dependent of the position of the atomizer. The matrix in Eqs. (22) and (23) describes the rotation of the atomizer, whereas the vector in Eq. (22) describes the inclination angle from the spray axis and the vector in Eq. (23) describes the position of the atomizer. 3.4. Calculating the shading function, n When modelling spray forming, the effect of shading has to be taken into account. Shading occurs due to the fact that the atomizer ‘‘cannot see’’ every point on the surface, Fig. 1. This results in a process step in which no material is deposited at those points. As mentioned earlier, the shading function n(x, y, z, t) has been introduced in order to describe this phenomenon and it is designed in such a way that n is equal to 0 when shading is present; otherwise it is 1, see Eq. (12). The establishment of the shading function is relatively complex and consists of several calculation steps. The main aspects in the shading procedure are described in the following. The billet surface is divided into triangles of which each is described by three surface points [(xA, yA, zA), (xB, yB, zB) and (xC, yC, zC)], see Fig. 5(a). Initially, when meshing the substrate surface, these triangles are rectangular. Later in the process, the triangles are distorted, Fig. 5(b). The three points form a plane, which can be described as follows: z ¼ C1x þ C2y þ C3:
ð24Þ
The constants C1, C2 and C3 are found by solving the three equations arising from applying the three coordinates of the three points of each triangle in Eq. (24). The planes, which are found, are not limited by the three points, but are infinite in space, Fig. 6. In order to find out whether a plane on the surface is creating shade for a point on the surface, lines are drawn from each point in the mesh (x1, y1, z1) to the atomizer (xatom, yatom, zatom). The equation for the drawn line in Fig. 6 is x x1 y y1 z z1 ¼ ¼ : ð25Þ xatom x1 y atom y 1 zatom z1 The intersection between the line and the plane (x0, y0, z0) is then found by combining Eqs. (24) and (25), and then solving for (x, y, z) = (x0, y0, z0). Once the intersection point (x0, y0, z0) is determined, the following two possibilities exist: (1) the line drawn from the atomization point intersects the plane within the triangle (which is the case in Fig. 6) leading to shading or (2) the line drawn from the atomization point intersects the plane outside the area of the triangle in which case no shading exists from that particular triangle.
Fig. 6. A line from an arbitrary point on the surface to the atomizer. This figure describes a special case in which the line passes through the triangle.
Fig. 5. Surface mesh composed of triangular elements. (a) Notation used in triangle. (b) Initial mesh and the new position of the surface after a period of time.
J.H. Hattel, N.H. Pryds / Acta Materialia 52 (2004) 5275–5288
At each time step, the deposited surface is composed of a large number of plane triangles and the above procedure is then carried out for the intersection point on each plane. If an intersection of a plane leads to no shading, the calculations are moved on to the next plane. In the case where shading is present, i.e. n is 0, the analysis of the remaining planes is bypassed and the calculation moves on to the next point (x2, y2, z2) on the surface where the entire procedure is repeated. Naturally, this procedure is very time consuming. For example in each time step for a triangle calculation grid of N · N, up to 2 · N4 intersections have to be determined and checked for shading. To reduce the calculation time, a procedure to eliminate triangles on which the shading algorithm has to be carried out, was designed. The procedure uses the following principle, see Fig. 7. If the length of the longest side in a triangle on the surface, |AB|, is shorter than the longest distance from the three corner points of the triangle to the intersection point of the plane composed by the corners of the triangle, the intersection point will always be outside the triangle and no shading will occur, hence the shading algorithm can be bypassed. The grey colour marks the area in which the intersection will never be. This area is found by constructing three circles with the radius of |AB| and centres in the three corner points of the triangle (A, B, C). The inter-
5281
Fig. 8. (a) Shape from Seok et al. [4]. (b) Calculated shape with the same process parameters. The process parameters are given in Table 1.
Fig. 9. (a) Shape from Seok et al. [4]. (b) Calculated shape with the same process parameters. The process parameters are given in Table 1.
section of the three circles marks the area in which the intersection of the plane and the line is further than |AB| away from one of the corner points of the triangle. That is, jABj > jBCj > jACj ^ jABj < maxðj AX 0 j; jBX 0 j; jCX 0 jÞ ) No shading possible; ð26Þ where X0 is the intersection between plane and line. Calculated billet profiles as predicted by the shape model have been compared with profiles predicted by the numerical model suggested Seok et al. [3,4], Figs. 8 and 9. The calculations of these shapes are based on the geometrical process parameters given in Table 1. As seen from these figures, the calculated heights are similar to the results obtained by Seok et al. but the shapes differ slightly, which probably is due to the different approaches regarding the shading function.
Fig. 7. If the longest side length in the triangle is shorter than the longest distance from the corner points to the intersection, the intersection is not within the grey area, and hence not within the triangle. Thus, no shading is present in point X0 from that particular triangle.
4. Thermal model in the preform The thermal model for the preform is based on a 3D cylindrical description, i.e. (r, h, z). Thus, the heat
Table 1 Geometrical parameters used to calculate the billet shapes in Figs. 8 and 9
Mass distribution parameter, a (m/s) Mass distribution parameter, b (m2) Withdrawal velocity, V (m/s) Angular velocity, x (deg/s) Inclination angle between the spray and rotation axes, / (deg) Initial distance from atomixer to substrate, d0 (m) Initial excentric distance, le (m)
Fig. 8
Fig. 9
0.01 1000 0.0004 700 35 0.5 0.05
0.025 1000 0.0004 700 60 0.5 0.05
5282
J.H. Hattel, N.H. Pryds / Acta Materialia 52 (2004) 5275–5288
conduction inside the preform is governed by the following equation: qcp
oT o oT 1 oT 1 o oT ¼ k þ 2 k þ k ot or or r or r oh oh
o oT 000 k þ þ Q_ ; oz oz
Rrþ ¼ Rr ¼
ð27Þ
ln
rijk þDr=2 rijk
DhDzk r þDr=2 rijk ln rijk Dr=2 þ ln i1jk ri1jk
DhDzrDrqcp tþDt T ijk T tijk Dt tþDt T tþDt i1jk T ijk ¼ RrðijkÞ þ Rrþði1jkÞ þ M rðijkÞ þ
tþDt T iþ1jk T tþDt ijk RrþðijkÞ þ Rrðiþ1jkÞ M rþðijkÞ
þ
tþDt T tþDt ijþ1k T ijk RhþðijkÞ þ Rhðijþ1kÞ þ M hþðijkÞ
þ
tþDt T tþDt ij1k T ijk RhðijkÞ þ Rhþðij1kÞ þ M hðijkÞ
þ
tþDt tþDt T ijkþ1 T ijk RzþðijkÞ þ Rzðijkþ1Þ þ M zþðijkÞ
þ
tþDt tþDt T ijk1 T ijk ; RzðijkÞ þ Rzþðijk1Þ þ M zðijkÞ
ð28Þ
where the thermal resistances are given by the following expressions [19]:
DhDzk
Rhþ ¼ Rh ¼
rijk Dh ; DrDzk
Rzþ ¼ Rz ¼
Dz ; rijk DhDrk
000
where Q_ is the source term (W m3) accounting for release of latent heat during continuation of solidification in the preform. In the present work, the numerical solution of Eq. (27) is based on an approach similar to the finite volume method, in which the governing differential equation is integrated over each of the control volumes. This way, the conservation principle, i.e. energy balance is expressed for the control volume. In the present case, however, the chosen formulation is slightly different in the sense that it bypasses the volume integration of the differential equation and uses directly the underlying conservation principle together with the concept of thermal resistance, resulting in the governing discretization equation. This is of course only possible due to the simple, structured mesh. In a general unstructured mesh framework, the integration of Eq. (27) must be performed together with application of GaussÕ divergence theorem in order to obtain the discretization equation. However, in the present case the more straightforward resistance formulation is employed. Integration of FourierÕs law yields the expressions for the thermal resistances and the final discretization equation using an implicit time integration scheme, can be written as
r
þ ln riþ1jk ijk Dr=2
;
; ð29Þ
where rijk is the radius of the centre of the actual volume. Note that the release of the remaining latent heat is accomplished by introducing an adjusted value for the specific heat, cp , according to the well-known expression 000 cp ¼ cp H f ðofs =oT Þ (this can be seen by setting Q_ equal to the released heat per unit volume from solidification in Eq. (27)) in combination with an algorithm pushing the temperature back to the liquidus temperature when entering the solidification interval [19]. This to some extent resembles the enthalpy method. Consequently, no iterations are needed for this non-linearity in order to ensure consistency between temperatures and material data. M (with the proper subscript, see Eq. (28)) is the interface resistance between the control volumes, for example due to heat transfer between the preform and the substrate. Another boundary condition which is required by the deposition model is the heat transfer between the preform substrate and the gas. For this purpose the gas should be assigned a temperature. The simplest way to do this is to use the temperature from the droplet module, which is dependent on the distance from the atomizer only, and assign it to all cells with the same distance to the atomizer. But this does not provide any radial dependence of the gas temperature owing to the fact that the atomization model is 1D. In reality, the gas temperature is dependent on the radial distance from the spray axis [20]. Fristching et al. [20] shows that the gas temperature decreases from maximum in the centre of the spray cone to around 298 K at a distance of 0.05 m. As a first approximation, the gas temperature has been modelled by a linear decrease from the centre to 298 K at a radial distance of 0.05 m. Finally, Eq. (28) written for all interior nodes together with boundary conditions linearized in a proper manner [9], results in an algebraic equation system for the temperatures to be solved in each timestep.
5. Coupling of the three models The three described models: the droplet model, the shape model and the deposition model are coupled in order to get a complete model for the spray forming proc-
J.H. Hattel, N.H. Pryds / Acta Materialia 52 (2004) 5275–5288
ess, which describes the formation of a billet shape. The coupling between the shape model and the deposition model is accomplished as follows: at each time step for the thermal calculation, the shape model calculates the z value for all nodal points in the cells in the cylindrical mesh. The z values are found by first finding the triangle plane of the surface which contains the x and y coordinates for the point and then using the equation for that particular plane to find the z value. In the deposition model, at each time step, the height of the preform as predicted by the shape model, is used. If the height of the preform reaches the height of the centre of a cell, the cell is filled. The temperature of the alloy, with which the cell is filled, is found by coupling the shape model with the droplet model. In reality, during the spray forming process, the powder particles, which arrive at the surface of the deposited material, are distributed according to their size; generally speaking, mainly large particles arrive at the centre of the deposited material while smaller particles at the periphery. To achieve this effect, local droplet size distributions are used [12]. The local droplet size distributions are logarithmic Gaussian distributions as is the total droplet size distribution. Their mean droplet size varies with large mean diameters for the distributions close to the spray axis and smaller mean diameters for distributions far from the spray axis [9]. When constructing these local droplet size distributions, the overall droplet size distribution has to be valid still; this is done by satisfying the following equation: P Global d
¼
kX max k¼1
R 2p R rk þDr=2 P local d;k
0
rk Dr=2
Dðrk Þrk dr dh
Total massflow
:
ð30Þ
That is, the probability for a droplet size, d, in the global droplet size distribution has to be the same as the sum of the probabilities for the same droplet size in the local distributions. In order to take into account the mass distribution too, the local distributions in the expression have been weighted in relation to the mass flow. The local distribution can be written as [9,12] P local d
½lnðdÞ 1 ¼ pffiffiffiffiffiffi exp r 2p
2 lnðd local r;50 Þ 2r2
:
ð31Þ
was determined experimenThe local mean size, tally by collecting droplets in a probe at different radial distances, r, and then analysing the mass mean size of each probe. The following relationship was then chosen as a best overall fit from experimental observations [8]: ¼
d 050
expðr=0:621Þ:
The droplet model provides an array of enthalpies for different sizes of droplets at different distances from the atomizer. These values along with the equation for the droplet size distribution are used in the deposition model to convert enthalpies into temperatures and fraction solidified as a function of radial distance from the spray axis, r, and the distance to the atomizer, d. When filling a cell in the deposition model, the knowledge about r and d from the shape model is used to ‘‘look up’’ the temperature and the fraction solidified, and assign them to the cell. The atomisation module provides the enthalpies for each droplet size. When coupling with the other modules, the local droplet size distributions are used. In order to obtain the average enthalpy in a cell, the following equation is used: P ðkÞ ¼ H
i
H ðd i ÞP local ðd i ÞD lnðd i Þ k ; P local P k ðd i ÞD lnðd i Þ
ð33Þ
i
ðkÞ is the average enthalpy for the cell k, H(di) is where H the enthalpy for the droplets with the diameter di and P local is the mass probability of the droplet sizes in the k local droplet size distribution. This calculation is carried out for each r- and z value of the centres of the (r, h, z) mesh. The spray cone itself is assumed to be rotational symmetric. When converting the enthalpies to temperatures, the calculations are divided into three intervals: (1) cooling of liquid, (2) solidification and (3) solid cooling. 5.1. T0 P T > Tliquidus (cooling of liquid) At the starting point at the liquidus temperature + superheat, i.e. T0, the enthalpy of the droplets is set to 0. As the droplets cool from the initial temperature to liquidus (T0 P T > Tliquidus), the change in enthalpy (per unit mass) is calculated as Z T DH ðT Þ ¼ C liquidus ðT ÞdT p T0
!
d local r;50 ,
d local 50 ðrÞ
5283
ð32Þ
This equation has a free parameter d 050 , which is chosen in such a way that the overall droplet size distribution is still valid.
ðT T 0 Þ; ) DH ðT Þ ¼ C liquidus p
ð34Þ
where T0 = Tliquidus + DTsuperheat. The temperature in this interval is found as T ¼ T0 þ
DH ðT Þ DH ðT Þ ¼ T liquidus þ DT superheat þ liquidus C liquidus C p p
ð35Þ
and the corresponding solid fraction is 0. Note that the change in enthalpy has a negative sign. The limit (Tliquidus) corresponds to an enthalpy change of DH ðT ¼ T liquidus Þ ¼ C liquidus ðT liquidus T 0 Þ: p
ð36Þ
5284
J.H. Hattel, N.H. Pryds / Acta Materialia 52 (2004) 5275–5288
5.2. Tliquidus P T > Tsolidus (solidification)
DH ðT Þ ¼ DH ðT solidus Þ þ
¼ DH ðT liquidus Þ þ
T
C liquidus–solidus ðT ÞdT H f fs ðT Þ p T liquidus
þ
ðT solidus C liquidus–solidus p
T liquidus Þ
T solidus Þ
DH C liquidus ðT liquidus T 0 Þ C liquidus–solidus ðT solidus T liquidus Þ þ H f p p C solidus þ T solidus p
:
ð43Þ
ðT T liquidus Þ H f fs ðT Þ þ C liquidus–solidus p ð37Þ and the fraction solid is expressed by the Scheil equation 1=ðke 1Þ T TM fs ðT Þ ¼ 1 ; ð38Þ T liquidus T M which results in DH ðT Þ ¼ C liquidus ðT liquidus T 0 Þ þ C liquidus–solidus ðT T liquidus Þ p p 1=ðke 1Þ T TM : ð39Þ Hf 1 T liquidus T M
As seen, this equation is a transcendental equation in T hence it is not possible to solve the temperature, T, from the equation, explicitly. It could of course be solved numerically, however, another more simple and robust approach is applied. For the temperature interval of Tliquidus P T > Tsolidus the change in enthalpy, DH(T*), is calculated by Eq. (39) for temperature steps of 1 K and this enthalpy is compared with the enthalpy from the droplet model, DH(T). If the enthalpy is within an interval, the temperature for this interval is used, i.e. if DH ðT Þ P DH ðT Þ > DH ðT 1Þ then T ¼ T : ð40Þ As a lower limit for the interval, the solidus temperature, Tsolidius, is used and the corresponding change in enthalpy is DH ¼ C liquidus ðT liquidus T 0 Þ p þ C liquidussolidus ðT solidus T liquidus Þ H f : p
C solidus ðT p
ð42Þ
and the corresponding temperature is obtained as T¼
ðT liquidus T 0 Þ ) DH ðT Þ ¼ C liquidus p
C solidus ðT ÞdT p
) DH ðT Þ ¼ C liquidus ðT liquidus T 0 Þ p DH f þ
Z
T T solidus
In the next interval during solidification, both cooling of the liquid/solid and the release of latent heat occur. The change in enthalpy is calculated as DH ðT Þ
Z
ð41Þ
It should be noted that when using ScheilÕs model for solidification in Eq. (39), not all latent heat is released at the solidus temperature. The numerical model ensures that by releasing the remaining latent heat at Tsolidus and setting fs equal to 1. 5.3. Tsolidus P T (solid cooling) For temperatures below Tsolidus the change in enthalpy is calculated as
The fraction solid is equal to 1 at all temperatures below Tsolidus. Apart from the already shown validation of the shape model Figs. 8 and 9, the sub-models applied in the present work, i.e. the atomization model the deposition model, have been thoroughly evaluated against analytical solutions and experiments previously in literature [10,12,18,21]. Moreover, the recent 2D unified model has been applied to numerically predict the relationship between the GMR and the evolving surface temperature for a Gaussian shape [9]. The unified model in the present work is a generalization of this 2D model in the sense that it is 3D and it takes shading into account.
6. Model prediction Fig. 10(a)–(c) show the results of the billet shape calculation for the three cases in Table 2 at different times and at GMR = 1. The chosen material was 100Cr16 steel and the atomizing gas was N2, see Tables 3 and 4. The height of the billets as seen from these figures depends on the mass flow rate and increases with increasing mass flow rate. This behaviour can be explained by observing Eq. (8) and realizing that the only geometrical parameter, which is varied, is the mass distribution parameter, a. Looking at the shape of the deposit as time increases, it is seen that the shape of the billet reaches a quasi-stationary state (‘‘billet-like’’ shape) already after 20 s for the sample deposited at 0.2 kg/s (Fig. 10(c)). However, for 0.05 kg/s (Fig. 10(a)) and 0.1 kg/s (Fig. 10(b)) the shape remains Gaussian even after longer deposition time. Figs. 11(a)–(d) and 12 show the temperature profiles of the centre (r, z) cross-section in the billet. The temperature profiles after 10, 20, 30 and 40 s for a billet sprayed with a mass flow of 0.2 kg/s and a GMR of 3 are shown in Fig. 11(a)–(d). As seen from these figures the temperature profiles indicate a heat reservoir in the middle of the billet where the melt is semi solid, i.e. the billet temperature is within the range of 1570–1724 K. However, the heat reservoir, i.e. the semi solid region, decreases in width as the billet grows due to heat transfer to the
J.H. Hattel, N.H. Pryds / Acta Materialia 52 (2004) 5275–5288
5285
0.06 10s 20s 30s 40s
Height [m]
0.04
0.02
0 -0.09
-0.06
-0.03
0
0.03
0.06
0.09
Radial distance [m]
(a)
Height [m]
0.1
10s 20s 30s 40s
0.05
0 -0.09
-0.06
-0.03
0
0.03
0.06
0.09
Radial distance [m]
(b)
0.2 10s 20s 30s 40s
Height [m]
0.15
0.1
0.05
0 -0.09
(c)
-0.06
-0.03
0
0.03
0.06
0.09
Radial distance [m]
Fig. 10. Shape of billet after 10, 20, 30 and 40 s for three calculation cases (GMR = 1) in Table 2. (a) Calculation case 1, (b) calculation case 2 and (c) calculation case 3.
surrounding billet side. The cooling from the bottom is not very significant; though the heat transfer coefficient between the substrate and the billet is twice as large as the heat transfer coefficient on the surface, see Table 5. The part of the substrate which is in contact with the billet, quickly reaches a temperature of around 700–800 K and since the heat transfer coefficient on the bottom of
the substrate is only 1/10 of the heat transfer coefficient on the top of the substrate, see Table 5, not much heat is removed from the substrate. This solidification behaviour of the billet naturally affects the structure of the solidified material and hence the properties of the billet. A possible way of increasing the cooling rate during the process is to enhance the removal of the heat from the
5286
J.H. Hattel, N.H. Pryds / Acta Materialia 52 (2004) 5275–5288
Table 2 Calculation case
1
2
3
_ (kg/s) Melt flow, M Mass distribution parameter, a (m/s) Mass distribution parameter, b (m2) Rotation velocity, x (rot/min)* Spray angle, / () Excentric distance, le (m) Distance to atomizer, D0 (m) Withdrawal velocity, (V) Gas flow G_ (kg/s) *** Constant in Eq. (5), Ka Area of gas delivery nozzles A (m2)
0.05 0.00204 1000 116 35 0.02 0.5
0.1 0.00408 1000 116 35 0.02 0.5
0.2 0.00816 1000 116 35 0.02 0.5
**
**
**
0.05–0.1 50 0.0003
0.1–0.25 50 0.0003
0.2–0.6 50 0.0003
*
The rotation velocity is chosen to 116 rot/min since literature suggests that this should be high enough to ensure rotational symmetry [1]. The withdrawal velocity is chosen in such a way that the distance from the center of the billet to the atomizer at all times is D0 in order to keep the solid fraction of the droplets arriving to the surface of the deposit constant. *** The GMR was varied between 1–2, 1–2.5 and 1–3 in the three cases, respectively. **
Table 3 Material data for 100Cr16 Solidus temperature, Tsolidus (K) Liquidus temperature, Tliquidus (K) Solvent melt temperature, TM (K) Eutectic temperature, Te (K) Density, q (kg/m3) Specific heat of solid, C solidus (J/kg K) p (J/kg K) Specific heat of solid, C liquidus p
Table 4 Material data for N2 1570 1724 1811 1419 7810 640 724
Density, q (kg/m3) Viscosity, m (N s/m2) Specific heat, Cp (J/kg K) Thermal conductivity, k (W/m K)
1.5 0.000032 1000 0.026
Fig. 11. Temperature profiles of an (r, z) cross-section of the billet. The mass flow is 0.2 kg/s and the GMR is 3 (calculation case 3). The profiles are shown after (a) 10, (b) 20, (c) 30 and (d) 40 s.
J.H. Hattel, N.H. Pryds / Acta Materialia 52 (2004) 5275–5288
5287
Fig. 12. Temperature profiles of billets after 40 s sprayed with a mass flow 0.05 kg/s (calculation case 1) and a GMR of (a) 1 and (b) 2.
Table 5 Heat transfer coefficients Heat transfer on the billet surface (W m2 K1) Heat transfer between billet and substrate (W m2 K1) Heat transfer on the bottom of the substrate (W m2 K1)
500 1000 100
substrate and/or from the surroundings by, e.g. watercooling the substrate or cooling the billet surface by a stream of gas. Fig. 12 shows temperature profiles of a billet after 40 s with at mass flow of 0.05 kg/s and a GMR of 1 and 2, respectively. As seen from these profiles the deposition temperatures of the preforms are lower for the process with the higher GMR, which is a direct consequence of the fact that the fraction solid for the case under consideration increases with increasing GMR [22]. Comparing Fig. 12(a)–(b) with Fig. 11(a), it is noticed that the billet height and shape are the same, however, the one in Fig. 11(a) is just achieved in lesser time (10 s as compared to 40 s) with a larger mass flow. It can be seen that the temperature of the billet in Fig. 11(a) (mass flow of 0.2 kg/s and GMR equal to 3) is higher than the temperature of the billets in Fig. 11(a)–(b) (mass flow of 0.05 kg/s and GMR equal to 1 and 2, respectively). This is expected and due to the high mass flow rate in Fig. 11(a). When more material is deposited in lesser time, the heat input in the billet is larger, leading to higher temperatures. Another thing worth noticing is the difference in substrate temperature, which is a natural consequence of the difference in time, i.e. 10 s for Fig. 11(a) and 40 s for Fig. 12(a)–(b).
7. Conclusion In the present work a new unified numerical model for simulating the spray forming process has been developed. The model is able to predict shape and temperature of a spray-formed preform and takes into account such factors as thermal coupling between the gas and the droplets, the change in droplet size distribution
along the r-axis in the spray cone, as well as the shading effect. The deposition analysis describes the evolution of the preform with time and a 3D model, which allows the atomizer to be placed asymmetrically over the substrate and also includes the withdrawal of the deposit, was developed. For this model, a coupled geometrical–thermal approach is used. This takes into account geometrical process parameters such as mass flow and mass distribution, position of the atomizer, distance between atomizer and the preform, withdrawal and rotation speed of the preform as well as thermal factors such as droplet enthalpy and heat transfer in the billet. An important factor in the geometrical analysis is the shading effect, which is the fact that the atomizer cannot see every point on the surface. The shading function n(x, y, z, t) has been introduced for this purpose. Finally, some predications of the model have been presented.
References [1] Frigaard IA. SIAM J Appl Math 1995;55(5):1161–203. [2] Frigaard IA. In: Proceedings of the international conference on spray depostion and melt atomization, SDMA I, Bremen, Germany; June 2000. p. 839–53. [3] Seok H-K, Lee HC, Oh KH, Lee J-C, Lee H-I, Ra HY. Metall Mater Trans A 2000;31:1479–88. [4] Seok H-K, Yeo D-H, Oh KH, Lee H-I, Ra HY. Metall Mater Trans B 1998;29:699–708. [5] Mathur PC, Annavarapu S, Appelian D, Lawley A. Mater Sci Eng A 1991;142:261–76. [6] Grant PS. Prog Mater Sci 1995;39:497. [7] Zhou Y, Wu Y, Lavernia EJ. Int J Non-Equilib Process 1997;10:95. [8] Pedersen TB, Hattel JH, Pryds NH, Pedersen AS, Buchholz M, Uhlenwinkel V. In: Proceedings of the international conference on spray depostion and melt atomization, SDMA I, Bremen, Germany; 2000. p. 813–23. [9] Pryds NH, Hattel JH, Pedersen TB, Thorborg J. Acta Mater 2002;50:4075–91. [10] Hattel JH, Pryds NH, J. Thorborg, T.B. Pedersen. In: Proceedings of the international conference on modelling of casting, welding and advanced solidification processes, MCWASP IX, Aachen, Germany; 2000. p. 167–74.
5288
J.H. Hattel, N.H. Pryds / Acta Materialia 52 (2004) 5275–5288
[11] Hattel, JH, Pryds, NH, Pedersen TB, Pedersen AS. In: Proceedings of the international conference on spray depostion and melt atomization, SDMA I, Bremen, Germany; 2000. p. 803–12. [12] Hattel J, Pryds N, Thorborg J, Ottosen P. Modelling Simulat Mater Sci Eng 1999;7(3):413–30. [13] Ranz WE, Marshal WR. Chem Eng Prog 1952;48:173. [14] Crowe CT, Sommerfeld M, Tsuji Y. Multiphase flows with droplets and particles Boca Raton (FL): CRC Press; 1997. [15] C.T. Crowe. In: K. Bauckhage, V. Uhlenwinkel, U. Fritsching, editors. Proceedings of the international conference on spray deposition and melt atomization, Bremen, Germany; 2000. [16] Pryds NH, Hattel JH, Thorborg J. Modelling Simulat Mater Sci Eng 1999;7(3):431–46. [17] Mathur P, Apelian D, Lawley A. Acta Metall 1989;37:429.
[18] Pedersen TB, Hattel JH, Pryds NH. In: Dinesen AR, Eldrup M, Juul Jensen D, Linderoth S, Pedersen TB, Pryds NH, Schrøder Pedersen A, Wert JA, editors. Proceedings of the 22nd Risø international symposium on materials science; 2001. p. 353–58. [19] Hattel JH. editor. Fundamentals of modelling of casting processes [to be published]. [20] Fritsching U, Zhang H, Bauckhage K. In: Wood JV, editor. Proceedings of the second international conference on spray forming Woodhead Publishing; 1993. p. 35–44. [21] Pedersen TB. Numerical modelling of the spray forming process. Ph.D. Thesis, Dept Man Eng and Manag, DTU, 2003. [22] Hattel JH, Pryds NH. Modelling the influence of the gas to melt ratio on the fraction solid of the surface in spray formed billets, SIMS2004. In: Proceedings of the 45th conference on simulation and modelling, Copenhagen, Denmark, 2004 [accepted].
A unified spray forming model for the prediction of billet shape geometry J.H. Hattel a
a,*
, N.H. Pryds
b
Process Modelling Group, Department of Manufacturing Engineering and Management, Technical University of Denmark, Building 425, DK-2800 Lyngby, Denmark b Materials Research Department, Risø National Laboratory, DK-4000 Roskilde, Denmark Received 15 April 2004; received in revised form 7 July 2004; accepted 10 July 2004
Abstract In the present work a unified model for simulating the spray forming process has been developed. Models for the atomization and the deposition processes have been coupled together in order to obtain a new unified description of the spray forming process. The model is able to predict the shape and the temperatures of a spray-formed billet and takes into account the thermal coupling between the gas and the droplets, the change in droplet size distribution along the r-axis in the spray cone and the shading effect. The deposition describes the evolution of the preform with time. For this stage a novel 3D model, which allows the atomizer to be placed asymmetrically over the substrate and also includes the withdrawal of the deposit, was developed. This makes it possible to model not only the growth of a Gaussian shaped preform in which case the spray axis and the rotation axis coincide, but also the surface evolution during billet growth. For this purpose, shading must be taken into account as a core part of the surface evolution algorithm. The unified model involves coupling of three sub models for the atomization, the deposition and the shape of the billet. This coupling, which is a central part of the present work, is also described. Results from the integrated model are presented and the potential for better process understanding as well as process optimization is evident. 2004 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. Keywords: Spray forming; Atomization; Modelling; Rapid solidification; Integrated model
1. Introduction The basic principle of the spray forming process is that molten metal is atomized in an inert or reactive atmosphere to give a spray of liquid particles. The particles are then collected onto a substrate situated below the point of atomization. In the spray forming process, it is possible to spray form in different geometries such as tubes, billets or strips. The most common geometry is a cylindrical billet. This is produced either horizontally or vertically. In or-
*
Corresponding author. Tel.: +45 4525 4710; fax: +45 4593 4570. E-mail address: [email protected] (J.H. Hattel).
der to account for the growth of the billet, the substrate is moved continuously to keep a constant distance between the atomizer and the substrate. Despite continuous contributions in literature towards fundamental understanding of the spray forming process, many important features of the process remain unexplored. An enhanced understanding of the process should permit development of predictive models elucidating the interrelationships between the process/structure/properties and performance. In the literature, different models for the prediction of the spray deposition have been described. Among those are the geometrical model formulated by Frigaard [1,2]. The model describes the dynamics of a spray-formed billet using a geometrical approach and predicts the shape
1359-6454/$30.00 2004 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.actamat.2004.07.016
5276
J.H. Hattel, N.H. Pryds / Acta Materialia 52 (2004) 5275–5288
of the billet produced by both single and multi atomizers. Seok et al. [3,4] have modelled the shape of the developing billet. Their approach is also geometrical and they use both a 2D and a 3D model. Both the models suggested by Frigaard [1,2] and Seok et al. [3,4] are geometrical models in the sense, that they do not take the thermal condition of the spray and other thermal factors into account when predicting the shape of the billet. Models that take both thermal and some geometrical factors into account can be found in literature. A classical model is the one formulated by Mathur et al. [5]. This model takes factors such as both geometrical and thermal sticking efficiency into consideration; however, no attention is paid to the effect of shading. A review of these models can also be found in the references [6,7]. Recently, a 2D model that couples both the thermal condition of the spray, sticking efficiency and mass distribution for a Gaussian shaped preform has been suggested by Hattel et al. [8–11]. In this model the atomizer is allowed to be placed symmetrically over the preform only, making the spray and rotation axes coincide. The objective of the present paper is to describe the latest developed spray forming model. This is a fully coupled, integrated 3D model describing the entire spray forming process of billets. The atomization model takes into account the interaction between the atomization gas and the different droplet sizes. The deposition model, which is coupled to the atomisation model, is able to predict both shape and temperature of the spray-formed billet, based on the coupling of the two models. In order to achieve this, a model describing the effect of shadowing as well as the thermal conditions of the billet has been developed.
2. Atomization model In order to solve correctly the two-phase flow problem, i.e. the gas phase and the liquid metal phase, during the atomisation process, a fully coupled momentum and energy formulation should be applied in general. In the present model, the momentum and energy conservation of the droplets during flight are coupled together via Eq. (1) and the momentum equation, which describes the motion of a spherical droplet in a 1D gas flow. The heat balance between the droplets and the surrounding gas is set up by assuming a 1D Eulerian frame, i.e. fixed finite control volumes along the centreline of the spray cone, assuming that the injected gas is only slightly expanded along the radial direction [12]. This assumption enables to simulate the process in one dimension, see Fig. 1. The total heat balance between the droplet group di and the surrounding gas is used
Fig. 1. Spray forming process modelled by a 1D atomization model and a 3D cylindrical deposition model, including the effect of shading.
to calculate the change in droplet temperature of group Þ, according to the following equation [12]: di, dðT droplet i;j _ global ðd i ÞDðln d i ÞC p dT MP ¼ hAðT d T g Þndroplets AerðT 4d T 4w Þndroplets dfs ndroplets ; þ mdi H f dt
ð1Þ
where ndroplets is the number of droplets having a diam_ the overeter in the interval [di (Dd/2); di + (Dd/2)], M all mass distribution of the droplets, P(di) is assumed to be a log normal distribution of the particles, and the heat transfer coefficient, h, between the droplets and the surrounding gas is given by the Ranz–Marshall equation [13] 2 3 !1=2 cg gg 1=3 kg 4 vr dq 5; h¼ 2 þ 0:6 ð2Þ gg d kg where vr is the relative velocity of the droplet in the gas flow and subscript ÔgÕ denotes gas. The gas velocity is given by the expression [15] vg ¼
Jg expðz=kv Þ; qg Ag
ð3Þ
where the decay constant varies approximately with the initial gas velocity as [15] kv ¼ avn0 ;
ð4Þ
and the two constants, a and n, have been experimentally determined as 3 · 104 and 1.24–1.25, respectively [15].
J.H. Hattel, N.H. Pryds / Acta Materialia 52 (2004) 5275–5288
The mean diameter of the droplet size distribution is directly affected by the gas to metal ratio (GMR) and several suggestions for this relationship have been given in literature. In the present model, the following empirical equation representing the mass median diameter has been used [13]: " #0:5 _ M Dm gl 1þ ¼ Ka ; ð5Þ D0 gg We G_ where Dm is the mass median diameter of the droplets, D0 is the liquid stream diameter, Ka is a constant _ and depending on the experimental conditions and M _G are mass flow rates of the liquid and gas, respectively, and We is the Weber number. From this expression it is readily seen that increasing the GMR decreases the mass median diameter and thereby the droplet size. The temperature of the arriving particles at the deposit surface depend both on the particle size, i.e. smaller droplets cool faster than large ones, as well as on the position of the deposit substrate along the axis of symmetry. This leads to an enthalpy of the spray cone, which is position dependent and influences the surface temperature of the billet. During the solidification of the droplets, the temperature-time/distance curve for each droplet size is characterized by five stages: (1) liquid cooling stage, (2) nucleation, (3) recalescence, (4) continuation of the solidification and (5) solid cooling. Eq. (1) is used in each of these stages (with proper modification according to the stage) to describe the energy balance of the droplets. A more thorough description of the numerical formulations of these stages is given in [4]. The temperature of the surrounding gas, Tg, in Eq. (1) is not prescribed empirically, but is calculated by setting up a heat balance between the heat that is released from different droplet sizes and the gas the following way: n X _ g DT g;j ¼ Qtotal ¼ GC Qgroup i : ð6Þ d i¼1
This enables the droplets to interact with each other via the gas. Recently, Crowe [14,15] suggested a parameter, which helps to evaluate whether the thermal coupling between the surrounding gas and the droplets is necessary. This parameter is based on the conservation of the thermal energy between the gas and the droplets and can be expressed as the change in gas temperature divided by the difference between the gas and the droplet temperature, i.e. _ gÞ _ d =GC ðMC DT g p p T g T d ðNu0 =NuÞðsT uc =LÞ þ 1 _ gÞ _ d =GC ðMC p p ¼ ; 0:5 StT =½ð1 þ 0:3Re Pr0:3 Þ þ 1
ð7Þ
5277
where the thermal Stokes number is defined as StT = sTuc/L, sT is the thermal response time and uc is _ g L2 Þ. the gas velocity in the spray ð G=q Based on this analysis, the choice of applying a thermal coupling between the droplets and the surrounding gas can be made; when the Stokes number is small, the droplets have sufficient time to achieve thermal equilibrium with the gas and in this case the thermal coupling should be included. On the other hand, if the Stokes number is large there is insufficient time for the droplets to transfer heat into the gas and they will exhibit little change in the temperature meaning that the thermal coupling can be neglected. Moreover, the relation between the heat capacity of the gas and the droplets controls the change of gas temperature related to the total temperature difference. As demonstrated [9], in order to estimate the solid fraction of the spray accurately, the whole size distribution must be taken into consideration in a dynamic way so that the different droplet size groups interact with each other via the dynamic calculation of the gas temperature. This thermal coupling between different droplets is often simulated in literature by the weighted average of the solid fraction at a specific flight distance and it is calculated based on a ‘‘summation’’ of discrete events. The discrete event corresponds to individual droplets not interacting with each other and travelling in a gas with a temperature field known a priori to the calculation. A more comprehensive description of the atomization model is given in [9,12,16] as well as an evaluation of the thermal coupling in [9]. For the present calculations, thermal coupling was included.
3. Deposition model In a previous study, the authors proposed a deposition model in which rotational symmetry and that the atomizer is placed perpendicular over the centre of the rotating substrate [8,10] was assumed. This model results in a Gaussian shaped geometry, and it takes into account the thermal component of the sticking efficiency as proposed by Mathur [17]. In the present study, the model was extended to 3D as well as being coupled to a shading algorithm, which allows the atomizer to be placed asymmetrically over the substrate. In order to predict the preform shape when modelling a billet, a geometrical approach is used. This geometric model takes into account geometrical process parameters such as: mass flow and mass distribution, position of the atomizer, distance between atomizer and the preform, preform withdrawal speed and rotation speed of the preform [3,4,18]. Finally, the results from the geometrical model are used in the 3D heat conduction analysis made for the preform.
5278
J.H. Hattel, N.H. Pryds / Acta Materialia 52 (2004) 5275–5288
3.1. Mass distribution The mass distribution used in the present model is an axis symmetrical Gaussian function proposed by Mathur [17] _ M ða; b; rÞ ¼ a expðbr2 Þ; ð8Þ A where r is the shortest distance from a point to the symmetry axis line (spray axis) and a (kg s1 m2) and b (m2) are spray distribution parameters which depend on the total mass flow, the design of the atomizer and the distance from the nearest point on the spray axis to the atomizer (the spray distance). The mass distribution as given by Eq. (8) is shown as function of the spray distances ds and dl in Fig. 2. Due to mass conservation, the mass fluxes through the areas As and Al will be equal. The mass fluxes are obtained by integration of the mass distribution over the proper area Z rs Z 2p _s¼ M as expðbs r2 Þr dh dr 0 0 as 2 ¼ 2p ðebs rs 1Þ; 2bs Z 2p rl _ M l ¼ nt0 al expðbl r2 Þr dh dr 0 al 2 ¼ 2p ðebl rl 1Þ; 2bl
ð9Þ
where rs and rl are the radii of the areas As and Al, respectively. The geometrical relations of ds/dl = rs/rl _s¼M _ l , lead to the foland the conservation of mass, M lowing relation: 2 as b s dl ¼ ¼ : ð10Þ al bl ds
mass distribution at all distances (d1, d2, etc.) by combining Eqs. (10) and (8). The mass flow is then used for calculation of the deposition rate. 3.2. Billet forming mechanism Fig. 3 shows a schematic of the spray forming process and the parameters used in the analyses. The spray axis is inclined an angle, /, from the rotation axis of the substrate. Normally, in the spray forming process the substrate is rotated with an angular velocity, x, and withdrawn with the withdrawal velocity, V, while the atomization point is fixed. In order to simplify the description of the process, the substrate is assumed to be fixed and the atomizer is assumed to move upwards with the velocity, V, and rotate counter clockwise with the angular velocity, x, for convenience. The origin of the reference Cartesian coordinate system, O, is the centre of the substrate; O 0 is the intersection of the spray axis and the rotation axis at time t; At is the position of the atomizer at time t; A0 is the position at the start of the spray forming process; S0 is the intersection point of the spray axis and the substrate at the start of the spray forming process, le is the initial eccentric distance from the point O to S0; ra is the distance from the atomizer to the rotation axis of the substrate; h0 is the initial vertical distance from the atomizer to the substrate; d0 is the axis distance from the atomizer, At, to the intersection, O 0 ; d is the spray distance from the atomizer, At, to point I, which is the nearest point on the spray axis from point P; and r is the distance from point P to point I, which is used to evaluate the droplet flow rate at point P.
Thus, if the mass distribution is known at a certain distance ds (the start distance), it is possible to calculate the
Fig. 2. Schematic diagram of the mass distribution function at different spray distances of ds and dl.
Fig. 3. Schematic diagram of the parameters used when the ! analysing! ! ! spray forming mechanism. d ¼ IAt ; r ¼ IP ; d 0 ¼ O0 At and le ¼ OS 0 .
J.H. Hattel, N.H. Pryds / Acta Materialia 52 (2004) 5275–5288
3.3. Surface evolution
en ¼
The surface at time t = 0 is initially meshed by assigning a number of grid points on the substrate surface. The position of these grid points is then advanced with time as seen in Fig. 5(b). In each time step of the calculation, the new position of all the points on the surface is calculated using the mass flow distribution and geometrical considerations. The new position is found as [4]
ð11Þ
znew ¼ z þ Dh en ez ; where the growth distance, Dh, at point P, during a time interval, Dt, along the surface normal direction en, can be given as follows: Z tþDt _ Dh ¼ nðx; y; z; tÞMða; b; rÞen ef dt; ð12Þ t
where ef is the unit vector along the line between the point P and the atomizer, see Fig. 4. ex, ey and ez are the unit basis vectors of the Cartesians coordinates. n(x, y, x, t) is a function which accounts for the effect of shading (zero in the case of shading, one otherwise). In order to calculate the surface evolution, i.e. Dh, the _ and n must be obtained. The vector four terms, ef, en, M ef is calculated from the following relationship [4]: ef ¼
ap ; ja pj
ð13Þ
where a is the position vector for the atomizer (Fig. 4) and p is the position vector for the point under consideration, P. The surface normal vector, en, can be expressed as
ð14Þ
where n1 and n2 are tangential vectors to the preform surface at point P. Numerically, the surface normal is found by firstly finding the surface composed by the grid point itself and its four neighbouring grid points. The surface is described by an equation of the following form: z ¼ ax2 þ by 2 þ cx þ dy þ e:
xnew ¼ x þ Dh en ex ; y new ¼ y þ Dh en ey ;
n1 n2 ; j n1 n2 j
5279
The normal vector in a point P is then given by oz oz ox ðP Þ; oy ðP Þ; 1 en ðP Þ ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : 2 oz 2 oz ðP Þ þ oy ðP Þ þ 1 ox
ð15Þ
ð16Þ
To ensure continuity between the surfaces, the normals in the neighbouring points are also used when finding the value of the normal in the numerical calculations. As weighting function the following has been found to yield good results: 1 1 en ði; jÞ ¼ en ði; jÞ þ ðen ði þ 1; jÞ þ en ði 1; jÞ 2 8 þ en ði; j þ 1Þ þ en ði; j 1ÞÞ:
ð17Þ
_ as given by Eq. (8) is dependThe mass distribution, M, ent on the radial distance to the spray axis. This distance r, the shortest distance from a point P on the surface to the spray axis can be obtained by vector subtraction as follows: r ¼j i p j;
ð18Þ
where i is the position vector of point I, which is the nearest point on the spray axis from point P. p is the position vector of point P, see Fig. 4. The position vector i can be written as i ¼ a þ des
ð19Þ
where a is the position vector of the atomizer, es is the unit vector for the spray axis and d is the spray distance for point I to the atomizer. Using the orthogonality between the vector es and (i p), ði pÞ es ¼ ða þ des pÞ es ¼ 0
ð20Þ
the spray distance, d, can be obtained as d ¼ es ðp aÞ:
Fig. 4. The geometrical configuration used to calculate surface evolution.
ð21Þ
The distance, d, is then used to find a and b in Eq. (8) using Eq. (10). The position of the atomizer, at, moves upward with the velocity V and it rotates counter clockwise with the angular velocity, x. At time t, the unit vector of the spray axis, es, and the position vector of the atomizer, a, can be calculated as follows [4]:
5280
J.H. Hattel, N.H. Pryds / Acta Materialia 52 (2004) 5275–5288
2
cos xt
6 es ¼ ½ sin /; 0; cos / 4 sin xt 0
sin xt
0
cos xt 0
3
7 0 5; 1 ð22Þ
2
cos xt
6 a ¼ ½ra ; 0; h0 þ V t 4 sin xt 0
sin xt cos xt 0
0
3
7 0 5: 1
ð23Þ
Both es and a are dependent of the position of the atomizer. The matrix in Eqs. (22) and (23) describes the rotation of the atomizer, whereas the vector in Eq. (22) describes the inclination angle from the spray axis and the vector in Eq. (23) describes the position of the atomizer. 3.4. Calculating the shading function, n When modelling spray forming, the effect of shading has to be taken into account. Shading occurs due to the fact that the atomizer ‘‘cannot see’’ every point on the surface, Fig. 1. This results in a process step in which no material is deposited at those points. As mentioned earlier, the shading function n(x, y, z, t) has been introduced in order to describe this phenomenon and it is designed in such a way that n is equal to 0 when shading is present; otherwise it is 1, see Eq. (12). The establishment of the shading function is relatively complex and consists of several calculation steps. The main aspects in the shading procedure are described in the following. The billet surface is divided into triangles of which each is described by three surface points [(xA, yA, zA), (xB, yB, zB) and (xC, yC, zC)], see Fig. 5(a). Initially, when meshing the substrate surface, these triangles are rectangular. Later in the process, the triangles are distorted, Fig. 5(b). The three points form a plane, which can be described as follows: z ¼ C1x þ C2y þ C3:
ð24Þ
The constants C1, C2 and C3 are found by solving the three equations arising from applying the three coordinates of the three points of each triangle in Eq. (24). The planes, which are found, are not limited by the three points, but are infinite in space, Fig. 6. In order to find out whether a plane on the surface is creating shade for a point on the surface, lines are drawn from each point in the mesh (x1, y1, z1) to the atomizer (xatom, yatom, zatom). The equation for the drawn line in Fig. 6 is x x1 y y1 z z1 ¼ ¼ : ð25Þ xatom x1 y atom y 1 zatom z1 The intersection between the line and the plane (x0, y0, z0) is then found by combining Eqs. (24) and (25), and then solving for (x, y, z) = (x0, y0, z0). Once the intersection point (x0, y0, z0) is determined, the following two possibilities exist: (1) the line drawn from the atomization point intersects the plane within the triangle (which is the case in Fig. 6) leading to shading or (2) the line drawn from the atomization point intersects the plane outside the area of the triangle in which case no shading exists from that particular triangle.
Fig. 6. A line from an arbitrary point on the surface to the atomizer. This figure describes a special case in which the line passes through the triangle.
Fig. 5. Surface mesh composed of triangular elements. (a) Notation used in triangle. (b) Initial mesh and the new position of the surface after a period of time.
J.H. Hattel, N.H. Pryds / Acta Materialia 52 (2004) 5275–5288
At each time step, the deposited surface is composed of a large number of plane triangles and the above procedure is then carried out for the intersection point on each plane. If an intersection of a plane leads to no shading, the calculations are moved on to the next plane. In the case where shading is present, i.e. n is 0, the analysis of the remaining planes is bypassed and the calculation moves on to the next point (x2, y2, z2) on the surface where the entire procedure is repeated. Naturally, this procedure is very time consuming. For example in each time step for a triangle calculation grid of N · N, up to 2 · N4 intersections have to be determined and checked for shading. To reduce the calculation time, a procedure to eliminate triangles on which the shading algorithm has to be carried out, was designed. The procedure uses the following principle, see Fig. 7. If the length of the longest side in a triangle on the surface, |AB|, is shorter than the longest distance from the three corner points of the triangle to the intersection point of the plane composed by the corners of the triangle, the intersection point will always be outside the triangle and no shading will occur, hence the shading algorithm can be bypassed. The grey colour marks the area in which the intersection will never be. This area is found by constructing three circles with the radius of |AB| and centres in the three corner points of the triangle (A, B, C). The inter-
5281
Fig. 8. (a) Shape from Seok et al. [4]. (b) Calculated shape with the same process parameters. The process parameters are given in Table 1.
Fig. 9. (a) Shape from Seok et al. [4]. (b) Calculated shape with the same process parameters. The process parameters are given in Table 1.
section of the three circles marks the area in which the intersection of the plane and the line is further than |AB| away from one of the corner points of the triangle. That is, jABj > jBCj > jACj ^ jABj < maxðj AX 0 j; jBX 0 j; jCX 0 jÞ ) No shading possible; ð26Þ where X0 is the intersection between plane and line. Calculated billet profiles as predicted by the shape model have been compared with profiles predicted by the numerical model suggested Seok et al. [3,4], Figs. 8 and 9. The calculations of these shapes are based on the geometrical process parameters given in Table 1. As seen from these figures, the calculated heights are similar to the results obtained by Seok et al. but the shapes differ slightly, which probably is due to the different approaches regarding the shading function.
Fig. 7. If the longest side length in the triangle is shorter than the longest distance from the corner points to the intersection, the intersection is not within the grey area, and hence not within the triangle. Thus, no shading is present in point X0 from that particular triangle.
4. Thermal model in the preform The thermal model for the preform is based on a 3D cylindrical description, i.e. (r, h, z). Thus, the heat
Table 1 Geometrical parameters used to calculate the billet shapes in Figs. 8 and 9
Mass distribution parameter, a (m/s) Mass distribution parameter, b (m2) Withdrawal velocity, V (m/s) Angular velocity, x (deg/s) Inclination angle between the spray and rotation axes, / (deg) Initial distance from atomixer to substrate, d0 (m) Initial excentric distance, le (m)
Fig. 8
Fig. 9
0.01 1000 0.0004 700 35 0.5 0.05
0.025 1000 0.0004 700 60 0.5 0.05
5282
J.H. Hattel, N.H. Pryds / Acta Materialia 52 (2004) 5275–5288
conduction inside the preform is governed by the following equation: qcp
oT o oT 1 oT 1 o oT ¼ k þ 2 k þ k ot or or r or r oh oh
o oT 000 k þ þ Q_ ; oz oz
Rrþ ¼ Rr ¼
ð27Þ
ln
rijk þDr=2 rijk
DhDzk r þDr=2 rijk ln rijk Dr=2 þ ln i1jk ri1jk
DhDzrDrqcp tþDt T ijk T tijk Dt tþDt T tþDt i1jk T ijk ¼ RrðijkÞ þ Rrþði1jkÞ þ M rðijkÞ þ
tþDt T iþ1jk T tþDt ijk RrþðijkÞ þ Rrðiþ1jkÞ M rþðijkÞ
þ
tþDt T tþDt ijþ1k T ijk RhþðijkÞ þ Rhðijþ1kÞ þ M hþðijkÞ
þ
tþDt T tþDt ij1k T ijk RhðijkÞ þ Rhþðij1kÞ þ M hðijkÞ
þ
tþDt tþDt T ijkþ1 T ijk RzþðijkÞ þ Rzðijkþ1Þ þ M zþðijkÞ
þ
tþDt tþDt T ijk1 T ijk ; RzðijkÞ þ Rzþðijk1Þ þ M zðijkÞ
ð28Þ
where the thermal resistances are given by the following expressions [19]:
DhDzk
Rhþ ¼ Rh ¼
rijk Dh ; DrDzk
Rzþ ¼ Rz ¼
Dz ; rijk DhDrk
000
where Q_ is the source term (W m3) accounting for release of latent heat during continuation of solidification in the preform. In the present work, the numerical solution of Eq. (27) is based on an approach similar to the finite volume method, in which the governing differential equation is integrated over each of the control volumes. This way, the conservation principle, i.e. energy balance is expressed for the control volume. In the present case, however, the chosen formulation is slightly different in the sense that it bypasses the volume integration of the differential equation and uses directly the underlying conservation principle together with the concept of thermal resistance, resulting in the governing discretization equation. This is of course only possible due to the simple, structured mesh. In a general unstructured mesh framework, the integration of Eq. (27) must be performed together with application of GaussÕ divergence theorem in order to obtain the discretization equation. However, in the present case the more straightforward resistance formulation is employed. Integration of FourierÕs law yields the expressions for the thermal resistances and the final discretization equation using an implicit time integration scheme, can be written as
r
þ ln riþ1jk ijk Dr=2
;
; ð29Þ
where rijk is the radius of the centre of the actual volume. Note that the release of the remaining latent heat is accomplished by introducing an adjusted value for the specific heat, cp , according to the well-known expression 000 cp ¼ cp H f ðofs =oT Þ (this can be seen by setting Q_ equal to the released heat per unit volume from solidification in Eq. (27)) in combination with an algorithm pushing the temperature back to the liquidus temperature when entering the solidification interval [19]. This to some extent resembles the enthalpy method. Consequently, no iterations are needed for this non-linearity in order to ensure consistency between temperatures and material data. M (with the proper subscript, see Eq. (28)) is the interface resistance between the control volumes, for example due to heat transfer between the preform and the substrate. Another boundary condition which is required by the deposition model is the heat transfer between the preform substrate and the gas. For this purpose the gas should be assigned a temperature. The simplest way to do this is to use the temperature from the droplet module, which is dependent on the distance from the atomizer only, and assign it to all cells with the same distance to the atomizer. But this does not provide any radial dependence of the gas temperature owing to the fact that the atomization model is 1D. In reality, the gas temperature is dependent on the radial distance from the spray axis [20]. Fristching et al. [20] shows that the gas temperature decreases from maximum in the centre of the spray cone to around 298 K at a distance of 0.05 m. As a first approximation, the gas temperature has been modelled by a linear decrease from the centre to 298 K at a radial distance of 0.05 m. Finally, Eq. (28) written for all interior nodes together with boundary conditions linearized in a proper manner [9], results in an algebraic equation system for the temperatures to be solved in each timestep.
5. Coupling of the three models The three described models: the droplet model, the shape model and the deposition model are coupled in order to get a complete model for the spray forming proc-
J.H. Hattel, N.H. Pryds / Acta Materialia 52 (2004) 5275–5288
ess, which describes the formation of a billet shape. The coupling between the shape model and the deposition model is accomplished as follows: at each time step for the thermal calculation, the shape model calculates the z value for all nodal points in the cells in the cylindrical mesh. The z values are found by first finding the triangle plane of the surface which contains the x and y coordinates for the point and then using the equation for that particular plane to find the z value. In the deposition model, at each time step, the height of the preform as predicted by the shape model, is used. If the height of the preform reaches the height of the centre of a cell, the cell is filled. The temperature of the alloy, with which the cell is filled, is found by coupling the shape model with the droplet model. In reality, during the spray forming process, the powder particles, which arrive at the surface of the deposited material, are distributed according to their size; generally speaking, mainly large particles arrive at the centre of the deposited material while smaller particles at the periphery. To achieve this effect, local droplet size distributions are used [12]. The local droplet size distributions are logarithmic Gaussian distributions as is the total droplet size distribution. Their mean droplet size varies with large mean diameters for the distributions close to the spray axis and smaller mean diameters for distributions far from the spray axis [9]. When constructing these local droplet size distributions, the overall droplet size distribution has to be valid still; this is done by satisfying the following equation: P Global d
¼
kX max k¼1
R 2p R rk þDr=2 P local d;k
0
rk Dr=2
Dðrk Þrk dr dh
Total massflow
:
ð30Þ
That is, the probability for a droplet size, d, in the global droplet size distribution has to be the same as the sum of the probabilities for the same droplet size in the local distributions. In order to take into account the mass distribution too, the local distributions in the expression have been weighted in relation to the mass flow. The local distribution can be written as [9,12] P local d
½lnðdÞ 1 ¼ pffiffiffiffiffiffi exp r 2p
2 lnðd local r;50 Þ 2r2
:
ð31Þ
was determined experimenThe local mean size, tally by collecting droplets in a probe at different radial distances, r, and then analysing the mass mean size of each probe. The following relationship was then chosen as a best overall fit from experimental observations [8]: ¼
d 050
expðr=0:621Þ:
The droplet model provides an array of enthalpies for different sizes of droplets at different distances from the atomizer. These values along with the equation for the droplet size distribution are used in the deposition model to convert enthalpies into temperatures and fraction solidified as a function of radial distance from the spray axis, r, and the distance to the atomizer, d. When filling a cell in the deposition model, the knowledge about r and d from the shape model is used to ‘‘look up’’ the temperature and the fraction solidified, and assign them to the cell. The atomisation module provides the enthalpies for each droplet size. When coupling with the other modules, the local droplet size distributions are used. In order to obtain the average enthalpy in a cell, the following equation is used: P ðkÞ ¼ H
i
H ðd i ÞP local ðd i ÞD lnðd i Þ k ; P local P k ðd i ÞD lnðd i Þ
ð33Þ
i
ðkÞ is the average enthalpy for the cell k, H(di) is where H the enthalpy for the droplets with the diameter di and P local is the mass probability of the droplet sizes in the k local droplet size distribution. This calculation is carried out for each r- and z value of the centres of the (r, h, z) mesh. The spray cone itself is assumed to be rotational symmetric. When converting the enthalpies to temperatures, the calculations are divided into three intervals: (1) cooling of liquid, (2) solidification and (3) solid cooling. 5.1. T0 P T > Tliquidus (cooling of liquid) At the starting point at the liquidus temperature + superheat, i.e. T0, the enthalpy of the droplets is set to 0. As the droplets cool from the initial temperature to liquidus (T0 P T > Tliquidus), the change in enthalpy (per unit mass) is calculated as Z T DH ðT Þ ¼ C liquidus ðT ÞdT p T0
!
d local r;50 ,
d local 50 ðrÞ
5283
ð32Þ
This equation has a free parameter d 050 , which is chosen in such a way that the overall droplet size distribution is still valid.
ðT T 0 Þ; ) DH ðT Þ ¼ C liquidus p
ð34Þ
where T0 = Tliquidus + DTsuperheat. The temperature in this interval is found as T ¼ T0 þ
DH ðT Þ DH ðT Þ ¼ T liquidus þ DT superheat þ liquidus C liquidus C p p
ð35Þ
and the corresponding solid fraction is 0. Note that the change in enthalpy has a negative sign. The limit (Tliquidus) corresponds to an enthalpy change of DH ðT ¼ T liquidus Þ ¼ C liquidus ðT liquidus T 0 Þ: p
ð36Þ
5284
J.H. Hattel, N.H. Pryds / Acta Materialia 52 (2004) 5275–5288
5.2. Tliquidus P T > Tsolidus (solidification)
DH ðT Þ ¼ DH ðT solidus Þ þ
¼ DH ðT liquidus Þ þ
T
C liquidus–solidus ðT ÞdT H f fs ðT Þ p T liquidus
þ
ðT solidus C liquidus–solidus p
T liquidus Þ
T solidus Þ
DH C liquidus ðT liquidus T 0 Þ C liquidus–solidus ðT solidus T liquidus Þ þ H f p p C solidus þ T solidus p
:
ð43Þ
ðT T liquidus Þ H f fs ðT Þ þ C liquidus–solidus p ð37Þ and the fraction solid is expressed by the Scheil equation 1=ðke 1Þ T TM fs ðT Þ ¼ 1 ; ð38Þ T liquidus T M which results in DH ðT Þ ¼ C liquidus ðT liquidus T 0 Þ þ C liquidus–solidus ðT T liquidus Þ p p 1=ðke 1Þ T TM : ð39Þ Hf 1 T liquidus T M
As seen, this equation is a transcendental equation in T hence it is not possible to solve the temperature, T, from the equation, explicitly. It could of course be solved numerically, however, another more simple and robust approach is applied. For the temperature interval of Tliquidus P T > Tsolidus the change in enthalpy, DH(T*), is calculated by Eq. (39) for temperature steps of 1 K and this enthalpy is compared with the enthalpy from the droplet model, DH(T). If the enthalpy is within an interval, the temperature for this interval is used, i.e. if DH ðT Þ P DH ðT Þ > DH ðT 1Þ then T ¼ T : ð40Þ As a lower limit for the interval, the solidus temperature, Tsolidius, is used and the corresponding change in enthalpy is DH ¼ C liquidus ðT liquidus T 0 Þ p þ C liquidussolidus ðT solidus T liquidus Þ H f : p
C solidus ðT p
ð42Þ
and the corresponding temperature is obtained as T¼
ðT liquidus T 0 Þ ) DH ðT Þ ¼ C liquidus p
C solidus ðT ÞdT p
) DH ðT Þ ¼ C liquidus ðT liquidus T 0 Þ p DH f þ
Z
T T solidus
In the next interval during solidification, both cooling of the liquid/solid and the release of latent heat occur. The change in enthalpy is calculated as DH ðT Þ
Z
ð41Þ
It should be noted that when using ScheilÕs model for solidification in Eq. (39), not all latent heat is released at the solidus temperature. The numerical model ensures that by releasing the remaining latent heat at Tsolidus and setting fs equal to 1. 5.3. Tsolidus P T (solid cooling) For temperatures below Tsolidus the change in enthalpy is calculated as
The fraction solid is equal to 1 at all temperatures below Tsolidus. Apart from the already shown validation of the shape model Figs. 8 and 9, the sub-models applied in the present work, i.e. the atomization model the deposition model, have been thoroughly evaluated against analytical solutions and experiments previously in literature [10,12,18,21]. Moreover, the recent 2D unified model has been applied to numerically predict the relationship between the GMR and the evolving surface temperature for a Gaussian shape [9]. The unified model in the present work is a generalization of this 2D model in the sense that it is 3D and it takes shading into account.
6. Model prediction Fig. 10(a)–(c) show the results of the billet shape calculation for the three cases in Table 2 at different times and at GMR = 1. The chosen material was 100Cr16 steel and the atomizing gas was N2, see Tables 3 and 4. The height of the billets as seen from these figures depends on the mass flow rate and increases with increasing mass flow rate. This behaviour can be explained by observing Eq. (8) and realizing that the only geometrical parameter, which is varied, is the mass distribution parameter, a. Looking at the shape of the deposit as time increases, it is seen that the shape of the billet reaches a quasi-stationary state (‘‘billet-like’’ shape) already after 20 s for the sample deposited at 0.2 kg/s (Fig. 10(c)). However, for 0.05 kg/s (Fig. 10(a)) and 0.1 kg/s (Fig. 10(b)) the shape remains Gaussian even after longer deposition time. Figs. 11(a)–(d) and 12 show the temperature profiles of the centre (r, z) cross-section in the billet. The temperature profiles after 10, 20, 30 and 40 s for a billet sprayed with a mass flow of 0.2 kg/s and a GMR of 3 are shown in Fig. 11(a)–(d). As seen from these figures the temperature profiles indicate a heat reservoir in the middle of the billet where the melt is semi solid, i.e. the billet temperature is within the range of 1570–1724 K. However, the heat reservoir, i.e. the semi solid region, decreases in width as the billet grows due to heat transfer to the
J.H. Hattel, N.H. Pryds / Acta Materialia 52 (2004) 5275–5288
5285
0.06 10s 20s 30s 40s
Height [m]
0.04
0.02
0 -0.09
-0.06
-0.03
0
0.03
0.06
0.09
Radial distance [m]
(a)
Height [m]
0.1
10s 20s 30s 40s
0.05
0 -0.09
-0.06
-0.03
0
0.03
0.06
0.09
Radial distance [m]
(b)
0.2 10s 20s 30s 40s
Height [m]
0.15
0.1
0.05
0 -0.09
(c)
-0.06
-0.03
0
0.03
0.06
0.09
Radial distance [m]
Fig. 10. Shape of billet after 10, 20, 30 and 40 s for three calculation cases (GMR = 1) in Table 2. (a) Calculation case 1, (b) calculation case 2 and (c) calculation case 3.
surrounding billet side. The cooling from the bottom is not very significant; though the heat transfer coefficient between the substrate and the billet is twice as large as the heat transfer coefficient on the surface, see Table 5. The part of the substrate which is in contact with the billet, quickly reaches a temperature of around 700–800 K and since the heat transfer coefficient on the bottom of
the substrate is only 1/10 of the heat transfer coefficient on the top of the substrate, see Table 5, not much heat is removed from the substrate. This solidification behaviour of the billet naturally affects the structure of the solidified material and hence the properties of the billet. A possible way of increasing the cooling rate during the process is to enhance the removal of the heat from the
5286
J.H. Hattel, N.H. Pryds / Acta Materialia 52 (2004) 5275–5288
Table 2 Calculation case
1
2
3
_ (kg/s) Melt flow, M Mass distribution parameter, a (m/s) Mass distribution parameter, b (m2) Rotation velocity, x (rot/min)* Spray angle, / () Excentric distance, le (m) Distance to atomizer, D0 (m) Withdrawal velocity, (V) Gas flow G_ (kg/s) *** Constant in Eq. (5), Ka Area of gas delivery nozzles A (m2)
0.05 0.00204 1000 116 35 0.02 0.5
0.1 0.00408 1000 116 35 0.02 0.5
0.2 0.00816 1000 116 35 0.02 0.5
**
**
**
0.05–0.1 50 0.0003
0.1–0.25 50 0.0003
0.2–0.6 50 0.0003
*
The rotation velocity is chosen to 116 rot/min since literature suggests that this should be high enough to ensure rotational symmetry [1]. The withdrawal velocity is chosen in such a way that the distance from the center of the billet to the atomizer at all times is D0 in order to keep the solid fraction of the droplets arriving to the surface of the deposit constant. *** The GMR was varied between 1–2, 1–2.5 and 1–3 in the three cases, respectively. **
Table 3 Material data for 100Cr16 Solidus temperature, Tsolidus (K) Liquidus temperature, Tliquidus (K) Solvent melt temperature, TM (K) Eutectic temperature, Te (K) Density, q (kg/m3) Specific heat of solid, C solidus (J/kg K) p (J/kg K) Specific heat of solid, C liquidus p
Table 4 Material data for N2 1570 1724 1811 1419 7810 640 724
Density, q (kg/m3) Viscosity, m (N s/m2) Specific heat, Cp (J/kg K) Thermal conductivity, k (W/m K)
1.5 0.000032 1000 0.026
Fig. 11. Temperature profiles of an (r, z) cross-section of the billet. The mass flow is 0.2 kg/s and the GMR is 3 (calculation case 3). The profiles are shown after (a) 10, (b) 20, (c) 30 and (d) 40 s.
J.H. Hattel, N.H. Pryds / Acta Materialia 52 (2004) 5275–5288
5287
Fig. 12. Temperature profiles of billets after 40 s sprayed with a mass flow 0.05 kg/s (calculation case 1) and a GMR of (a) 1 and (b) 2.
Table 5 Heat transfer coefficients Heat transfer on the billet surface (W m2 K1) Heat transfer between billet and substrate (W m2 K1) Heat transfer on the bottom of the substrate (W m2 K1)
500 1000 100
substrate and/or from the surroundings by, e.g. watercooling the substrate or cooling the billet surface by a stream of gas. Fig. 12 shows temperature profiles of a billet after 40 s with at mass flow of 0.05 kg/s and a GMR of 1 and 2, respectively. As seen from these profiles the deposition temperatures of the preforms are lower for the process with the higher GMR, which is a direct consequence of the fact that the fraction solid for the case under consideration increases with increasing GMR [22]. Comparing Fig. 12(a)–(b) with Fig. 11(a), it is noticed that the billet height and shape are the same, however, the one in Fig. 11(a) is just achieved in lesser time (10 s as compared to 40 s) with a larger mass flow. It can be seen that the temperature of the billet in Fig. 11(a) (mass flow of 0.2 kg/s and GMR equal to 3) is higher than the temperature of the billets in Fig. 11(a)–(b) (mass flow of 0.05 kg/s and GMR equal to 1 and 2, respectively). This is expected and due to the high mass flow rate in Fig. 11(a). When more material is deposited in lesser time, the heat input in the billet is larger, leading to higher temperatures. Another thing worth noticing is the difference in substrate temperature, which is a natural consequence of the difference in time, i.e. 10 s for Fig. 11(a) and 40 s for Fig. 12(a)–(b).
7. Conclusion In the present work a new unified numerical model for simulating the spray forming process has been developed. The model is able to predict shape and temperature of a spray-formed preform and takes into account such factors as thermal coupling between the gas and the droplets, the change in droplet size distribution
along the r-axis in the spray cone, as well as the shading effect. The deposition analysis describes the evolution of the preform with time and a 3D model, which allows the atomizer to be placed asymmetrically over the substrate and also includes the withdrawal of the deposit, was developed. For this model, a coupled geometrical–thermal approach is used. This takes into account geometrical process parameters such as mass flow and mass distribution, position of the atomizer, distance between atomizer and the preform, withdrawal and rotation speed of the preform as well as thermal factors such as droplet enthalpy and heat transfer in the billet. An important factor in the geometrical analysis is the shading effect, which is the fact that the atomizer cannot see every point on the surface. The shading function n(x, y, z, t) has been introduced for this purpose. Finally, some predications of the model have been presented.
References [1] Frigaard IA. SIAM J Appl Math 1995;55(5):1161–203. [2] Frigaard IA. In: Proceedings of the international conference on spray depostion and melt atomization, SDMA I, Bremen, Germany; June 2000. p. 839–53. [3] Seok H-K, Lee HC, Oh KH, Lee J-C, Lee H-I, Ra HY. Metall Mater Trans A 2000;31:1479–88. [4] Seok H-K, Yeo D-H, Oh KH, Lee H-I, Ra HY. Metall Mater Trans B 1998;29:699–708. [5] Mathur PC, Annavarapu S, Appelian D, Lawley A. Mater Sci Eng A 1991;142:261–76. [6] Grant PS. Prog Mater Sci 1995;39:497. [7] Zhou Y, Wu Y, Lavernia EJ. Int J Non-Equilib Process 1997;10:95. [8] Pedersen TB, Hattel JH, Pryds NH, Pedersen AS, Buchholz M, Uhlenwinkel V. In: Proceedings of the international conference on spray depostion and melt atomization, SDMA I, Bremen, Germany; 2000. p. 813–23. [9] Pryds NH, Hattel JH, Pedersen TB, Thorborg J. Acta Mater 2002;50:4075–91. [10] Hattel JH, Pryds NH, J. Thorborg, T.B. Pedersen. In: Proceedings of the international conference on modelling of casting, welding and advanced solidification processes, MCWASP IX, Aachen, Germany; 2000. p. 167–74.
5288
J.H. Hattel, N.H. Pryds / Acta Materialia 52 (2004) 5275–5288
[11] Hattel, JH, Pryds, NH, Pedersen TB, Pedersen AS. In: Proceedings of the international conference on spray depostion and melt atomization, SDMA I, Bremen, Germany; 2000. p. 803–12. [12] Hattel J, Pryds N, Thorborg J, Ottosen P. Modelling Simulat Mater Sci Eng 1999;7(3):413–30. [13] Ranz WE, Marshal WR. Chem Eng Prog 1952;48:173. [14] Crowe CT, Sommerfeld M, Tsuji Y. Multiphase flows with droplets and particles Boca Raton (FL): CRC Press; 1997. [15] C.T. Crowe. In: K. Bauckhage, V. Uhlenwinkel, U. Fritsching, editors. Proceedings of the international conference on spray deposition and melt atomization, Bremen, Germany; 2000. [16] Pryds NH, Hattel JH, Thorborg J. Modelling Simulat Mater Sci Eng 1999;7(3):431–46. [17] Mathur P, Apelian D, Lawley A. Acta Metall 1989;37:429.
[18] Pedersen TB, Hattel JH, Pryds NH. In: Dinesen AR, Eldrup M, Juul Jensen D, Linderoth S, Pedersen TB, Pryds NH, Schrøder Pedersen A, Wert JA, editors. Proceedings of the 22nd Risø international symposium on materials science; 2001. p. 353–58. [19] Hattel JH. editor. Fundamentals of modelling of casting processes [to be published]. [20] Fritsching U, Zhang H, Bauckhage K. In: Wood JV, editor. Proceedings of the second international conference on spray forming Woodhead Publishing; 1993. p. 35–44. [21] Pedersen TB. Numerical modelling of the spray forming process. Ph.D. Thesis, Dept Man Eng and Manag, DTU, 2003. [22] Hattel JH, Pryds NH. Modelling the influence of the gas to melt ratio on the fraction solid of the surface in spray formed billets, SIMS2004. In: Proceedings of the 45th conference on simulation and modelling, Copenhagen, Denmark, 2004 [accepted].