EnergyExamples

Introduction

This Mathematica Notebook contains examples of the EnergyWorker package being used to solve problems in heat-transfer and energy engineering, as applied to common structures and materials.

Typical problems with which EnergyWorker can help are:

* calculation of heatsink sizes
* assessing optimal pipe-lagging thickness
* predicting air-conditioning load
* cost-saving from insulation installation
* solar-gain
* degree-day calculations
* temperature distribution in structures

The primary audience for the package includes architects, building designers, energy consultants, engineers, teachers and students who are involved in energy studies, in thermal calculations for small engineering projects, or in scenario testing for energy-saving strategies or designs.

Review of functions

Temperature-at-node in structures: TemperatureSolver[]

Black-body functions: BlackBodySpectral[], BlackBodyTotal[]

Radiative transfer: NetRadiativeTransfer[], ViewFactor[]

Convective transfer for free and forced convection around plates and 
cylinders: ThermalConvectiveResistance[]

Conductive transfer through slabs, hollow cylinders and spheres:
	ThermalConductiveResistance[], SlabHeatFlow[], BiotNumber[],	 
	HollowCylinderHeatFlow[], and UValue[]

Solar radiation: SunInSky[], DirectInsolation[], Airmass[], and ExtinctionLoss[]

Conversions between *C, *F, and K: CToF[], CToK[], FToC[], FToK[], 
KToC[], and	KToF[]
	
Material and fluid properties: ThermalExpansionCoefficient[], ThermalConductivity[], ThermalDiffusivity[], Density[],  ThermalCapacity[], KinematicViscosity[] (defined for more than 80 common materials)

Surface radiative properties:   Emissivity[], SolarAbsorptivity[]  (defined for more than 40 common surfaces)

Examples

unit conversion

EnergyWorker's has built-in constants that you can use to convert between different systems of units. The constants behave as multiplicative qualifiers when converting into the SI (MKS) unit system. For example, to convert one Calorie to Joules, you enter

1 Calorie

To convert from Joules to the chosen unit, you divide by the unit name:

8.36/Calorie

Here are some other examples.

1 TonneCoalEquivalent

3 Therm

1. Erg

Conversions into MKS for length, time, etc., are calculated in a similar manner. Note that, where possible, Mathematica keeps numbers in their most exact form; floating-point numbers are only used where a floating-point number already exists somewhere in the calculation, otherwise integers or rational numbers are used. Distances are converted into meters, clock time into seconds, and calendar time into days.

3 cm

3.0 cm

Mile

Foot

Time units are converted into seconds. For example, how many seconds are there in a day, a year, and a minute?

Day

Year

Minute

Calendar information is calculated in days: the mid-point of June is approximately day 166 in the year.

June

temperature scale conversion

All the functions in EnergyWorker use Kelvin for the temperature scale. Conversion functions are provided to convert between Centigrade, Fahrenheit, and Kelvin.

CToK[0]

KToC[0]

FToC[32]

CToF[0]

FToK[212]

U-Values and thermal conductive resistance

The U-Value of a material gives the rate at which heat will pass through a slab or the material, for a given temperature difference between two slab surfaces (e.g. the inner and outer sides of a wall). For example, how many Watts will flow through a 6 mm thick, 1 m2 pane of window glass,  for every 1*C temperature difference between the inside and outside surfaces.

UValue[6 mm,"GlassWindow"]

If EnergyWorker does not have a predefined numerical value for the thermal conductivity of the material for which you wish to calculate the U-Value then an expression is returned which leaves the unknown thermal conductivity unevaluated.

UValue[6 mm,"CrownFlint"]

Alternatively, if you know the numeric value for the material's thermal conductivity (in W/m/K), then you can use that instead of the material's name.

UValue[6 mm, 1.4]

The U-Value of a material is closely linked with the thermal conductive resistance, which gives the reciprocal of the heat-flow through a structure per 1*C temperature differential, rather than the heat-flow per square meter. For example, to calculate the thermal conductive resistance of a 2m high, 1m wide, 6mm thick glass window pane:

ThermalConductiveResistance["Slab",
               2,1,6 mm,"GlassWindow"]

The heat flow for a 15*C temperature difference would be:

15/%

(Note that this ignores the effect of the convective resistance at both of the air-glass interfaces.)

thermal convective resistance

The above example assumes that the glass surfaces are 15*C different in temperature. For a normal domestic window with air on both sides this is unlikely to be true because there is a much higher thermal resistance due to convection on either side of the window.

EnergyWorker provides a function for calculating the thermal convective resistance for a number of geometrical situations. Just like Mathematica's built-in functions, you can obtain help on functions within EnergyWorker either by using the function browser in versions of Mathematica that have a Notebook interface, or by using '?'.

?ThermalConvectiveResistance

So, to use ThermalConvectiveResistance, you have to supply four arguments. The first tells EnergyWorker what kind of geometry is being used - horizontal plates have a different convective resistance to vertical plates, for instance. The geometry is entered as a string, in quotes. The second is a two-long list containing the surface and fluid temperatures (in Kelvin). The third and fourth arguments are the characteristic length for the surface, and its area (both in meters.) 

For the above-mentioned glass window pane, the free convective resistance on a cool day is much higher than the conductive resistance of the glass itself.

ThermalConvectiveResistance["VFlatPlate",
                    {CToK[5],CToK[20]},2,2]

So the heat flow (in W) for a 15*C temperature difference is

15/%

This example used free convection: the only fluid motion present is that created by the convection itself. You can also calculate convective resistances for surfaces in a forced flow (perhaps created by a fan) by specifying the flow speed (in m/s) as an option.

For example, the same temperature differential for a flat plate (with its heated side uppermost) in an airflow of 1.5m/s would result in a much higher heatflow:

15/ThermalConvectiveResistance["HFlatPlateUp",
                               {CToK[20],CToK[5]},2,2,
                               Flow->1.5]

EnergyWorker uses dry Air at standard atmospheric pressure as the default fluid. In this version of the package, Water (saturated, at atmospheric pressure) can be specified as an option:

15/ThermalConvectiveResistance["HFlatPlateUp",
                               {CToK[20],CToK[5]},2,2,
                               Flow->1.5,Fluid->"Water"]

Returning to the U-Value problem, considering the effect of convection on both sides of the  glass, the U-Value resulting from the convective resistance is approximately

1/(2 ThermalConvectiveResistance["VFlatPlate",{CToK[5],CToK[20]},1,1])

heatsink example

A small electronics project is being build into an anodized aluminum box whose upper surface is to act as the heatsink. Assuming there will be no forced convection, the box lid is 5 cm by 10 cm, the consumption is 25 W, and the surroundings are at 20*C, what will be the temperature of the box lid?

In equilibrium, the heatflow by free convection and radiation will equal the consumption (25 W). One way to tackle this problem is to define a function, flow, that returns the sum of the convective and radiative flows when the surface is at a given temperature, ts. The Mathematica function FindRoot can then be used to find at what surface temperature the flow equals 25W.

Clear[flow];
flow[ts_]:=((ts-293)/ThermalConvectiveResistance["HFlatPlateUp",
                         {ts,293},7.5 cm,5 cm 10 cm]  +
             NetRadiativeHeatFlow[{ts,293}, 5 cm 10 cm,"AluminumAnodized"])

FindRoot[flow[CToK[ts]]==25,{ts,200,210}]

As another exampl of convective loss, what will be maximum dissipation from a small electronic component (with a 2 cm square vertical surface) if its maximum surface temperature is 100*C and it only losses heat by convection?

Clear[flow];
flow[ts_]:=(ts-293)/ThermalConvectiveResistance["VFlatPlate",
                         {ts,293},1.4 cm,2 cm 2 cm]

flow[373]

Hence, the dissipation can be just over 400mW if the device is to remain within limits.

insolation

EnergyWorker contains functions for evaluating the amount of solar radiation (insolation) incident on surfaces at given times and locations. The amount of insolation on a surface depends on the direction of both the Sun and the surface. The SunInSky function allows you to calculate the approximate altitude and azimuth (in radians) of the Sun - it takes three arguments: the time-in-day (seconds), the day-in-year (days), and the longitude and latitude of site (radians, in a two-long list):

sunDirection=SunInSky[Noon,June,{0,49 Degree}]

Once the Sun's direction is known, the insolation (W/m2) on a zenith-facing Asphalt surface can be found using the DirectInsolation function:

in=DirectInsolation[ZenithFacing,
                    sunDirection]*
                      SolarAbsorptivity["AsphaltPavement"]

Note that any direction of surface can be used. ZenithFacing is just one of the surface directions defined and named within EnergyWorker; arbitrary surfaces are defined using a two-long list that gives the surface normal vector's altitude and azimuth.

black-body functions

The total irradiance of a black body at 5762 K is just over 62 MW/m2.

BlackBodyTotal[5762]

As well as the total energy emitted by a blackbody, you an use BlackBodySpectral to calculate the irrdiance at a given wavelength. For example, a black body at 609 K emits 995 W/m2/micron at a wavelength of 4 microns.

BlackBodySpectral[609,4 um]

You can use these functions to perform related calculations. For example, what fraction of radiation emitted by a black body at 2600 K falls between 400 and 700 nm (i.e. in the visible).

NIntegrate[ 10^6 BlackBodySpectral[2600,w] ,{w, 0.4 um ,0.7 um }]/
BlackBodyTotal[2600]

temperatures in structures

With EnergyWorker's TemperatureSolver function, you can calculate the temperature at nodes throughout structures that contain temperature sources, heat flows, and thermal resistances. For example,  using an elemental model, you can calculate the temperature distribution within a concrete wall, 10 cm thick, which has one face maintained at a temperature of 310 K, the other face at 300 K, and, from electric heating wires embedded in the concrete, an internal heat generation of 10 kW per cubic meter.

One way to tackle this problem is to divide the wall into four 1 m square slabs, each 2.5 cm in depth. For each slab, the thermal resistance and the generated heat-flow can be calculated.

r1=ThermalConductiveResistance["Slab",
                                1,
                                1, 
                                0.025,
                                "Concrete"];

flow1=1 1 0.025 (10^4);

The next step is to generate the network that describes the thermal components. The network has seven nodes: node 1 is the zero Kelvin reference node, nodes 2 and 7 are the surfaces of the wall, and nodes 3, 4, 5, and 6 are internal nodes, each at the centre of a slab 1 m by 1 m by 2.5 cm. The surface temperatures are defined using Temperature components that (like the other thermal components) have the nodes between which they are connected as a two-long list for their first argument and the value for the second argument. The ThResistance components between  nodes 2 and 3 and between 6 and 7 are half the value of the other internode resistances because there is only half a slab's thickness between those nodes and the next (surface) node. The thermal flows are connected to each of the internal nodes, and all flow from the zero Kelvin node.

You need to declare the network, to which you can assign a name, inside a ThermalNet function:

wall=ThermalNet[Temperature[{1,2},310],
                ThResistance[{2,3},r1/2],
                ThFlow[{1,3},flow1],
                ThResistance[{3,4},r1],
                ThFlow[{1,4},flow1],
                ThResistance[{4,5},r1],
                ThFlow[{1,5},flow1],
                ThResistance[{5,6},r1],
                ThFlow[{1,6},flow1],
                ThResistance[{6,7},r1/2],
                Temperature[{1,7},300.]
               ];

Once wall is defined, you can call TemperatureSolver to return a list of node temperatures; node 1 is zero Kelvin, node 4 is nearly 319 Kelvin.

TemperatureSolver[wall]

These node temperatures can be converted to Centrigrade by  mapping the KToC function over each element of the temperature list. You can then plot the temperature distribution through the wall (we've removed node 1 by using Drop but note that the node numbers are reassigned), as a function of node number.

KToC[#]& /@ %

tempThruWall=Drop[%,1]

ListPlot[tempThruWall,PlotRange->{20,50}]

Product information

Further information about the package EnergyWorker is available from

Goth, Goth & Chandleri Ltd.,
The Old Sweet Factory,
Gardner Street,
Herstmonceux,
East Sussex  BN27 4LB,
United Kingdom

tel: +44-1323-832-683, fax:+44-1323-832-588
email: CompuServe 100070,3267

or, in the USA and Canada,

Macallan Consulting,
1583 Pinewood Way,
Milpitas,
CA 95035,
USA

(408)-262-3575