This section provides an introduction to network concepts, techniques, and terminology. Network analog methods are discussed for thermal problems which involve conduction, convection, gray-body radiation, wavelength-dependent radiation, advection (mass flow), heat sources, or temperature sources. Recommendations for certain problem types are also presented.

Thermal analysts have for some time been trained to formulate heat transfer problems in terms of thermal networks which are analogous to resistor-capacitor electrical networks. The basic formula for current flow across an electrical resistor is:

where I is current, E is voltage, and R is electrical resistance.

In heat transfer, the analogous equation for heat flow across a thermal resistor is:

For transient current or heat flows, flow into an electrical or thermal capacitor must be included. For current flow I into an electrical capacitor, the equation is:

where C is the value of the electrical capacitor and dE/dt is the time derivative of the voltage across the capacitor. In heat transfer, the analogous equation for heat flow into a thermal capacitor is:

The value of C in this equation is taken from:

The parallels between electrical and heat transfer equations have led to the development of the thermal network concept. The remainder of this section defines more completely thermal network concepts, especially as they apply to the QTRAN thermal network program. QTRAN has taken the conventional thermal network concepts and has added to them vis-a-vis wavelength-dependent thermal radiation resistors, the three-node convective resistors, LaGrange cubic finite element resistors, and others. For those who are more comfortable with finite element concepts, it might be easier to think of a thermal resistor as a 1-D finite element. Another viewpoint is to simply think of the thermal resistors as computational molecules or numerical models, especially where normal thermal resistor concepts do not strictly apply (e.g., the three-node convective resistors and wavelength- dependent resistors).

This section briefly outlines the fundamentals of the electrical network analog approach for conduction networks. It includes multidimensional transient and steady-state network models and the associated mathematical relations which govern them.

Steady-state heat conduction networks in Cartesian coordinates will be the first network model discussed in this section due to their relative simplicity. Steady- state networks in general do not require thermal capacitors and are simple to derive from first principle considerations. Remember, however, that the conductive network is simply another method of arriving at a finite difference approximation to the problem. This implies that many of the numerical considerations which apply to more classical finite difference or finite element discretization schemes will also apply to thermal networks. There are two major reasons for resorting to a thermal network as the method of solving a problem: (1) the ease with which various combinations of boundary conditions may be handled as compared to the more classical finite difference or finite element schemes and (2) computational efficiency for large, strongly nonlinear problems. Network methods in general will allow more flexibility for the user (indeed, the finite element method may be thought of as a very powerful subset of the more general thermal network approach).

1-D conduction is the cornerstone for both 2-D and 3-D conduction. For starters, a Cartesian 1-D conductive resistor is defined as follows:

and the heat conducted from node 1 to node 2 is computed from the following expression:

Transient conduction involves the added complications of capacitors, where a single capacitor value is defined as follows:

If the node in question lies at a material interface, it may be convenient to assign one capacitor to the node for each material at the interface. The heat equation being solved for transient conduction at a given node [i] is then seen to be:

which states the sum of all the capacitances multiplied by the time derivative of node [i] is equal to all of the heat flows across all resistances into the node.

While the maximum stable time step in the system does not affect QTRAN’s stability, it may have an effect on the execution time. In general, the smaller the stable time step the more time a given problem may take to execute. A network with tiny stable time steps tends to generate stiff systems of equations which are more difficult to solve than non-stiff systems. However, a network whose nodes have predominantly small stable time steps will also give better resolution in the time domain to fast transients, especially at very early times in the simulation. A rule of thumb for the mesh spacing (which directly affects the explicit stable time step) is to compute the approximate mesh spacing from the following formula:

This is an extremely approximate formula and you can typically violate it by a factor of several. If violated by orders of magnitude, however, the results may be in serious error. The formula given can be arrived at from several approaches. One approach is to look at the most significant term of a 1-D conduction solution which is an infinite series. Another approach is to look at response times and decay rates of exponentially damped dynamic systems. Try it out and see what results are obtained. Review past analyses that have been performed and compare the mesh spacing to this criteria.

Although this mesh spacing criteria is important in regions undergoing rapid change, this criteria can typically be grossly violated in regions where the transient is very mild.

Cylindrical, spherical, polar, triangular, skew, or other networks all use the same underlying principles as Cartesian meshes for both steady-state and transient runs. However, there are special formulas for the conductive resistance in these coordinate systems that must be used to account for the differences in effective cross-sectional areas. For hollow cylinders and a 360-degree circumference, the effective conductive resistance is as follows:

And, for spheres:

Heat conduction with phase change introduces an extra term into the network capacitors. Specifically, the capacitors must now be able to properly account for latent heat effects at a given temperature as well as the normal specific heat effects. Furthermore, to be perfectly general a code must be able to account for multiple phase change temperatures and latent heat effects within each capacitor. In addition, simultaneous phase changes at multiple nodes with the transition direction going in either direction (e.g., either melting or freezing) must be allowed.

To be mathematically ideal, at least some phase changes are supposed to occur at a single temperature, not over a temperature band. This causes the internal energy of a material undergoing a phase transition to cease being a function of temperature. This can cause some rather significant numerical headaches, and the introduction of a narrow phase transition temperature range (say, 0.1 to 1.0 degrees in width) will not usually significantly affect the answers to most problems. The introduction of this artificial phase transition range has the benefit of causing the internal energy to be a function of temperature at all points, and this is a great advantage from a computational standpoint.

The advantage of smearing this phase change over a very narrow band is not, however, enough of an advantage to allow even QTRAN’s SNPSOR algorithm to converge if it relied only upon Newton’s Second Order Method. The “S” shape that the phase transition introduces into the nodal heat flow rate curve is a classic way to kill Newton’s methods. QTRAN, however, detects when the old time step’s temperature and the new time step’s temperature straddle or intrude into a phase transition region. At these times and for only those nodes involved, QTRAN’s SNPSOR algorithm reverts to a bisection algorithm technique which is immune to “S” curve effects. Thus, the SNPSOR algorithm remains convergent even for implicitly solved phase transition problems.

It is also possible to solve some problems involving phase transitions over a large temperature range by putting a bump into the specific heat curve and then being careful not to allow large temperature changes in a single time step. This has the advantage of sometimes being much faster than the built-in phase change algorithm, but it also carries some risk with regard to accuracy.

Phase transition problems tend to consume considerably more CPU time than normal problems. It is, therefore, an advantage to specify that a problem does not involve phase changes if it is known a-priori that this is the case. This is done by specifying 0 for all capacitor PHIDs. When using PATQ to translate a Patran neutral file, specify 0 for the PHID on all MID templates.

Convection networks can be a source of great anxiety. They tend to embody an empirical convection coefficient  h  in the convective resistors, and the formulas used to compute this h value are usually highly complex. Transition regimes also frequently exist, where the empirical correlation used for a particular configuration may be a function of Reynolds, Grashoff, Rayleigh, Prandtl, Graetz, and add infinitum number ranges.

For example, the pipe flow configuration in QTRAN uses five different correlations with seven factorial transition regimes based on everything from Reynolds number ranges to viscosity ratio ranges. Also, since the correlations are not in general continuous through the transition regimes, interpolation must be used when in a transition regime to maintain continuity and avoid numerical convergence problems. Further complications arise when mass flow through tubes or packed beds occur, as three temperatures are then needed. The following sections expand on these concepts.

The majority of the QTRAN convection configurations are in keeping with the conventional thermal resistor concept in that there are two nodes associated with each resistor. The conventional convective resistance is given by the following expression:

R is the thermal resistance, h is the convection coefficient, and Area is the surface area from which heat is being convected. The heat flow across the resistor is then given by:

The tricky detail of this operation, of course, is to compute an h value for the resistance expression. In general, empirical correlations must be relied upon. These correlations are usually somewhat tedious to evaluate, since they themselves are in general a function of temperature and/or temperature difference. Furthermore, each correlation is usually only applicable to a range of dimensionless parameters, and hence these parameters must be evaluated to determine which correlation to use for a given configuration. The result is a very painful (if done by hand) iterative solution for even a relatively simple convection problem.

QTRAN alleviates the pain somewhat by offering a convection library from which the user may select any one of 37 configurations which are in turn supported by 61 correlations (several generic correlations also exist). QTRAN will automatically and dynamically select the appropriate correlation for both the configuration type and the dimensionless parameter range for the resistor. A catalogue of available configurations is contained in Convection Library (Ch. 9) along with a catalogue of supporting correlations. All correlations are fully documented and referenced such that the user may easily locate the references from which the correlations were taken.

As stated previously, the majority of the QTRAN convection configurations are in keeping with the conventional thermal resistor concept in that two temperature nodes are associated with each resistor. However, certain configurations (i.e., configurations CFIG = 1, 2, 21, and 25) require the association of a third temperature node with the convective resistor. Configurations CFIG = 1, 2, and 21 are pipe flow configurations and CFIG = 25 is a flow through porous media configuration. For all four configurations, three temperature values are needed. The three temperatures needed are:

It is at this point that the computational hydrodynamics concept of “upwind differencing” must be introduced. Upwind differencing is a numerical technique that is used to maintain stability for a numerical solution which involves advection. Simply stated, upwind differencing assumes that information (in this case, heat) is carried only from upstream to downstream. Information (or heat) is never carried from downstream to upstream. The impact of upwind differencing on the three-node convective resistors is that the upstream node NEVER gains or loses heat through the convective resistor so long as the flow velocity associated with the resistor remains positive.

Thus, it is seen that all heat transfer through the three-node convective resistors occurs between the wall (or media) node and the downstream node. The upstream node is used solely as a reference temperature node for LMTD calculations (if appropriate) so long as the flow velocity is positive. The LMTD is used for those flow configurations for which it is appropriate, and a differential between the wall temperature and the average bulk fluid temperature is used for other configurations.

From this discussion, it is obvious to the casual observer that the three-node convective resistors do not account for advective heat transfer (i.e., the heat carried along on the mass flow). In order to complete a thermal model involving a three-node convective resistor, it is frequently necessary (depending on the boundary conditions of the problem) to use an advection resistor to properly account for the heat balance.

Another point that may not be quite so obvious is that the three-node convective resistors may not be directly used to compute entrance temperatures given that the wall (or media) and exit temperatures are known. Remember that the heat transfer to the upstream node is not a function of a the three-node convective resistor since the upstream node is used only as a reference temperature. The upstream node temperature must be taken from boundary conditions or else calculated from other thermal resistor, capacitor, and heat source influences.

Gray body radiation networks also differ slightly from the normal thermal resistor concept in that the potential across the radiative resistor is evaluated as σ * T 4, where σ is the Stefan-Boltzmann constant. QTRAN automatically uses the σ * T 4 potential when evaluating gray-body resistors, and hence the difference between radiative and conductive resistors is essentially transparent to the user. Note that the σ * T 4 potential is also applicable to non-surface radiosity nodes such as are encountered in gray networks.

The formula for radiative heat transfer used in thermal networks in QTRAN is as follows:

The formula for R is dependent upon whether R is a surface resistor, a view factor resistor with or without a participating media, and so on. The actual formulation of the specific resistor types is defined after Table 8‑1 and Table 8‑2 for gray body and wavelength dependent radiation respectively.

For more information on radiative networks, consult any radiative heat transfer book that describes the network approach. See Ref. 6. in Appendix A.

The gap radiation creates radiation coupling without the radiosity network or performing form factor calculations. Some assumption inherent in the formulation are that the area radiating is large relative to the distance between the radiating surfaces. A form factor may be included, but if it is significantly less than one, consideration should be given to using the radiosity network and calculating the appropriate form factors. Type 5 radiation resistors are used thus the emissivities must be constant. Type 5 resistor are constant and are treated internally as conductances. The last major assumption is that the two surfaces radiating to each other are of equal value. Since a form factor is included, small variation can be compensated for by adjusting emissivity or form factor values. The radiation conductance is determined with (7‑16)

Wavelength-dependent thermal radiation networks are a rather significant extension of the gray-body radiation network theory. In essence, the normal radiosity network is divided up into discrete frequency bands such that the emissivities, transmissivities, and absorptivities of the radiating materials can be assumed to be gray within each frequency band. For each frequency band that is used, a distinct but frequently parallel radiosity network is generated between the radiating surfaces and participating media. Each frequency band network is typically coupled to any of the other bands only at the surfaces of the radiating materials or at participating media nodes. This treatment allows reflections from dissimilar surfaces to be handled correctly, and is directly analogous to the multigroup treatment used in certain high-energy physics codes. This approach does not assume an average emissivity over the entire spectrum (e.g., the Air Force Weapons Lab’s “TRAP” code), but instead inherently allows emissivities, absorptivities, and transmissivities to be functions of frequency band, temperature, and/or time.

Advection occurs whenever some quantity is carried along on a mass flow stream. For example, if hot fluid enters a container, heat (energy) is carried with the fluid into the container. Similarly, if a cold fluid enters a hot container, cold (or negative energy) is carried with the fluid into the container. This transfer of positive or negative energy by mass flow must be accounted for when doing thermal simulations if such mass flows exist. The mathematical equation for the heat added to node [i] by such a mass flow stream is as follows:

 

One facet of numerical stability should be mentioned at this point because of its influence on advective resistors. It is necessary to limit the transport of energy by advective resistors to one direction only for the purposes of the QTRAN program. Specifically, energy is carried by advective resistors only from the upstream node to the downstream node (which node is upstream or downstream is determined by the sign of the MDOT mass flow term). This scheme is immediately seen to be equivalent to one-sided or “upwind” differencing which is a common practice among numerical fluid flow analysts, and is commonly used because of its stability enhancing characteristics.

By introducing upwind differencing, we can guarantee that the system of equations has a stable solution. This does not mean, however, that the SNPSOR algorithm in QTRAN can always solve the system of equations. This is due to the point iterative nature of the SNPSOR algorithm. As long as the system of equations is diagonally dominant, point iterative schemes tend to work very well. However, the introduction of advection tends to weaken this diagonal dominance. As the diagonal dominance weakens, point iterative methods have more and more difficulty converging, until they eventually fail. One way of aiding convergence (at the expense of speed) is to under-relax. With a small enough relaxation parameter, a solution will usually converge. However, be aware that CPU times increase rapidly as the relaxation parameter approaches 0.0.

The user has the option of applying the relaxation parameters to the entire system of equations, by node group types, or on a node by node basis. To under relax, apply the under relaxation only to the advection nodes by using a relaxation multiplier less than one. Apply this to the advection nodes. This will direct the under relaxation only to those nodes that need to be under relaxed and still give some flexibility, by allowing the solution module to search for an optimum relaxation parameter.

In virtually all thermal problems either temperature sources or heat sources exist. A temperature source is defined to be simply an externally controlled boundary temperature. A heat source is defined to be anything that contributes heat to the network, not including the resistors and capacitors.

Microfunctions are a unique QTRAN feature. Frequently, there is a need to apply very complicated temperature or heat sources that may consist of tabular data, sines or cosines, exponentials, constants, or terms which may or may not be multiplied together or divided by one another. QTRAN lets each discrete component of a complicated temperature or heat source be defined as a “microfunction”. Once all of these component microfunctions have been defined, then a library is available from which to build the more complicated functions (macrofunctions) that are needed for temperature or heat sources. Many macrofunctions are allowed to reference the same microfunction. Thus, for a very complicated tabular data microfunction, it is only necessary to build it once and then reference it as many times as needed. Other microfunctions can then be used to multiply or scale this complicated microfunction. Currently, the argument of a microfunction can be either time, temperature, a temperature difference, a temperature average, or a σ * (T[1] 4 - T[2] 4) radiation potential. If the radiation potential is specified, QTRAN will convert the T[1] and T[2] values to degrees absolute before computing the radiation potential value.

As stated above, a macrofunction consists of one or more microfunctions that are either multiplied or divided by one another. Quite complicated macrofunctions can be expressed in this manner. Note that, when building temperature or heat sources, it is quite allowable to assign more than one macrofunction to the same node. Multiple macrofunctions assigned to the same node are simply added together to compute the final source value.

This section briefly outlines the fundamentals of fluid flow networks and the associated equations.

Flow networks are incorporated into the Patran Thermal system to analyze the complex fluid flow system coupled with the heat transfer problem. The analysis is focused on average flow parameters, like mass flow rate, and pressure. Users requiring a detailed analysis of the flow field have to resort to CFD solvers.

The pressure drop in a fluid flowing in a duct from point I to J can be expressed as:

casting equation 1 in terms of mass flow rate.

linearizing the above equation.

or:

adding the head losses in flow with the loss coefficient, K.

Adding gravity and pump/turbine heads, the basic equation in matrix form can be written as:

The flow resistance is a function of pressure; therefore, the problem is nonlinear and an iterative procedure is required for solution. The resistances are computed, assembled and the matrix equations are solved for nodal pressures.