p_heat

Contents Functions

p_heat

Purpose

Formulation and material support for the heat equation.

Syntax

il = p_heat('default')

Description

This help starts by describing the main commands : p_heat Database and Dbval. Supported p_heat subtypes and their formats are then described. For theory see section 6.1.13.

Database,Dbval] ...

Element properties are described by the row of an element property matrix or a data structure with an .il field containing this row (see section 7.4). Element property functions such as p_solid support graphical editing of properties and a database of standard properties.

p_heat database

 il=p_heat('database');

Accepted commands for the database are

d3 Integ SubType : Integ integration rule for 3D volumes (default -3).
d2 Integ SubType : Integ integration rule for 2D volumes (default -3).

For fixed values, use p_heat('info').

Example of database property construction

  il=p_heat([100 fe_mat('p_heat','SI',1) 0 -3 3],...
             'dbval 101 d3 -3 2');

Heat equation element properties

Element properties are described by the row of an element property matrix or a data structure with an .il field containing this row. Element property functions such as p_beam support graphical editing of properties and a database of standard properties.

1 : Volume element for heat diffusion (dimension DIM)

  [ProId fe_mat('p_heat','SI',1) CoordM Integ DIM]

ProID	element property identification number
type	identifier obtained with fe_mat('p_beam','SI',1)
Integ	is rule number in integrules
DIM	is problem dimension 2 or 3 D

2 : Surface element for heat exchange (dimension DIM-1)

   [ProId fe_mat('p_heat','SI',2) CoordM Integ DIM]

ProID	element property identification number
type	identifier obtained with fe_mat('p_beam','SI',2)
Integ	is rule number in integrules
DIM	is problem dimension 2 or 3 D

SetFace

This command can be used to define a surface exchange and optionally associated load. Surface exchange elements add a stiffness term to the stiffness matrix related to the exchange coefficient Hf defined in corresponding material property. One then should add a load corresponding to the exchange with the source temperature at T₀ through a convection coefficient Hf which is Hf.T_0. If not defined, the exchange is done with source at temperature equal to 0.

model=p_heat('SetFace',model,SelElt,pl,il);

SelElt is a findelt command string to find faces that exchange heat (use 'SelFace' to select face of a given preselected element).
pl is the identifier of existing material property (MatId), or a vector defining an m_heat property.
il is the identifier of existing element property (ProId), or a vector defining an p_heat property.

Command option -load T can be used to defined associated load, for exchange with fluid at temperature T. Note that if you modify Hf in surface exchange material property you have to update the load.

Following example defines a simple cube that exchanges with thermal source at 55 deg on the bottom face.

 
model=femesh('TestHexa8'); % Build simple cube model
model.pl=m_heat('dbval 100 steel'); % define steel heat diffusion parameter
model.il=p_heat('dbval 111 d3 -3 1'); % volume heat diffusion (1)
model=p_heat('SetFace-load55',... % exchange at 55 deg
    model,...
    'SelFace & InNode{z==0}',... % on the bottom face
    100,... % keep same matid for exchange coef
    p_heat('dbval 1111 d3 -3 2')); % define 3d, integ-3, for surface exchange (2)
cf=feplot(model); fecom colordatapro
def=fe_simul('Static',model); % compute static thermal state
mean(def.def)

2D validation

Consider a bi-dimensional annular thick domain Ω with radii r_e=1 and r_i=0.5. The data are specified on the internal circle Γ_i and on the external circle Γ_e. The solid is made of homogeneous isotropic material, and its conductivity tensor thus reduces to a constant k. The steady state temperature distribution is then given by

− k Δθ(x,y) = f(x,y) in Ω. (9.13)

The solid is subject to the following boundary conditions

Γ_i (r=r_i) : Neumann condition

∂ θ

∂ n

(x,y) = g(x,y) (9.14)
Γ_e (r=r_e) : Dirichlet condition

θ(x,y)=θ_ext(x,y) (9.15)

In above expressions, f is an internal heat source, θ_ext an external temperature at r=r_e, and g a function. All the variables depend on the variable x and y.

The OpenFEM model for this example can be found in ofdemos('AnnularHeat').
Numerical application : assuming k=1, f=0, Hf=1e⁻¹⁰, θ_ext(x,y) = exp(x) cos(y) and g(x,y)= −exp(x)/r_i ( cos(y) x − sin(y) x ), the solution of the problem is given by

θ(x,y) = exp(x) cos(y)