COMSOL 6.4 - Spherically Symmetric Transport

Many models of industrial-transport problems allow the assumption that the problem is spherically symmetric. This assumption is of great importance because it eliminates two space coordinates and leaves a 1D problem that is computationally fast and has very small memory requirements. Some applications where spherical symmetry assumptions are useful include:

•

Reaction and diffusion in catalytic pellets in chemical reactors

•

Heat and mass transfer in the processing of upgraded iron-ore pellets

•

Any other process that takes place in beads that are nearly spherical

For spherical symmetry to be valid, the following assumptions must apply:

•

The computational domain has a spherical shape

•

The outer-perimeter boundary condition does not change with the position on the surface, that is, it does not vary with the space angles θ and φ

•

At any given time for a time-dependent problem, the material properties depend only on the radial distance from the center, r, and not on the space angles θ and φ

•

For a time-dependent problem, the initial condition depends only on the radial distance from the center, r, and not on the space angles θ and φ

Model Definition

The following example simulates the initial transient heating process of a pelletized piece of magnetite ore. This is the first step in the process of making hematite ore pellets, an important raw material for the steel industry.

During the initial heating of a magnetite pellet the temperature is in a range that allows you to disregard any phase change of moisture. Thus it is possible to use a transient heat-conduction equation with constant properties in spherical symmetry. You can also scale the equation for easy parameterization of the radius.

Figure 1 depicts some pellets together with a push pin as a scale reference.

Figure 1: Hematite pellets after drying and oxidation (end product).

The figure shows that these pellets are indeed not perfectly spherical. Nonetheless, this model takes advantage of the assumption of spherical symmetry.

Domain Equations

Starting with the time-dependent heat conduction equation

and expanding it in spherical polar coordinates, the result is the equation

where ρ is the density (kg/m3), cp gives the heat capacity (J/(kg·K)), k denotes the thermal conductivity (W/(m·K)), and Q is an internal heat source (W/m3). Further, r, θ, and φ are the spatial coordinates.

Assuming a perfect sphere of radius Rp and no change in temperature with differing space angles, or ∂T/∂θ = ∂T/∂φ = 0, gives

To avoid division by zero at r = 0, which causes numerical problems, multiply this equation by r2:

(1)

Using a dimensionless radial coordinate

by scaling the equation provides the option to quickly change the pellet’s radius without changing or parameterizing the geometry size¹. Introducing the dimensionless coordinate

and substituting in Equation 1 leads to

(2)

on the following domain:

In a similar manner, it is possible to derive equations similar to Equation 2 for porous media flow, diffusion-reaction problems, and so on.

The model uses the following material data:


Symbol	Name	value
ρ	Density	2000 kg/m3
cp	Heat capacity	300 J/(kg·K)
k	Conductivity	0.5 W/(m·K)
Rp	Pellet radius	0.005 m
Q	Heat source	0 W/m3

Boundary Conditions and Initial Conditions

Because of symmetry about r = 0, there is zero flux through this point, meaning

At the surface,

, you use a convective heating expression with a heat transfer coefficient, hs (W/(m2·K)), for the influx of heat (W/m2):

This expression describes a hot gas with a temperature Text flowing around the pellet. Text is chosen at 95°C. The heat transfer coefficient is set to 1000 W/(m2·K). The initial condition is set to 25°C.

Results

Figure 2: Temperature profiles from t = 0 to t = 10 s.

Figure 2 shows the temperature profiles from 0 to 10 seconds. Each line represents an increment of 0.5 s from the preceding line. From the topmost line it is clear that the center of the pellet has not reached steady state at 10 s. You can also plot the time evolution of the temperature at the center of the pellet.

Figure 3: Time evolution of temperature in the center of a pellet with radius Rp = 5 mm.

Figure 3 shows even more clearly how long the process must yet go before it reaches steady state. An interesting next step is to experiment with different particle radii, Rp, and different heating times.

Figure 4: Time evolution in the center of a pellet with radius Rp = 2.5 mm.

Simply reducing the radius somewhat lets the model reach steady state within 7 s.

Notes About the COMSOL Implementation

To implement Equation 2 and the boundary conditions of this problem, use the 1D time-dependent version of the General Form PDE interface:

The space coordinate in the model is

. For typographical reasons we use rh in this model for “r-hat.” Identifying the general form with Equation 2, the following settings generate the correct equation:


Coefficient	Expression
ea	0
da
Γ (flux vector)
F (source term)	0

You must take special care when setting the heat-influx boundary condition on the pellet surface. hs( Text − T) = k ∂T/∂r on the surface, so you need to compensate G accordingly:

COMSOL_Multiphysics/Equation_Based/spherically_symmetric_transport

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

In the table, enter the following settings:


T

Click

In the tree, select > .

Click

Root

In the window, click the root node.

In the root node’s window, locate the section.

From the list, choose .

The equations in this model are given using the dimensionless radial coordinate. Turning off unit support avoids warnings about inconsistent use of units.

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
rho	2000	2000	Density (kg/m^3)
cp	300	300	Heat capacity (J/(kg*K))
k	0.5	0.5	Thermal conductivity (W/(m*K))
Rp	0.005	0.005	Pellet radius (m)
Qs	0	0	Heat source (W/m^3)
hs	1000	1000	Heat transfer coefficient (W/(m^2*K))
Text	368	368	External temperature (K)
Tinit	298	298	Initial value (K)