enable_heat_transfer command

Purpose

Command for calculating particle-particle and particle-wall heat transfer and temperature update for the particle(s) and/or wall (via mesh_module_heattransfer).

Note

This command is supported by Aspherix GPU.

Syntax

enable_heat_transfer keyword value

Keywords:

Keyword	Description
initial_particle_temperature	obligatory, initial temperature of the particles units: [temperature]
id	user-assigned name for the command call
particle_group	ID of the group of particles on which the command is applied default: all
store_contact_data	available options: `yes` or `no` default: `no`
maximum_particle_temperature	maximum temperature of the particles units: [temperature]
number_of_shells	number of shells used to discretize each particle default: 1
enable_radiation	available options: `yes` or `no` default: `no`
add_external_heatflux	available options: `yes` or `no` default: `no`
add_external_heatflux_from_grid	available options: `yes` or `no` default: `no`
output_detailed_heatflux	available options: `yes` or `no` default: `no`

Keyword	Description
view_factors_file	filename of the view factors file
emissivity	emissivity of the particles units [power/length^2]
max_relative_radius	relative cutoff radius for the radiation omdel units [-]
interpolation_resolution	number of interpolation points in the view factors data default: 20
apparent_particle_radius	particle radius taken for area calculation units [length]

Examples

enable_heat_transfer initial_particle_temperature 273.15
enable_heat_transfer particle_group my_particles id my_transfer initial_particle_temperature 273.15 &
  maximum_particle_temperature 373.15
enable_heat_transfer initial_particle_temperature 273.15 store_contact_data yes
enable_heat_transfer initial_particle_temperature 400 number_of_shells 20
enable_heat_transfer initial_particle_temperature 273.15 enable_radiation yes view_factors_file vf.csv emissivity 0.65 max_relative_radius 4 interpolation_resolution 15 add_external_heatflux yes

Associated material properties

Material properties

thermalConductivity ( $k$ ): thermal conductivity of a material [power / length temperature]
thermalCapacity ( $c$ ): thermal (specific) capacity of a material [energy / mass temperature]
youngsModulusOriginal ( $Y^*_{orig}$ ): original Young’s modulus of each material [pressure] (required if area_correction yes is used in the particle_contact_model command)

Material interaction properties

ht_modification ( $\alpha$ ): modifier for contact area calculation [-] (required if modified_area_correction yes is used in the particle_contact_model command)

Description

This command computes the particle-particle and particle-wall heat transfer due to interaction and updates the particle temperature. To control the behaviour of this model several keywords can be used and their behaviour is detailed below.

The obligatory initial_particle_temperature defines the initial temperature of the particles.

The temperature inside the particle can be updated using two different models. The default uniform temperature model assumes a constant temperature throughout the particle and can be chosen by using number_of_shells 1. The alternative shell model discretizes the particle into a number of shells specified by number_of_shells N (N > 1) and computes the corresponding heat transfer inside the particle.

This model chosen with number_of_shells 1 assumes a constant temperature inside a particle and updates the temperature based on the following model [1]

$m_i c_{p,i} \frac{d T_i}{d t} = \sum_{j} q_{ij} + q_{src,i},$

where $m_i$ is the mass, $c_{p,i}$ is the heat capacity, $T_i$ is the temperature and $q_{src,i}$ is the heat source (or sink) for the i-th particle. The heat flux $q_{ij}$ stemming from conduction via inter-particle contacts is summed over all contacts (index j) of particle i.

This model discretizes the interior of each particle into N radial shells and computes the heat transfer between these shells inside each particle.

Heat conduction and temperature update are solved according to following 1D heat equation

$\frac{\partial T}{\partial t} = D \left(\frac{2}{r}\frac{\partial T}{\partial r} + \frac{\partial^2 T}{\partial r^2}\right)$

where the temperature inside a particle does depend only on the radial distance, but not on the azimuthal and polar angles. D is the thermal diffusivity:

$D = \frac{\lambda}{\rho C_p}$

Dividing the spherical particle of radius $r$ into $N$ spherical shells and discretizing the heat equation using the implicit integration scheme gives us

$\frac{T_k^{n+1} - T_k^n}{\Delta t} = D\left(\frac{T_{k+1}^{n+1} - T_k^{n+1}}{r \Delta r} + \frac{T_{k+1}^{n+1} - 2 T_k^{n+1} + T_{k+1}^{n+1}}{\Delta r^2}\right),\; k = 1,\ldots,N \quad (*)$

where the temperature field is defined at the border points between shells (grid points) and $k$ is the number of the grid point. Note that this scheme is valid only for internal grid points ( $0<k<N$ ).

In order to avoid singularity problem for the particle center ( $r=k=0$ ) the L’Hopital’s rule is applied to the original heat equation, leading to

$\frac{\partial T}{\partial t} = 3D \frac{\partial^2 T}{\partial r^2},\; r = 0 \quad (**)$

Applying the auxiliary point method:

$T_{k+1}^{n+1} = T_{k-1}^{n+1},\;k = 0$

the discretization of Equation (**) for $r = 0$ leads to

$\frac{T_k^{n+1} - T_k^n}{\Delta t} = 6D \frac{T_{k+1}^{n+1} - T_k^{n+1}}{\Delta r^2},\; k = 0$

For the outer grid point ( $k=N$ ), Neumann boundary condition is applied and discretized by introducing the auxiliary point $k=N+1$ :

$-\lambda \frac{T_{k+1}^{n+1} - T_k^{n+1}}{\Delta r} = \sum_{j} q_{ij} + q_{src,i},\; k = N$

where $q_{src,i}$ is the heat source (or sink) for the i-th particle. The heat flux $q_{ij}$ stemming from conduction via inter-particle contacts is summed over all contacts (index j) of particle i. The value for the auxiliary point ( $k+1$ ) can be expressed and substituted in Equation (*).

As a result, the above mentioned equations provide a system of $N+1$ linear equations with $N+1$ unknown values for temperature inside a particles at the next time step:

$b_0 T_0^{n+1} + c_0 T_0^{n+1} = d_0,\; k = 0 a_k T_{k-1}^{n+1} + b_k T_k^{n+1} + c_k T_{k+1}^{n+1} = d_k,\; 1 \leq k \leq N-1 a_N T_N^{n+1} + b_N T_N^{n+1} = d_N,\; k = N$

This system is then solved using the tridiagonal matrix algorithm.

To make particles adiabatic (so they do not change the temperature), do not include them in the fix group. However, heat transfer is calculated between particles in the group and particles not in the group (but temperature update is not performed for particles not in the group). Thermal conductivity and specific thermal capacity must be defined for each atom type used in the simulation by means of fix property/global commands:

fix id all property/global thermalConductivity peratomtype value_1 value_2 ...
(value_i=value for thermal conductivity of atom type i)

fix id all property/global thermalCapacity peratomtype value_1 value_2 ...
(value_i=value for thermal capacity of atom type i)

To set the temperature distribution for a group of particles, you can use the set command with keyword property/atom and values temp_array T_0 T_1 .. T_N. T_0 T_1 .. T_N are the temperature values for $N+1$ grid points respectively you want the particles to have. To set heat sources (or sinks) for a group of particles, you can also use the set command with the set keyword: property/atom and the set values: heatSource h where $h$ is the heat source value you want the particles to have (in Energy/time units). A negative value means it is a heat sink. Examples would be:

set region halfbed property/atom temp_array 800. 800. 800. 800. 800. 800. 800. 800. 800. 800. 800. 800. 800. 800. 800. 800. 800. 800. 800. 800. 800.
set region srcreg property/peratom heatSource 0.5

The heat flux between two particles or a particle and a wall in contact is defined as

$q_{ij} = h_{ij} \left( T_j - T_i \right),$

where the heat transfer coefficient $h_{ij}$ is calculated from the individual thermal conductivities $k$ and a contact area $A_{c,ij}$ as

$h_{ij} = H(k_i, k_j, A_{c,ij}),$

where $H$ depends on whether a particle-particle or particle-wall contact is considered.

The heat transfer coefficient for particle-particle interactions $h_{ij}$ is calculated from the individual thermal conductivities $k$ as follows

$h_{ij} = 2 \frac{2 k_i k_j}{k_i + k_j} \sqrt{\frac{A_{c,ij}}{\pi}},$

where $A_{c,ij}$ is the contact area between the i-th and the j-th particle. Note that for these relations to be valid, in general the temperatures must be defined far away from the point of contact.

Similarly, the particle-wall heat flux can be written as

$q_{iw} = h_{iw} \left( T_{wall} - T_i \right),$

where the particle-wall heat transfer coefficient $h_{iw}$ for particle i and a wall is calculated as follows:

$h_{iw} = 4 k_i \sqrt{\frac{A_{c,iw}}{\pi}},$

where $k_i$ is the particle thermal conductivity and $A_{c,iw}$ is the contact area between the i-th and the wall. The difference between these two equations for the heat transfer coefficient comes from the definition of the (discrete) temperature locations. For particle-particle interactions these are both far away for the point of contact. However, for particle-wall interactions only one the particle temperature is defined far away, while the wall temperature is defined at the point of contact. Therefore, if the particle and wall would have identical thermal conductivities, these two equations will lead to different heat transfer coefficients ( $2 h_{ij} = h_{iw}$ ). The factor 2 compensates for the distance change between the (discrete) temperature locations, which for particle-particle interactions is 2 times the distance for particle-wall interactions. Another way to look at it is to consider a pile of 2 particles squeezed between two flat horizontal walls with different temperatures. To obtain (vertically) a linear temperature gradient (which implies a constant heat flux value between the walls), the particle-particle and particle-wall temperature differences should be equal to $\left( T_j - T_i \right) = 2 \left( T_{wall} - T_i \right)$ . Therefore, the shown heat transfer coefficients and temperature differences will lead to a constant heat flux value across the domain (it is not be affected by the interaction type: particle-particle or particle-wall): $q_{ij} = q_{iw}$ .

To make particles adiabatic (so they do not change the temperature), do not include them into the group addressed by group-ID. Independently from that, heat transfer is calculated between particles in the group and particles not in the group (but temperature update is not performed for particles not in the group). Thermal conductivity and specific thermal capacity must be defined for each atom type used in the simulation by means of material_properties command:

material_properties glass thermalConductivity 100 thermalCapacity 1 ...

To set the temperature for a group of particles, you can use the set command with keyword property/atom and values Temp T. T is the temperature value you want the particles to have. To set heat sources (or sinks) for a group of particles, you can also use the set command with the set keyword: property/atom and the set values: heatSource $q_{src}$ , where $q_{src}$ is the heat source value you want the particles to have (in Energy/time units). A negative value means it is a heat sink. Examples would be:

set region halfbed property/peratom Temp 800.
set region srcreg property/peratom heatSource 0.5

The calculation of the contact area $A_{c,ij}$ can be influenced by several keyword settings of the particle_contact_model command to adjust the particle shapes, surface properties and other relevant effects.

Please see the documentation of the particle contact model and wall contact model commands. They have several keywords to control the details of the contact area calculations. Additionally, the area_shape keyword is described here below.

Area shape

If area_shape is set to circle [default] then the particle-particle ( $h_{ij}$ ) and particle-wall ( $h_{iw}$ ) heat transfer coefficients are calculated as indicated in the equations above. Note that circle is in general appropriate, but for certain reasons (backward compatibility) another setting might be chosen. If this keyword is set to square, then $\pi$ is replaced by 4 in the equation. If this keyword is set to legacy, then $\pi$ is omitted in the equation (replaced by 1). As the thermal conductivity is a calibration parameter this is not a big issue, as switching from the legacy representation (default up to Aspherix 5.3.x) to the physically more correct circle representation (default from Aspherix 5.4) can be achieved by dividing the current thermal conductivities by $\sqrt(\pi)$ .

In case enable_radiation is set to yes a radiation model based on local void fractions and particle distance is employed using view factors read from a file.

The model that was implemented is based on the work of Johnson et al. [2]. It defines the heat transferred from a particle $j$ to a particle $i$ as

$q_{r,ij} = \epsilon_j A_{ij} \sigma D_{ij} (T_j^4 - T_i^4),$

where $\epsilon$ is the emissivity of the particles which can be specified with the corresponding keyword. $A_{ij}$ is the average surface area of the two particles and $\sigma$ the Stefan Boltzmann constant. $T_i$ and $T_j$ are the absolute temperatures of the two praticles and $D_{ij}$ the view factor which depends on the average local volume fraction $(\alpha_i + \alpha_j)/2$ and the distance between the two particles centroids $r_{ij}$ .

The area $A_{ij}$ can also be set as constant for all contacts using the apparent_particle_radius keyword, which defines the area based on this radius value.

The total heat transferred to a particle is calculated as the sum of all particle-particle heat rates where the particle centroid distance is less than a certain value $r_{max}$ . This is either defined via the max_relative_radius keyword or taken from the largest radius listed in the view_factors_file. This file is a csv file with comma (or space) separated values. The first line must be a header which must define the first column as ‘r’ and all subsequent columns via their corresponding volume fraction. The rows of the data that will follow below specify the relative radius in the first row and then the view factor corresponding to that radius and the volume fraction. As very basic file could look like:

r, 0.2, 0.4
, 1e-1, 1
, 5e-2, 1e-1
, 1e-1, 1e-2

So for a local volume fraction of 0.2 and a particle centroid distance of 3 particle radii, the view factor would be 5e-2. The data is read in at the beginning of each simulation and then linearly interpolated into a lookup table. This table has $n*n$ entries, where $n$ is given by the interpolation_resolution (default value 20). If a value is outside the table the nearest value will be chosen.

Finally, the local volume fraction is computed as:

$\alpha_i = \frac{3}{4 \pi r_{max}^3} \sum_{j: r_{ij} \leq r_{max}} V_j,$

where $V_j$ is the particle volume.

External heat flux

If add_external_heatflux is set to yes an external heat flux can be imposed on particles directly, as described below.

The external application must be able to read the output created with the output_settings command. This command will write an empty file output_*_written.txt after it has completed its writing indicating that the external application can start reading the VTK data. The * indicates the time step that has just been written. It is expected that the external application deletes this file once it has completed its writing, otherwise a deadlock will occur. After this file deletion Aspherix will start reading the external heat flux present in the particle data.

The heatFluxExternal field that is part of the particle data will be read after being written by an external program and will be imposed directly onto the particle.

In case add_external_heatflux_from_grid is set to yes an external heat flux is distributed onto particles.

This model works in combination with the calculate spatial_average command and allows using an external application to write a heat flux onto this Eulerian grid that is then distributed to all the particles in a cell of this grid.

See add_external_heatflux details above for how to write the data.

The additional heat flux a particle will see is thus

$h_i = h_g/n_g,$

where $h_g$ is the heat flux of the grid cell in which particle $i$ is located and $n_g$ is the number of particles in this grid cell.

In case output_detailed_heatflux is set to yes additional per particle properties are included in the output which contain the heat flux coming from both conduction and radiation. The heatFluxParticleConduction and heatFluxParticleRadiation particle properties are created and output using the output_settings command which each contain the instanteaneous heat flux a particle experiences due to conduction and radiation, respectively.

Coarse-graining information:

Using coarsegraining in combination with this command should lead to statistically equivalent dynamics and system state.

Output information

This command computes a scalar which can be accessed by various output commands. This scalar is the total thermal energy of the particles in the group associated with this command. It can also be accessed using id_command-ID.total_thermal_energy.

The temperature, heat flux and heat source of an atom with index i can be obtained by id_Temp[i], id_heatFlux[i] and id_heatSource[i], respectively. This can be used to define variables, using the variable command, that can be accessed by output commands.

Additional information

The particle temperature and heat source is written to binary restart files so simulations can continue properly.

Restrictions

Warning

The following functionalities do not have full GPU support: enable_radiation, area_correction yes, modified_area_correction yes, store_contact_data yes and dot access for the total thermal energy.