COMSOL 6.4 - Continuous Casting — Apparent Heat Capacity Method

Continuous Casting — Apparent Heat Capacity Method

This example simulates the process of continuous casting of a metal rod from a molten state (Figure 1). To optimize the casting process in terms of casting rate and cooling, it is helpful to model the thermal and fluid dynamics aspects of the process. To get accurate results, you must model the melt flow field in combination with the heat transfer and phase change. The model includes the phase transition from melt to solid, both in terms of latent heat and the varying physical properties. Continuous Casting — Arbitrary Lagrangian–Eulerian Method is a variant of this model using the boundary condition.

Figure 1: Continuous metal-casting process with a view of the modeled section.

This example simplifies the rod’s 3D geometry in Figure 1 to an axisymmetric 2D model in the rz-plane. Figure 2 shows the dimensions of the 2D geometry.

Figure 2: 2D axisymmetric model of the casting process.

As the melt cools down in the mold it solidifies. The phase transition releases latent heat, which the model includes. Furthermore, for metal alloys, the transition is often spread out over a temperature range. As the material solidifies, the material properties change considerably. Finally, the model also includes the “mushy” zone — a mixture of solid and melted material that occurs due to the rather broad transition temperature of the alloy and the solidification kinetics.

This example models the casting process as being stationary using the Heat Transfer in Fluids interface combined with the Laminar Flow interface.

Model Definition

The process operates at steady state, because it is a continuous process. The heat transport is described by the equation:

where k, Cp, and Q denote thermal conductivity, specific heat, and heating power per unit volume (heat source term), respectively.

As the melt cools down in the mold, it solidifies. During the phase transition, a significant amount of latent heat is released. The total amount of heat released per unit mass of alloy during the transition is given by the change in enthalpy, ΔH. In addition, the specific heat capacity, Cp, also changes considerably during the transition.

As opposed to pure metals, an alloy generally undergoes a broad temperature transition zone, over several kelvin, in which a mixture of both solid and molten material co-exist in a “mushy” zone. To account for the latent heat related to the phase transition, the Apparent Heat capacity method is used through the Heat Transfer with Phase Change domain condition. The half-width of the transition interval, ΔT, is set to 10 K in this case, and represents half the transition temperature span.

This example models the laminar flow by describing the fluid velocity, u, and the pressure, p, according to the equations

where ρ is the density (in this case constant), μ is the viscosity, and κ is the dilatational viscosity (here assumed to be zero). Here, the role of the source term, F, is to dampen the velocity at the phase-change interface so that it becomes that of the solidified phase after the transition. The source term follows from the equation (see Ref. 1)

where α can be seen as the volume fraction of the liquid phase; Amush and ε represent arbitrary constants (Amush should be large and ε small to produce a proper damping); and ucast is the velocity of the cast rod.

Table 1 reviews the material properties in this model.

Table 1: Material Properties.
Property	Symbol	Melt	Solid
Density	ρ (kg/m3)	8500	8500
Heat capacity at constant pressure	Cp (J/(kg·K))	530	380
Thermal conductivity	k (W/(m·K))	200	200
Dynamic viscosity	μ (N·s/m2)	0.0434	-

Furthermore, the melting temperature, Tm, and enthalpy of phase change, ΔH, are set to 1356 K and 205 kJ/kg, respectively.

The model uses the parametric solver in combination with adaptive meshing to solve the problem efficiently. In particular, using an adaptive mesh makes it possible to resolve the steep gradients in the mushy zone at a comparatively low computational cost.

Results and Discussion

The plots in Figure 3 display the temperature and phase distributions, showing that the melt cools down and solidifies in the mold region. Interestingly, the transition zone stretches out toward the center of the rod because of poorer cooling in that area.

Figure 3: Temperature distribution (top) and fraction of liquid phase (bottom) near the inlet part of the cast at a casting rate of 1.6 mm/s.

With the modeled casting rate, the rod is fully solidified before leaving the mold (the first section after the die). This means that the process engineers can increase the casting rate without running into problems, thus increasing the production rate.

The phase transition occurs in a very narrow zone although the model uses a transition half width, ΔT, of 10 K. In reality it would be even more distinct if a pure metal were being cast but somewhat broader if the cast material were an alloy with a wider ΔT.

It is interesting to study in detail the flow field in the melt as it exits the die.

Figure 4: Velocity field with streamlines near the inlet part of the process.

In Figure 4, notice the disturbance in the streamlines close to the die wall resulting in a vortex. This eddy flow could create problems with nonuniform surface quality in a real process. Process engineers can thus use the model to avoid these problems and find an optimal die shape.

To help determine how to optimize process cooling, Figure 5 plots the conductive heat flux. It shows that the conductive heat flux is very large in the mold zone. This is a consequence of the heat released during the phase transition, which is cooled by the water-cooling jacket of the mold. An interesting phenomenon of the process is the peak of conductive heat flux appearing in the center of the flow at the transition zone.

Figure 5: The cooling viewed as conductive heat flux in the domains (top), and through the outer boundary (the cooling zones) after the die (bottom).

Furthermore, by plotting the conductive heat flux at the outer boundary for the process as in the lower plot in Figure 5, you can see that a majority of the process cooling occurs in the mold. More interestingly, the heat flux varies along the mold wall length. This information can help in optimizing the cooling of the mold (that is, the cooling rate and choice of cooling method).

You solve the model using a built-in adaptive meshing technique. This is necessary because the transition zone — that is, the region where the phase change occurs — requires a fine discretization. Figure 6 depicts the final mesh of the model. Notice that the majority of the elements are concentrated to the transition zone.

Figure 6: Close-up of the final computational mesh, resulting from the built-in adaptive technique.

The adaptive meshing technique allows for fast and accurate calculations even if the transition width is brought down to a low value, such as for pure metals.

Reference

1. V.R. Voller and C. Prakash, “A fixed grid numerical modeling methodology for convection — diffusion mushy region phase-change problems,” Int.J.Heat Mass Transfer, vol. 30, pp. 1709–1719, 1987.

Heat_Transfer_Module/Thermal_Processing/continuous_casting_apparent_heat_capacity

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

Click

In the tree, select > .

Click

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file continuous_casting_apparent_heat_capacity_parameters.txt.

Note, in particular, the value of the parameter dT, which represents the parameter ΔT in the Model Definition section. It will apply when you solve with adaptive mesh refinement because that solution stage is not related to a parametric study. It is then crucial that the value of dT matches that of the final parameter step for the parametric solution that is used as the initial solution.

Definitions

Variables 1

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Define the variables by loading the corresponding text file provided.

Locate the section. Click

Browse to the model’s Application Libraries folder and double-click the file continuous_casting_apparent_heat_capacity_variables.txt.

Geometry 1

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.065.

In the text field, type 0.1.

Locate the section. In the text field, type -0.1.

In the toolbar, click

Rectangle 2 (r2)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.0625.

In the text field, type 0.025.

Locate the section. In the text field, type -0.125.

In the toolbar, click

Rectangle 3 (r3)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.11575.

In the text field, type 1.4075.

Locate the section. In the text field, type -1.5725.

Click to expand the section. In the table, enter the following settings:


Layer name	Thickness (m)
Layer 1	0.6
Layer 2	0.4
Layer 3	0.3675

In the toolbar, click

Click the

button in the toolbar.

Polygon 1 (pol1)

In the toolbar, click

In the window for , locate the section.

In the table, enter the following settings:


r (m)	z (m)
0	-0.125
0	-0.165
0.11575	-0.165
0.0625	-0.125
0	-0.125

In the toolbar, click

Click the

button in the toolbar.

This completes the geometry modeling stage.

Materials

Now, add the following two materials to the model, labeled and . The solid metal alloy is used in the feature for the solid phase, while the liquid metal alloy is used for the liquid phase. The liquid metal alloy also defines fluid properties used in the interface.

Solid Metal Alloy

In the toolbar, click

In the window for , type Solid Metal Alloy in the text field.

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Dynamic viscosity	mu	0.0434[Pa*s]	Pa·s	Basic
Heat capacity at constant pressure	Cp	Cp_s	J/(kg·K)	Basic
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	200[W/(m*K)]	W/(m·K)	Basic
Density	rho	8500[kg/m^3]	kg/m³	Basic

Liquid Metal Alloy

In the toolbar, click

In the window for , type Liquid Metal Alloy in the text field.

Locate the section. From the list, choose .

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Dynamic viscosity	mu	0.0434[Pa*s]	Pa·s	Basic
Heat capacity at constant pressure	Cp	Cp_l	J/(kg·K)	Basic
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	200[W/(m*K)]	W/(m·K)	Basic
Density	rho	8500[kg/m^3]	kg/m³	Basic

Laminar Flow (spf)

Initial Values 1

In the window, under > click .

In the window for , locate the section.

Specify the u vector as


0	r
v_cast	z

Inlet 1

In the toolbar, click

and choose .

Select Boundary 15 only.

In the window for , locate the section.

From the list, choose .

Outlet 1

In the toolbar, click

and choose .

Select Boundary 2 only.

In the window for , locate the section.

From the list, choose .

Locate the section. Click the button.

Specify the u0 vector as


0	r
v_cast	z

Wall 2

In the toolbar, click

and choose .

Select Boundaries 20–22 only.

In the window for , locate the section.

From the list, choose .

Volume Force 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. Specify the F vector as


Fr	r
Fz	z

Definitions

Ambient Properties 1 (ampr1)

In the toolbar, click

and choose .

In the window for , locate the section.

In the Tamb text field, type 300[K].

This defines the ambient temperature for heat transfer between the outer surfaces and the surroundings.

Heat Transfer in Fluids (ht)

Fluid 1

In the window, under > click .

Phase Change Material 1

In the toolbar, click

and choose .

In the window for , locate the section.

In the T1 → 2 text field, type T_m.

In the ΔT1 → 2 text field, type 2*dT.

The parameter dT is multiplied by 2 because it is only the half width of the phase change interval.

In the L1 → 2 text field, type dH.

Locate the section. From the list, choose .

Initial Values 1

In the window, under > click .

In the window for , locate the section.

In the T text field, type T_in.

Inflow 1

In the toolbar, click

and choose .

Select Boundary 15 only.

In the window for , locate the section.

In the Tustr text field, type T_in.

Heat Flux 1

In the toolbar, click

and choose .

Select Boundary 23 only.

In the window for , locate the section.

From the list, choose .

In the h text field, type h_br.

From the Text list, choose .

Heat Flux 2

In the toolbar, click

and choose .

Select Boundary 22 only.

In the window for , locate the section.

From the list, choose .

In the h text field, type h_mold.

From the Text list, choose .

Heat Flux 3

In the toolbar, click

and choose .

Select Boundaries 20 and 21 only.

In the window for , locate the section.

From the list, choose .

In the h text field, type h_air.

From the Text list, choose .

Outflow 1

In the toolbar, click

and choose .

Select Boundary 2 only.

Surface-to-Ambient Radiation 1

In the toolbar, click

and choose .

Select Boundaries 20 and 21 only.

In the window for , locate the section.

From the ε list, choose . In the associated text field, type eps_s.

From the Tamb list, choose .

Mesh 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Size

In the toolbar, click

Boundary Layers 1

In the window, right-click and choose .

Click to confirm.

Deleting the boundary layers is necessary in order to use the adaptive mesh functionality.

Size 1

In the window, under > click .

Select Boundaries 16–21 only.

In the window for , locate the section.

From the list, choose .

In the window, right-click and choose .

Study 1

Compute the solution using a three-step process. First, solve the problem using dT as a continuation parameter with the parametric solver on the default mesh, gradually decreasing the value of dT. Then, use the adaptive solver to adapt the mesh. Finally, use the parametric solver again to decrease dT further down to a value of 10 K.

In the window, click .

In the window for , locate the section.

Clear the checkbox.

The default plots are disabled for this study because they will be added from the last study.

Step 1: Stationary

In the window, under click .

In the window for , click to expand the section.

Select the checkbox.

Click

In the table, enter the following settings:


Parameter name	Parameter value list	Parameter unit
dT (Temperature transition zone half width)	300 100 50 30	K

Step 2: Stationary 2

In the toolbar, click

In the window for , click to expand the section.

From the list, choose .

Solution 1 (sol1)

In the toolbar, click

In the window, expand the node.

In the window, expand the > > > node, then click .

In the window for , locate the section.

Clear the checkbox.

In the toolbar, click

Level 2 Adapted Mesh 2

Before proceeding with the final solution stage, inspect the adapted mesh. You find it under the automatically created branch in the model tree.

In the window, expand the > node, then click .

In the window, expand the node, then click .

Click the button in the toolbar and then use the mouse to zoom in on the transition zone where the mesh is the densest.

The mesh should look like that in Figure 6.

Add a second study for the second parametric study step.

Add Study

In the toolbar, click

to open the window.

Go to the window.

Find the subsection. In the tree, select > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Study 2

In order to get faster convergence, you need to use the previous solution as the initial value for this study.

Step 1: Stationary

In the window, under click .

In the window for , click to expand the section.

Find the subsection. From the list, choose .

From the list, choose .

Locate the section. Select the checkbox.

Click

In the table, enter the following settings:


Parameter name	Parameter value list	Parameter unit
dT (Temperature transition zone half width)	25 20 16 13 10	K

Notice that , the adapted mesh, is the default selection in the mesh list. Keep this setting.

Again, a fully coupled solver is more robust for this model. Tweak the solver sequence accordingly with the instructions below.

Solution 7 (sol7)

In the toolbar, click

In the window, expand the node.

Right-click > > > and choose .

In the toolbar, click

Results

Velocity (spf)

To reproduce the plot in Figure 4, plot the velocity field as a combined surface and streamline plot.

Surface

In the window, expand the node, then click .

In the window for , locate the section.

From the list, choose .

Velocity (spf)

In the window, click .

Streamline 1

In the toolbar, click

In the window for , locate the section.

From the list, choose .

In the text field, type 12.

Locate the section. Find the subsection. From the list, choose .

Click to expand the section. In the toolbar, click

Pressure (spf)

The second default plot shows the pressure profile in the 2D slice.

Velocity, 3D (spf)

The third default plot shows the velocity magnitude in 3D obtained by revolution of the 2D axisymmetric dataset.

Temperature (ht)

This default plot shows the temperature in the 2D slice.

Proceed to reproduce the lower plot in Figure 3, showing the fraction of liquid phase.

Fraction of Liquid Phase

In the toolbar, click

In the window for , type Fraction of Liquid Phase in the text field.

Locate the section. From the list, choose .

Surface 1

In the toolbar, click

In the window for , click in the upper-right corner of the section. From the menu, choose > > > .

In the toolbar, click

Notice, in particular, the narrow transition zone between the two phases.

Temperature (ht)

To reproduce the upper plot in Figure 3, which visualizes the temperature and velocity fields, proceed as follows.

Duplicate the default temperature plot to create a new plot combining temperature and flow streamlines.

Temperature and Streamlines

In the window, right-click and choose .

In the window for , type Temperature and Streamlines in the text field.

Streamline 1

In the toolbar, click

In the window for , locate the section.

From the list, choose .

In the text field, type 2.2.

In the text field, type 11.4.

Locate the section. Find the subsection. From the list, choose .

In the toolbar, click

Proceed to reproduce the heat flux plots shown in Figure 5.

Conductive Heat Flux

In the toolbar, click

In the window for , type Conductive Heat Flux in the text field.

Locate the section. From the list, choose .

Surface 1

In the toolbar, click

In the window for , click in the upper-right corner of the section. From the menu, choose > > > .

In the toolbar, click

Conductive Heat Flux

In the window, click .

Streamline 1

In the toolbar, click

In the window for , locate the section.

From the list, choose .

In the text field, type 2.2.

In the text field, type 11.4.

Locate the section. Find the subsection. From the list, choose .

In the toolbar, click

The following steps reproduce the lower plot in the same figure, showing the conductive heat flux through the outer boundaries.

Conductive Heat Flux through Outer Boundaries

In the toolbar, click

In the window for , type Conductive Heat Flux through Outer Boundaries in the text field.

Locate the section. From the list, choose .

From the list, choose .

Click to expand the section. From the list, choose .

In the text area, type Conductive heat flux through outer boundaries.

Locate the section.

Select the checkbox. In the associated text field, type z-coordinate (m).

Select the checkbox. In the associated text field, type Normal conductive heat flux (W/m^2).

Line Graph 1

In the toolbar, click

Select Boundaries 20–23 only.

In the window for , click in the upper-right corner of the section. From the menu, choose > > > .

Click in the upper-right corner of the section. From the menu, choose > > > .

Click to expand the section. From the list, choose .

From the list, choose .

Click to collapse the section. In the toolbar, click

Compare the result with the lower plot of Figure 5.

Finally, verify that the final mesh is sufficiently fine to resolve the temperature-dependence of the latent heat.

Cut Line 2D 1

In the toolbar, click

In the window for , locate the section.

In row , set to 0.045, and to -0.42.

In row , set to 0.085, and to -0.43.

These values are chosen such that the two points are on opposite sides of and approximately perpendicular to the transition zone.

Alternatively, you can select the two end points and create the Cut Line 2D dataset with the help of the plot. Select the plot node, then click in the window after first selecting, in turn, and in the main toolbar.

Locate the section. From the list, choose .

Temperature Dependence, Latent Heat

In the toolbar, click

In the window for , type Temperature Dependence, Latent Heat in the text field.

Locate the section. From the list, choose .

Line Graph 1

In the toolbar, click

In the window for , click in the upper-right corner of the section. From the menu, choose > > > .

Click to expand the section. Select the checkbox.

In the toolbar, click

As you can see, the curves for the lower ΔT values, in particular at 10 K, are not entirely smooth. Thus, if you were to reduce ΔT further to model the casting of some pure metal, you would need to increase the mesh resolution.