As a simple example of the coding of user subroutine
VUMATHT, consider a material with isotropic conductivity and a
constant specific heat. The heat transfer equations are developed here, and the
example
VUMATHT is given. The problem can also be solved by directly
specifying thermal conductivity, specific heat, density, and internal heat
generation.
First, the heat transfer equations are outlined.
The basic energy balance is
where V is the volume of solid material with surface
area S,
is the density of the material,
is the material time rate of the internal thermal energy,
q is the heat flux per unit area of the body flowing into
the body, and r is the heat supplied externally into the
body per unit volume. q and r can
depend on the strains or displacements of the body.
In
Abaqus/Explicit
the heat transfer and mechanical solutions are obtained simultaneously by an
explicit coupling. Thus, the material time rate of the internal thermal energy
can be written as
A heat flux vector is defined such that
where is the unit outward
normal to the surface S. Introducing the above relation
into the energy balance equation and using the divergence theorem, the
following relation is obtained:
The corresponding weak form is given by
where
is the temperature gradient and
is an arbitrary variational field satisfying the essential boundary conditions.
In
Abaqus/Explicit
this nonlinear system is solved using the explicit central-difference time
integration rule with a lumped capacitance matrix as described in
Fully Coupled Thermal-Stress Analysis.
The thermal constitutive behavior for this example is now defined.
Isotropic conductivity and a constant specific heat for the material are
assumed. The heat conduction in the material is assumed to be governed by
Fourier's law.
The internal thermal energy per unit mass is defined as
and its variation with respect to temperature is
where c is the specific heat of the material and
is the value of absolute zero on the temperature scale being used.
Fourier's law for heat conduction is given as
where is the thermal
conductivity matrix.
An approximation to the stability limit for the forward-difference operator
in the thermal solution response is given by
where
is the smallest element dimension in the mesh and the parameters
and
represent the material's effective specific heat and effective thermal
conductivity, respectively. For a conservative estimate of the time increment
size of the thermal solution response, we can specify
where
are components of the thermal conductivity matrix .
No state variables are needed for this material, so the allocation of space
for them is not necessary.
A thermal user material definition can be used to read in the two constants
for our simple case, namely the specific heat, c, and the
coefficient of isotropic thermal conductivity, k, so that
subroutine vumatht (
C Read only (unmodifiable) variables -
* nblock, nElem, nIntPt, nLayer, nSectPt,
* ntgrad, nstatev, nfieldv, nprops,
* cmname, stepTime, totalTime, dt,
* coordMp, density, props,
* tempOld, fieldOld, stateOld, enerThermOld,
* tempNew, tempgradNew, fieldNew,
C Write only (modifiable) variables -
* stateNew, fluxNew, enerThermNew, dEnerThDTemp, condEff )
C
include 'vaba_param.inc'
C
dimension nElem(nblock)
C
dimension coordMp(nblock,*), density(nblock), props(nprops),
* tempOld(nblock), fieldOld(nblock, nfieldv),
* stateOld(nblock, nstatev), enerThermOld(nblock),
* tempNew(nblock), tempgradNew(nblock, ntgrad),
* fieldNew(nblock, nfieldv), stateNew(nblock, nstatev),
* fluxNew(nblock, ntgrad), enerThermNew(nblock),
* dEnerThDTemp(nblock,2), condEff(nblock)
C
character*80 cmname
*
parameter ( zero = 0.d0 )
parameter ( tempZero = 0.d0)
* properties
* 1. conductivity
* 2. specific heat
*
cond = props(1)
specHeat = props(2)
*
do km = 1, nblock
*
* update heat flux vector
do i = 1, ntgrad
fluxNew(km, i) = -cond*tempgradNew(km,i)
end do
condEff(km) = cond
*
* update internal thermal energy and its derivatives
if (totalTime .eq. zero) then
dEnerTh = specHeat*(tempNew(km) - tempZero)
else
dEnerTh = specHeat*(tempNew(km) - tempOld(km))
end if
enerThermNew(km) = enerThermOld(km) + dEnerTh
dEnerThDTemp(km,1) = specHeat
dEnerThDTemp(km,2) = zero
end do
*
return
end