Buildings.HeatTransfer.Conduction

Package with models for heat conduction

Information

This package provides component models to compute heat conduction.

Implementation

The package declares the constant nSupPCM, which is equal to the number of support points that are used to approximate the specific internal energy versus temperature relation. This approximation is used by Buildings.HeatTransfer.Conduction.SingleLayer to replace the piece-wise linear function by a cubic hermite spline, with linear extrapolation, in order to avoid state events during the simulation.

Extends from Modelica.Icons.VariantsPackage (Icon for package containing variants).

Package Content

Name Description
Buildings.HeatTransfer.Conduction.MultiLayer MultiLayer Model for heat conductance through a solid with multiple material layers
Buildings.HeatTransfer.Conduction.SingleLayer SingleLayer Model for single layer heat conductance
Buildings.HeatTransfer.Conduction.SingleLayerCylinder SingleLayerCylinder Heat conduction in a cylinder
nSupPCM=6 Number of support points to approximate u(T) releation, used only for phase change material
Buildings.HeatTransfer.Conduction.BaseClasses BaseClasses Package with base classes for Buildings.HeatTransfer.Conduction

Types and constants

  constant Integer nSupPCM = 6
    "Number of support points to approximate u(T) releation, used only for phase change material";

Buildings.HeatTransfer.Conduction.MultiLayer Buildings.HeatTransfer.Conduction.MultiLayer

Model for heat conductance through a solid with multiple material layers

Buildings.HeatTransfer.Conduction.MultiLayer

Information

This is a model of a heat conductor with multiple material layers and energy storage. The construction has at least one material layer, and each layer has at least one temperature node. The layers are modeled using an instance of Buildings.HeatTransfer.Conduction.SingleLayer. See this model for an explanation of the equations that are applied to each material layer.

Important parameters

The construction material is defined by a record of the package Buildings.HeatTransfer.Data.OpaqueConstructions. This record allows specifying materials that store energy, and material that are a thermal conductor only with no heat storage. To assign the material properties to this model, do the following:

  1. Create an instance of a record of Buildings.HeatTransfer.Data.OpaqueConstructions, for example by dragging the record into the schematic model editor.
  2. Make sure the instance has the attribute parameter, which may not be assigned automatically when you drop the model in a graphical editor. For example, an instanciation may look like
     parameter Data.OpaqueConstructions.Insulation100Concrete200 layers
       "Material layers of construction"
       annotation (Placement(transformation(extent={{-80,60},{-60,80}})));
    
  3. Assign the instance of the material to the instance of the heat transfer model as shown in Buildings.HeatTransfer.Examples.ConductorMultiLayer.

The parameters stateAtSurface_a and stateAtSurface_b determine whether there is a state variable at these surfaces, as described above. Note that if stateAtSurface_a = true, then there is temperature state on the surface a with prescribed value, as determined by the differential equation of the heat conduction. Hence, in this situation, it is not possible to connect a temperature boundary condition such as Buildings.HeatTransfer.Sources.FixedTemperature as this would yield to specifying the same temperature twice. To avoid this, either set stateAtSurface_a = false, or place a thermal resistance between the boundary condition and the surface of this model. The same applies for surface b. See the examples in Buildings.HeatTransfer.Examples.

Extends from Buildings.HeatTransfer.Conduction.BaseClasses.PartialConductor (Partial model for heat conductor), Buildings.HeatTransfer.Conduction.BaseClasses.PartialConstruction (Partial model for multi-layer constructions).

Parameters

TypeNameDefaultDescription
AreaA Heat transfer area [m2]
ThermalResistanceRsum(lay[i].R for i in 1:nLay)Thermal resistance of construction [K/W]
Genericlayersredeclare parameter Building...Construction definition from Data.OpaqueConstructions
Initialization
BooleansteadyStateInitialfalse=true initializes dT(0)/dt=0, false initializes T(0) at fixed temperature using T_a_start and T_b_start
TemperatureT_a_start293.15Initial temperature at port_a, used if steadyStateInitial = false [K]
TemperatureT_b_start293.15Initial temperature at port_b, used if steadyStateInitial = false [K]
Dynamics
BooleanstateAtSurface_atrue=true, a state will be at the surface a
BooleanstateAtSurface_btrue=true, a state will be at the surface b

Connectors

TypeNameDescription
HeatPort_aport_aHeat port at surface a
HeatPort_bport_bHeat port at surface b

Modelica definition

model MultiLayer "Model for heat conductance through a solid with multiple material layers" extends Buildings.HeatTransfer.Conduction.BaseClasses.PartialConductor( final R=sum(lay[i].R for i in 1:nLay)); Modelica.SIunits.Temperature T[sum(layers.nSta)]( each nominal = 300) "Temperature at the states"; Modelica.SIunits.HeatFlowRate Q_flow[sum(layers.nSta)+nLay] "Heat flow rate from state i to i+1"; extends Buildings.HeatTransfer.Conduction.BaseClasses.PartialConstruction; parameter Boolean stateAtSurface_a=true "=true, a state will be at the surface a"; parameter Boolean stateAtSurface_b=true "=true, a state will be at the surface b"; protected Buildings.HeatTransfer.Conduction.SingleLayer[nLay] lay( final nSta2={layers.nSta[i] for i in 1:nLay}, each final A=A, final stateAtSurface_a = {if i == 1 then stateAtSurface_a else false for i in 1:nLay}, final stateAtSurface_b = {if i == nLay then stateAtSurface_b else false for i in 1:nLay}, material = {layers.material[i] for i in 1:size(layers.material, 1)}, T_a_start = { T_b_start+(T_a_start-T_b_start) * 1/R * sum(layers.material[k].R for k in i:size(layers.material, 1)) for i in 1:size(layers.material, 1)}, T_b_start = { T_a_start+(T_b_start-T_a_start) * 1/R * sum(layers.material[k].R for k in 1:i) for i in 1:size(layers.material, 1)}, each steadyStateInitial = steadyStateInitial) "Material layer"; equation // This section assigns the temperatures and heat flow rates of the layer models to // an array that makes plotting the results easier. for i in 1:nLay loop for j in 1:layers.nSta[i] loop T[sum(layers.nSta[k] for k in 1:(i-1)) +j] = lay[i].T[j]; end for; for j in 1:layers.nSta[i]+1 loop Q_flow[sum(layers.nSta[k] for k in 1:i-1)+(i-1)+j] = lay[i].Q_flow[j]; end for; end for; connect(port_a, lay[1].port_a); for i in 1:nLay-1 loop connect(lay[i].port_b, lay[i+1].port_a); end for; connect(lay[nLay].port_b, port_b); end MultiLayer;

Buildings.HeatTransfer.Conduction.SingleLayer Buildings.HeatTransfer.Conduction.SingleLayer

Model for single layer heat conductance

Buildings.HeatTransfer.Conduction.SingleLayer

Information

This is a model of a heat conductor for a single layer of homogeneous material that computes transient or steady-state heat conduction.

Main equations

Transient heat conduction in materials without phase change

If the material is a record that extends Buildings.HeatTransfer.Data.Solids and its specific heat capacity (as defined by the record material.c) is non-zero, then this model computes transient heat conduction, i.e., it computes a numerical approximation to the solution of the heat equation

ρ c (∂ T(s,t) ⁄ ∂t) = k (∂² T(s,t) ⁄ ∂s²),

where ρ is the mass density, c is the specific heat capacity per unit mass, T is the temperature at location s and time t and k is the heat conductivity. At the locations s=0 and s=x, where x is the material thickness, the temperature and heat flow rate is equal to the temperature and heat flow rate of the heat ports.

Transient heat conduction in phase change materials

If the material is declared using a record of type Buildings.HeatTransfer.Data.SolidsPCM, the heat transfer in a phase change material is computed. The record Buildings.HeatTransfer.Data.SolidsPCM declares the solidus temperature TSol, the liquidus temperature TLiq and the latent heat of phase transformation LHea. For heat transfer with phase change, the specific internal energy u is the dependent variable, rather than the temperature. Therefore, the governing equation is

ρ (∂ u(s,t) ⁄ ∂t) = k (∂² T(s,t) ⁄ ∂s²).

The constitutive relation between specific internal energy u and temperature T is defined in Buildings.HeatTransfer.Conduction.BaseClasses.temperature_u by using cubic hermite spline interpolation with linear extrapolation.

Steady-state heat conduction

If material.c=0, or if the material extends Buildings.HeatTransfer.Data.Resistances, then steady-state heat conduction is computed. In this situation, the heat flow between its heat ports is

Q = A   k ⁄ x   (Ta-Tb),

where A is the cross sectional area, x is the layer thickness, Ta is the temperature at port a and Tb is the temperature at port b.

Spatial discretization

To spatially discretize the heat equation, the construction is divided into compartments (control volumes) with material.nSta ≥ 1 state variables. Each control volume has the same material properties. The state variables are connected to each other through thermal resistances. If stateAtSurface_a = true, a state is placed at the surface a, and similarly, if stateAtSurface_b = true, a state is placed at the surface b. Otherwise, these states are placed inside the material, away from the surface. Thus, to obtain the surface temperature, use port_a.T (or port_b.T) and not the variable T[1].

As an example, we assume a material with a length of x and a discretization with four state variables.

To build multi-layer constructions, use Buildings.HeatTransfer.Conduction.MultiLayer instead of this model.

Important parameters

The parameters stateAtSurface_a and stateAtSurface_b determine whether there is a state variable at these surfaces, as described above. Note that if stateAtSurface_a = true, then there is temperature state on the surface a with prescribed value, as determined by the differential equation of the heat conduction. Hence, in this situation, it is not possible to connect a temperature boundary condition such as Buildings.HeatTransfer.Sources.FixedTemperature as this would yield to specifying the same temperature twice. To avoid this, either set stateAtSurface_a = false, or place a thermal resistance between the boundary condition and the surface of this model. The same applies for surface b. See the examples in Buildings.HeatTransfer.Examples.

Extends from Buildings.HeatTransfer.Conduction.BaseClasses.PartialConductor (Partial model for heat conductor).

Parameters

TypeNameDefaultDescription
AreaA Heat transfer area [m2]
ThermalResistanceRif (material.R < Modelica.Co...Thermal resistance of construction [K/W]
Materialmaterialredeclare parameter Data.Bas...Material from Data.Solids, Data.SolidsPCM or Data.Resistances
Initialization
BooleansteadyStateInitialfalse=true initializes dT(0)/dt=0, false initializes T(0) at fixed temperature using T_a_start and T_b_start
TemperatureT_a_start293.15Initial temperature at port_a, used if steadyStateInitial = false [K]
TemperatureT_b_start293.15Initial temperature at port_b, used if steadyStateInitial = false [K]
Dynamics
BooleanstateAtSurface_atrue=true, a state will be at the surface a
BooleanstateAtSurface_btrue=true, a state will be at the surface b
Advanced
IntegernSta2material.nStaNumber of states in a material (do not overwrite, used to work around Dymola 2017 bug)

Connectors

TypeNameDescription
HeatPort_aport_aHeat port at surface a
HeatPort_bport_bHeat port at surface b

Modelica definition

model SingleLayer "Model for single layer heat conductance" extends Buildings.HeatTransfer.Conduction.BaseClasses.PartialConductor( final R=if (material.R < Modelica.Constants.eps) then material.x/material.k/A else material.R/A); // if material.R == 0, then the material specifies material.k, and this model specifies x // For resistances, material.k need not be specified, and hence we use material.R // The value T[:].start is used by the solver when finding initial states // that satisfy dT/dt=0, which requires solving a system of nonlinear equations // if the convection coefficient is a function of temperature. Modelica.SIunits.Temperature T[nSta](start= if stateAtSurface_a then cat(1, {T_a_start}, {(T_a_start + (T_b_start - T_a_start)*UA*sum(RNod[k] for k in 1:i-1)) for i in 2:nSta}) else {(T_a_start + (T_b_start - T_a_start)*UA*sum(RNod[k] for k in 1:i)) for i in 1:nSta}, each nominal=300) "Temperature at the states"; Modelica.SIunits.HeatFlowRate Q_flow[nSta+1](each start=0) "Heat flow rates to each state"; Modelica.SIunits.SpecificInternalEnergy u[nSta]( each start=2.7E5, each nominal=2.7E5) "Definition of specific internal energy"; parameter Boolean stateAtSurface_a=true "=true, a state will be at the surface a"; parameter Boolean stateAtSurface_b=true "=true, a state will be at the surface b"; replaceable parameter Data.BaseClasses.Material material "Material from Data.Solids, Data.SolidsPCM or Data.Resistances"; parameter Boolean steadyStateInitial=false "=true initializes dT(0)/dt=0, false initializes T(0) at fixed temperature using T_a_start and T_b_start"; parameter Modelica.SIunits.Temperature T_a_start=293.15 "Initial temperature at port_a, used if steadyStateInitial = false"; parameter Modelica.SIunits.Temperature T_b_start=293.15 "Initial temperature at port_b, used if steadyStateInitial = false"; parameter Integer nSta2=material.nSta "Number of states in a material (do not overwrite, used to work around Dymola 2017 bug)"; protected final parameter Integer nSta= max(nSta2, if stateAtSurface_a or stateAtSurface_b then 2 else 1) "Number of state variables"; final parameter Integer nR=nSta+1 "Number of thermal resistances"; parameter Modelica.SIunits.ThermalResistance RNod[nR]= if (stateAtSurface_a and stateAtSurface_b) then if (nSta==2) then {(if i==1 or i==nR then 0 else R/(nSta-1)) for i in 1:nR} else {(if i==1 or i==nR then 0 elseif i==2 or i==nR-1 then R/(2*(nSta-2)) else R/(nSta-2)) for i in 1:nR} elseif (stateAtSurface_a and (not stateAtSurface_b)) then {(if i==1 then 0 elseif i==2 or i==nR then R/(2*(nSta-1)) else R/(nSta-1)) for i in 1:nR} elseif (stateAtSurface_b and (not stateAtSurface_a)) then {(if i==nR then 0 elseif i==1 or i==nR-1 then R/(2*(nSta-1)) else R/(nSta-1)) for i in 1:nR} else {R/(if i==1 or i==nR then (2*nSta) else nSta) for i in 1:nR} "Thermal resistance"; parameter Modelica.SIunits.Mass m[nSta]= (A*material.x*material.d) * (if (stateAtSurface_a and stateAtSurface_b) then if (nSta==2) then {1/(2*(nSta-1)) for i in 1:nSta} elseif (nSta==3) then {1/(if i==1 or i==nSta then (2*(nSta-1)) else (nSta-1)) for i in 1:nSta} else {1/(if i==1 or i==nSta or i==2 or i==nSta-1 then (2*(nSta-2)) else (nSta-2)) for i in 1:nSta} elseif (stateAtSurface_a and (not stateAtSurface_b)) then {1/(if i==1 or i==2 then (2*(nSta-1)) else (nSta-1)) for i in 1:nSta} elseif (stateAtSurface_b and (not stateAtSurface_a)) then {1/(if i==nSta or i==nSta-1 then (2*(nSta-1)) else (nSta-1)) for i in 1:nSta} else {1/(nSta) for i in 1:nSta}) "Mass associated with the temperature state"; final parameter Real mInv[nSta]= if material.steadyState then zeros(nSta) else {1/m[i] for i in 1:nSta} "Inverse of the mass associated with the temperature state"; final parameter Modelica.SIunits.HeatCapacity C[nSta] = m*material.c "Heat capacity associated with the temperature state"; final parameter Real CInv[nSta]= if material.steadyState then zeros(nSta) else {1/C[i] for i in 1:nSta} "Inverse of heat capacity associated with the temperature state"; parameter Modelica.SIunits.SpecificInternalEnergy ud[Buildings.HeatTransfer.Conduction.nSupPCM]( each fixed=false) "Support points for derivatives (used for PCM)"; parameter Modelica.SIunits.Temperature Td[Buildings.HeatTransfer.Conduction.nSupPCM]( each fixed=false) "Support points for derivatives (used for PCM)"; parameter Real dT_du[Buildings.HeatTransfer.Conduction.nSupPCM]( each fixed=false, each unit="kg.K2/J") "Derivatives dT/du at the support points (used for PCM)"; initial equation // The initialization is only done for materials that store energy. if not material.steadyState then if steadyStateInitial then if material.phasechange then der(u) = zeros(nSta); else der(T) = zeros(nSta); end if; else if stateAtSurface_a then T[1] = T_a_start; for i in 2:nSta loop T[i] =T_a_start + (T_b_start - T_a_start)*UA*sum(RNod[k] for k in 1:i-1); end for; else // stateAtSurface_a == false for i in 1:nSta loop T[i] = T_a_start + (T_b_start - T_a_start)*UA*sum(RNod[k] for k in 1:i); end for; end if; end if; end if; if material.phasechange then (ud, Td, dT_du) = Buildings.HeatTransfer.Conduction.BaseClasses.der_temperature_u( c = material.c, TSol=material.TSol, TLiq=material.TLiq, LHea=material.LHea, ensureMonotonicity=material.ensureMonotonicity); else ud = zeros(Buildings.HeatTransfer.Conduction.nSupPCM); Td = zeros(Buildings.HeatTransfer.Conduction.nSupPCM); dT_du = zeros(Buildings.HeatTransfer.Conduction.nSupPCM); end if; equation port_a.Q_flow = +Q_flow[1]; port_b.Q_flow = -Q_flow[end]; port_a.T-T[1] = if stateAtSurface_a then 0 else Q_flow[1]*RNod[1]; T[nSta]-port_b.T = if stateAtSurface_b then 0 else Q_flow[end]*RNod[end]; for i in 1:nSta-1 loop // Q_flow[i+1] is heat flowing from (i) to (i+1) // because T[1] has Q_flow[1] and Q_flow[2] acting on it. T[i]-T[i+1] = Q_flow[i+1]*RNod[i+1]; end for; // Steady-state heat balance if material.steadyState then for i in 2:nSta+1 loop Q_flow[i] = port_a.Q_flow; end for; for i in 1:nSta loop if material.phasechange then // Phase change material T[i]=Buildings.HeatTransfer.Conduction.BaseClasses.temperature_u( ud=ud, Td=Td, dT_du=dT_du, u=u[i]); else // Regular material u[i]=0; // u is not required in this case end if; end for; else // Transient heat conduction if material.phasechange then // Phase change material for i in 1:nSta loop der(u[i]) = (Q_flow[i]-Q_flow[i+1])*mInv[i]; // Recalculation of temperature based on specific internal energy T[i]=Buildings.HeatTransfer.Conduction.BaseClasses.temperature_u( ud=ud, Td=Td, dT_du=dT_du, u=u[i]); end for; else // Regular material for i in 1:nSta loop der(T[i]) = (Q_flow[i]-Q_flow[i+1])*CInv[i]; end for; for i in 1:nSta loop u[i]=0; // u is not required in this case end for; end if; end if; end SingleLayer;

Buildings.HeatTransfer.Conduction.SingleLayerCylinder Buildings.HeatTransfer.Conduction.SingleLayerCylinder

Heat conduction in a cylinder

Buildings.HeatTransfer.Conduction.SingleLayerCylinder

Information

Model for radial heat transfer in a hollow cylinder.

If the heat capacity of the material is non-zero, then this model computes transient heat conduction, i.e., it computes a numerical approximation to the solution of the heat equation

ρ c ( ∂ T(r,t) ⁄ ∂t ) = k ( ∂² T(r,t) ⁄ ∂r² + 1 ⁄ r   ∂ T(r,t) ⁄ ∂r ),

where ρ is the mass density, c is the specific heat capacity per unit mass, T is the temperature at location r and time t and k is the heat conductivity. At the locations r=ra and r=rb, the temperature and heat flow rate are equal to the temperature and heat flow rate of the heat ports.

If the heat capacity of the material is set to zero, then steady-state heat flow is computed using

Q = 2 π k (Ta-Tb)⁄ ln(ra ⁄ rb),

where ra is the internal radius, rb is the external radius, Ta is the temperature at port a and Tb is the temperature at port b.

Implementation

To spatially discretize the heat equation, the construction is divided into compartments with material.nSta ≥ 1 state variables. The state variables are connected to each other through thermal conductors. There is also a thermal conductor between the surfaces and the outermost state variables. Thus, to obtain the surface temperature, use port_a.T (or port_b.T) and not the variable T[1].

Parameters

TypeNameDefaultDescription
Genericmaterialredeclare parameter Building...Material thermal properties
Heighth Height of the cylinder [m]
Radiusr_a Internal radius [m]
Radiusr_b External radius [m]
IntegernSta Number of state variables
RealgriFac2Grid factor for spacing
Initialization
TemperatureTInt_start293.15Initial temperature at port_a, used if steadyStateInitial = false [K]
TemperatureTExt_start293.15Initial temperature at port_b, used if steadyStateInitial = false [K]
BooleansteadyStateInitialfalsetrue initializes dT(0)/dt=0, false initializes T(0) at fixed temperature using T_a_start and T_b_start

Connectors

TypeNameDescription
HeatPort_aport_aHeat port at surface a
HeatPort_bport_bHeat port at surface b

Modelica definition

model SingleLayerCylinder "Heat conduction in a cylinder" replaceable parameter Buildings.HeatTransfer.Data.Soil.Generic material "Material thermal properties"; parameter Modelica.SIunits.Height h "Height of the cylinder"; parameter Modelica.SIunits.Radius r_a "Internal radius"; parameter Modelica.SIunits.Radius r_b "External radius"; parameter Integer nSta(min=1) "Number of state variables"; parameter Modelica.SIunits.Temperature TInt_start=293.15 "Initial temperature at port_a, used if steadyStateInitial = false"; parameter Modelica.SIunits.Temperature TExt_start=293.15 "Initial temperature at port_b, used if steadyStateInitial = false"; parameter Boolean steadyStateInitial=false "true initializes dT(0)/dt=0, false initializes T(0) at fixed temperature using T_a_start and T_b_start"; parameter Real griFac(min=1) = 2 "Grid factor for spacing"; Modelica.SIunits.TemperatureDifference dT "port_a.T - port_b.T"; Modelica.Thermal.HeatTransfer.Interfaces.HeatPort_a port_a "Heat port at surface a"; Modelica.Thermal.HeatTransfer.Interfaces.HeatPort_b port_b "Heat port at surface b"; Modelica.SIunits.Temperature T[nSta](start= {TInt_start+ (TExt_start-TInt_start)/Modelica.Math.log(r_b/r_a) *Modelica.Math.log((r_a+(r_b-r_a)/(nSta)*(i-0.5))/r_a) for i in 1:nSta}) "Temperature of the states"; Modelica.SIunits.HeatFlowRate Q_flow[nSta+1] "Heat flow rate from state i to i+1"; // Modelica.SIunits.TemperatureSlope der_T[nSta] // "Time derivative of temperature (= der(T))"; protected parameter Modelica.SIunits.Radius r[nSta+1](each fixed=false) "Radius to the boundary of the i-th domain"; parameter Modelica.SIunits.Radius rC[nSta](each fixed=false) "Radius to the center of the i-th domain"; final parameter Modelica.SIunits.SpecificHeatCapacity c= material.c "Specific heat capacity"; final parameter Modelica.SIunits.ThermalConductivity k= material.k "Thermal conductivity of the material"; final parameter Modelica.SIunits.Density d = material.d "Density of the material"; parameter Modelica.SIunits.ThermalConductance G[nSta+1](each fixed=false) "Heat conductance between the temperature nodes"; parameter Modelica.SIunits.HeatCapacity C[nSta](each fixed=false) "Heat capacity of each state"; initial equation assert(r_a < r_b, "Error: Model requires r_a < r_b."); assert(0 < r_a, "Error: Model requires 0 < r_a."); r[1] = r_a; for i in 2:nSta+1 loop r[i]= r[i-1] + ( r_b - r_a) * (1-griFac)/(1-griFac^(nSta)) * griFac^(i-2); end for; // Radii to the center of the domain for i in 1:nSta loop rC[i] = (r[i]+r[i+1])/2; end for; // Conductance between nodes (which are in the center of the domain) G[1] = 2*Modelica.Constants.pi*k*h/Modelica.Math.log(rC[1]/r_a); G[nSta+1] = 2*Modelica.Constants.pi*k*h/Modelica.Math.log(r_b/rC[nSta]); for i in 2:nSta loop G[i] = 2*Modelica.Constants.pi*k*h/Modelica.Math.log(rC[i]/rC[i-1]); end for; // Heat capacity of each segment for i in 1:nSta loop C[i] = (d*Modelica.Constants.pi*c*h*((r[i+1])^2-(r[i])^2)); end for; // The initialization is only done for materials that store energy. if not material.steadyState then if steadyStateInitial then der(T) = zeros(nSta); else for i in 1:nSta loop T[i]=TInt_start+ (TExt_start-TInt_start)/Modelica.Math.log(r_b/r_a) * Modelica.Math.log(rC[i]/r_a); end for; end if; end if; equation dT = port_a.T - port_b.T; port_a.Q_flow = +Q_flow[1]; port_b.Q_flow = -Q_flow[nSta+1]; Q_flow[1] = G[1] * (port_a.T-T[1]); Q_flow[nSta+1] = G[nSta+1] * (T[nSta]-port_b.T); for i in 2:nSta loop Q_flow[i]=G[i]*(T[i-1]-T[i]); // Q_flow[i] represents the heat flowing between two nodes end for; if material.steadyState then for i in 2:nSta+1 loop Q_flow[i]=Q_flow[1]; end for; else for i in 1:nSta loop der(T[i])= (Q_flow[i]-Q_flow[i+1])/C[i]; end for; end if; end SingleLayerCylinder;