Transient analysis

Initialization

Initial values are assigned to all the element temperatures, conductances, and heat inputs. The initial temperatures default to zero unless a file TEMPF, which contains calculated temperatures, is present, or initial temperatures have been specified, in which case these define the initial temperatures.

For multilayer elements, the layer temperatures are initialized to the temperature of the middle layer.

The temperatures of uninitialized Oppenheim elements are initialized at the average temperature of the elements the Oppenheim element is connected to.

If the element CG method is used, the boundary element temperatures are initialized to the weighted average temperatures of the elements they are connected to.

For more information on initialization of hydraulic elements, see Hydraulic network analysis.

Using the user-subroutine, you can add new conductances to the conductance matrix.

Time t at the start of the run is defined by the TST parameter.

If the element CG method is used, large negative conductances are modified according to the MODCOND values.

Conductive conductances for hydraulic elements are calculated from the CG to the boundary:

where:

G is the conductance from the CG to the boundary.
A_d is the duct cross-section.
L is the length of the duct.
K_therm is the thermal conductivity of the fluid.

Solving hydraulic network

Total pressures of hydraulic elements are calculated for the current time step. For more information, see Hydraulic network analysis.

Calculation of conductances

Convective and free convection conductances associated with hydraulic elements are updated as described in the Hydraulic network analysis section. One-way conductances for hydraulic networks are updated by multiplying the mass flow values flowing through DUCT and FLOWRES elements with the specific heat. Free convection conductances and not associated with hydraulic elements are calculated with:

where:

G_ij,t is the value of the conductance between elements i and j.
i,t refers to the property of element i at time t.
G_{F_ij} is the free convection conductance parameter.
T_i,t is the temperature of element i.
T_j,t is the temperature of element j.
EXP is the specified exponent.

Integration time step calculation

Next step is calculation of integration time step, dt, the time value at which the heat loads and temperature boundary conditions are evaluated, t_calc, and the smallest value of in the model, RC_MIN.

The integration time step dt is calculated, according to specified the integration time step parameter DT.

dt is checked to ensure that t + dt does not exceeded the time of the next printout interval, defined by the elapsed time between printouts, DTP, and the final time of the run ,TF. If t + dt is too large, dt is appropriately reduced.

The parameter t_calc is defined as the time value at which the heat loads and temperature boundary conditions are evaluated. For implicit methods when calculating the temperatures at time t + αdt, t_calc is t + αdt where α is the explicit-implicit weighting factor, between 0 and 1.

For Crank-Nicolson α = 0.5, for forward and exponential forward differencing α = 0 and for backward differencing α = 1.

When the initial values are to be printed, t_calc is set to the starting time of a transient run, TST.

Table Interpolation

Linear table interpolations are performed at t = t_calc. Heat inputs, sink temperatures, and table- dependent conductances are evaluated.

If the independent variable is greater than the largest value in the table, the last dependent variable is used, and a warning message is issued. If the independent variable is greater than the smallest value in the table, the first dependent variable is used. In other words, the first and last values of the table are extended beyond the table limits.

If electrical elements are present, the electrical network is solved for voltage and current. The power dissipation in the electrical network is computed and summed with other heat inputs specified for the electrical elements.

If series or interface conductances are present, the appropriate conductance values are modified.

If the default element CG formulation is used and restriction and the CAPDIST OFF parameter is not present, then the heat loads on solid elements are redistributed from the CG to the boundary elements. For tetrahedra, a uniform weighting factor is used. For wedges and hexahedra, the weighting factors are proportional to the areas of the boundary elements.

Thermostat-controlled heat inputs

If thermostats or PID controllers have been specified, the heat load to each heater element is checked as to whether it should be zero or not. For PID controllers, the heat input Q to the heaters is calculated using the following equations:

for

where:

Q_htr is the maximum heat input to the heaters.
T₂ is the controller’s set point.
T₇ is the controller’s gain.
T₈ is the integral constant.
T₉ is the derivative constant.
T₁₀ is the bias.

For thermostats, a heater element is turned on if the temperature T_i,t of the sensor element falls below the cut-in temperature of the thermostat. It is turned off, if the sensor element temperature is above the cut-off temperature. If the sensor element temperature is between the cut-in and cut-off temperatures, the heater's on/off status is not altered.

For proportional controllers, the heat input Q to the heaters follows the following relationship:

where:

T(N₁) is the sensor temperature.
Q_htr is the maximum heater power.
T₂ is the cut-in temperature of the heater.
T₃ is the cut-out temperature for the proportional controller.

Modification temperature-dependent orthotropic conductances

If orthotropic materials are present, all conductances associated with orthotropic elements are recalculated, provided the maximum temperature change of an orthotropic element has not exceeded 20 degrees.

Adding series conductances

If series or interface conductances are present the conductances the interface element is connected to are modified.

Using the user–written subroutine to modify heat inputs

Next, the optional user–written subroutine USER1 is entered. You can modify conductances, heat loads, etc. before calculating temperatures.

Initial printout at start time of transient run

If t_calc = TST, initial printout is performed. t_calc is then incremented to TST + αdt, and the program loops to table interpolation for heat inputs to recalculate boundary conditions.

Solving for temperatures

The temperatures at time t+dt are calculated:

Explicit techniques

For the exponential forward differencing technique:

where:

τ_i,t is the time constant of the element.
T_{EQ_i,t} is the equilibrium temperature the element would reach if all its neighboring element temperatures were constant.
G_{SUM_i,t} is the sum of conductances G_ij,t for element i at time t.
GT_{SUM_i,t} is the sum of the T_j,tG_ij,t for element i at time t.
G_SUMCi,t is the sum of non-radiative conductances for elementi at time t.
G_SUMRi,t is the sum of all radiative conductances elementi at time t.
GT_SUMCi,t is the sum of GR_ij's for all non-radiative conductances for elementi at time t.
GT_SUMRi,t is the sum of the T_j,t⁴GR_ij for all radiative conductances elementi at time t.

For the forward differencing technique:

Implicit Techniques

For the implicit techniques, all temperatures are calculated iteratively. The implicit techniques are unconditionally stable if the explicit–implicit weighting factor α > 0.5.

If the default conjugate gradient solver is used, T_{i,est, t+dt,n+1} is calculated using the technique described below, where T_{i,est, t+dt,n+1} is the estimated temperature of element i at iteration n+1 without damping calculated with the conjugate-gradient technique. All parameters used at time t+dt are calculated using temperature and conductance values at iteration n.

Next, damping is used to calculate the new set of temperatures:

where:

T_i,t+dt,n+1 is the temperature estimate of element i at iteration n+1.
DP_t is the transient damping parameter.

During iteration, only the radiative conductances are updated; all other conductances are evaluated at t. The maximum number of iterations defaults to 100 and may be overridden. If the maximum number of iterations is reached, a warning message is written to the verbose log file.

If the TRACE parameter is enabled, the highlights of the current iteration are printed to the screen as well as to the report log file.

Convergence is achieved when the maximum temperature difference between two subsequent iterations is less than the product TDIFS * DP_t, where TDIFS is the temperature iteration convergence criterion and its default value = 0.001.

Initial temperature estimation for the implicit techniques

To accelerate convergence, the first iteration is estimated from the past behavior of the equilibrium temperature of the element.

where dtp is the previous integration time step.

Phase change calculations

Next, the temperatures of phase change elements are appropriately adjusted to account for the energy absorbed by the latent heat of absorption. For the explicit techniques, this adjustment occurs at the completion of the first iteration only. For the implicit techniques, it occurs at the completion of each iteration.

The adjustment is as follows:

For each phase change element, a quality QUAL_i,t+dt is defined, which is the fraction of the element in its hotter phase. Initially, QUAL_i,t+dt is set to QUAL_i,t. If T_i,t is below the phase change temperature, QUAL_i,t = 0. If T_i,t is above the phase change temperature, QUAL_i,t = 1.
A check is performed to see if phase change effects should be considered. Phase change effects are not considered if both T_i,t and T_i,t+dt both lie on the same side of the phase change temperature.
If phase change is considered, the energy flowing into the element over the integration time step is partitioned into three components: energy that affects its temperature below the phase change temperature, energy that affects its quality at the phase change temperature, and energy that affects its temperature above the phase change temperature.

Using these three values, the temperature and the quality of the phase change element are appropriately adjusted.

Ablative elements are a subset of phase change elements. Once an ablative element enters its higher temperature phase, a determination is made whether it is considered to have burned up or merely changed properties.

If the change in material properties occurs, the thermal conductivity and capacitance are appropriately adjusted. If, instead, the element is considered to have burned up, then the capacitance is reduced by a factor of 1.E-6. For solid elements, thermal conductivity is increased by a factor of 1.E5; for shell elements, the in-plane thermal conductivity is set to zero, and the out-of-plane thermal conductivity is increased by a factor of 1.E5.

Incrementing time to t+dt

The time is incremented to t+dt.

Modifying temperatures

Once temperatures are calculated, USER1 is entered. This gives the option to alter some of the calculated temperatures.

Recovering substructures temperature and print

Next, a check is made to determine whether the printout is enabled.

If the printout is enabled and substructuring are performed, the temperatures of the eliminated elements are recovered.

Next, printout is performed. At every printout a set of results are written on files QNODEF, <simulation/model name>-<solution/analysis name>_report.log, PRESSF, TEMPF, tmgrslt.dat, and <simulation/model name>-<solution/analysis name>_verbose.log.

Run completion check

Next, a run completion check is performed. A transient run is considered completed if t = TF parameter, the final time of the run.

Using the user–written subroutine for multiple runs

Once the run is completed, USER1 is entered. At this time you can specify additional steady-state or transient runs by setting a flag in the array IR. For further details, please refer to Card 10 User-Written Subroutines.

If an additional run is specified, the program continues at solving hydraulic network of either the steady state or transient algorithms.

Multiple runs check

If an additional run is specified, the program re-enters at either Step 2 of the steady state or transient algorithm.

Calculation completed

If there are subsequent runs, a printout on files TEMPFn, PRESSFn, and QNODEFn is created, where n is the run number.