Thermo-fluid-dynamics, also known as thermal-fluid sciences, is an interdisciplinary field that combines principles from thermodynamics, fluid mechanics, and heat transfer to study and analyze the behavior of thermal systems. These systems involve the transfer of energy in the form of heat and the movement of fluids such as air and water so it could be interesting to analyze also small thermal systems from a user perspective.
Here, a small domestic thermal system application is analyzed with some perspective of thermodynamics rising up during the personalization of the irrespective control system to optimize comfort reducing energy consumption.
Furthermore, the following is mainly a personal reflection on these systems however I write it down with the intent of dedicating it to the professor of the Thermal Physics course who was more worried that I remembered the Guoj-Stodola theorem as well as his publications instead of teaching something really useful and applicable.
Thermal systems encompass a wide range of applications, including power generation, heating, ventilation, and air conditioning (HVAC), refrigeration, and aerospace engineering. Thermo-fluid-dynamics provides the tools and methods to model, analyze, and optimize these systems, allowing engineers and researchers to make informed decisions to improve their performance, energy efficiency, and environmental sustainability.
Power
It is well known that to heat water as well as the environment we need power. In small thermal systems, generally, power comes from a boiler that uses gas combustion in order to transfer energy to water. In the following, a condensing boiler is considered. It could heat both radiators or a puffer that works as technical water for primary exchangers.
The heat rate represents the power given by the source. The heat exchange in a hot water loop is related to the mass flow, specific heat capacity, and differences in temperature, as:
\dot Q=\dot mC_p\Delta T
Where m is the mass flow rate, cp is the heat capacity and ΔT the temperature difference.
\dot m\left[\text{kg}\over\text{s}\right]\\ \:\\c_p\left[\text{kJ}\over\text{kg K}\right]\\ \: \\\Delta T=T_{\operatorname{cold}}-T_{\operatorname{hot}}Then it is possible to compute the energy by the integration. In reality, due to the discrete temperature sampling, it is computed as the discrete left integration by summing all the i samples.
Q=\int_{t_{on}}^{t_{off}}\dot Q dt\approx \sum_{i=1}^{Q_i}\dot mC_p[\Delta T\Delta t]_iBoiler heat exchange
A condensing boiler is a type of heating system that is designed to maximize energy efficiency by recovering heat from the flue gases that would typically be lost in traditional boilers. It achieves this by utilizing the principles of condensation. When a fuel, such as natural gas or oil, is burned in a boiler, heat is produced as a result of the combustion process.
In a conventional boiler, the hot combustion gases are released through the flue and into the atmosphere, carrying away a significant amount of heat energy. In contrast, a condensing boiler is equipped with a secondary heat exchanger that is designed to extract as much heat as possible from the flue gases. These flue gases, which still contain a substantial amount of heat, are directed through the heat exchanger. The heat exchanger is usually made of a corrosion-resistant material, such as stainless steel or aluminum, that allows for efficient heat transfer.
Condensation
As the flue gases pass through the heat exchanger, they cool down significantly. The cooling causes the water vapor present in the flue gases to condense into liquid form. This process of condensation releases additional latent heat, which is captured and transferred to the system’s heating circuit. The condensed water, along with any impurities and acidic components from the flue gases, is then drained away through a condensate pipe. The condensate is typically acidic due to the presence of carbon dioxide, and measures are taken to neutralize its acidity to prevent damage to the boiler and the drainage system.
By recovering the latent heat from the flue gases through condensation, a condensing boiler can achieve higher energy efficiency compared to conventional boilers. This increased efficiency translates into lower fuel consumption, reduced greenhouse gas emissions, and cost savings for the user.
So, I can use thermal probes to read the temperature difference between supply and return pipes to compute the exchanged heat. I know that the heat transfer fluid is water with a heat capacity of 4.184 kJ kg−1K−1 . However, I need to discover the mass flow rate. Basically, I could use a flow sensor but it requires a sensor capable of supporting high temperature. Furthermore, it is a relatively constant quantity that doesn’t need to be continuously monitored as the temperature.
In this system under investigation, there are basically two states: pump ON with the corresponding mass flow rate and pump OFF with no flow, no heating, and no heat exchange. There is not a modulation control.
Pump and system curve
The heating loop is driven by the boiler pump through a 20mm copper pipe. By comparing the pump curve and the empirical system curve we can theoretically find the working point and so the current flow. This allows me to implement the computation of the energy obtained from the boiler in the smart home control system. The following images came from the Home Assistant dashboard but this is described in another article.
The pump needs to be able to satisfy the energy balance to pump the water from the desired height against friction losses. Following Bernoulli’s equation between the input (1) and the output (2) of the pipe:
\frac{p_1}{\rho g}+\frac{1}{2} \frac{v_1^2}{g}+h_1+H_{\text {Pump }}=\\\:\\\frac{p_2}{\rho g}+\frac{1}{2} \frac{v_2^2}{g}+h_2+H_{\text {friction }}In this closed loop, we can only consider the head losses due to friction. The head loss could be computed through empirical equations, related to fluid and flow properties, otherwise, we can use tabulated data. By considering the copper pipe pressure losses, an approximative loop length, and an appropriate factor to account for curve as well as restrictions we can find the working point of the pump, and so the current flow.
It requires some re-arrangements to fit the log-log curve and compare it with the pump one. Based on the specific curve of the products that constitute the system under investigation the fitting leads to:
Adaptive pump
While for the boiler I can rely on a constant flow rate pump and estimate the current flow from the curves for the solar panel control unit is a little different. It is designed to adjust the speed of the pump to meet the varying demands of a system in order to optimize energy efficiency and maintain desired system performance.
As described in a previous article, it relies on Triac dimming for the alternate current load as the 230V pump. I designed a little circuit to convert the 310 Vpk AC into a signal readable with Arduino. By measuring the signal amplitude in time I can compute the apparent and real power as well as the phase shift of the Triac dimming.
The analog-digital converter of Arudino allows sampling of the signal and computing the root mean square voltage with an appropriate frequency.
If you are interested in more details have a look at this article:
With an appropriate computation, it is possible to obtain the pump percentage of power as well as the corresponding pump curve.
Arduino NANO Code
float prop=900.000;
float power=0;
unsigned long previousMillis = 0;
const short sizeAV=10.00;
const short timeToSample=240; // it samples two times the pattern
float v=0;
unsigned long sum = 0;
#include<SoftwareSerial.h>
SoftwareSerial mySUART(4,5);
int readElSens() {
int n = 0;
unsigned long start_time = millis();
v=0;
sum=0;
while (millis() - start_time < 100) {
v = 310.000*analogRead(A1)/prop;
sum += v * v;
n++;
}
return ((double) (100*sum/(230.00*230.00*n)));
}
void setup() {
Serial.begin(9600);
mySUART.begin(9600);
}
void loop() {
power=0;
for (int i = 0; i < sizeAV; i++) {
power+=readElSens();
delay(20);
}
power=power/sizeAV;
if(power<10){
power=0;
}
if(power>10&&power<30){
power=30;
}
if(power>95){
power=100;
}
Serial.println(power);
mySUART.println(power);
delay(10000);
}
Thermal power estimation
The Arduino computation of the pump power percentage is sent through a serial port to a web server so that Home Assistant can read the signal thanks to the REST integration. However, the actual pump power is only the starting point. I have to convert it into the current flow by the intersection with the system curve. In addition, the pump curve varies depending on the current power.
I made a linear interpolation between the minimum and maximal pump curves corresponding to the minimal and maximal power, respectively 30% and 100% of the nominal power.
Fit pump data
Both curves are well represented through a parabolic interpolation of several points obtained from the manufacturer data.
parabola2 = Fit[pompa2, {1, x, x^2}, x];
parabola = Fit[pompa, {1, x, x^2}, x];
GraphicsRow[{Show[ListPlot[pompa2, PlotStyle -> Red],
Plot[parabola2, {x, 0, 5}]],
Show[ListPlot[pompa, PlotStyle -> Red], Plot[parabola, {x, 0, 5}]]}]
System estimation
The system could be estimated through the Darcy-Weisbach equation which describes the flow of a fluid through a porous conduct. The frictional loss is:
h=f{L\over D}{v^2\over 2g}Where v is the mean flow velocity, L is the conduct length, D is the diameter, and f is the Moody factor. The Moody friction factor is related to the pipe roughness through empirical relation. It is based on the Reynolds numbers:
\operatorname{Re}={\rho v D\over \mu}With the Reynold number the range of completely turbulent flow (Re>4000) reads:
f=\left[1.14+2\log_{10 }\left(D\over\epsilon\right)\right]^{-2}metri=40;
\[Epsilon] = 0.045;
f = (1.14 + 2 Log10[20/\[Epsilon]])^-2;
d = 0.02;
v [p_] := (p/1000)/(3600*Pi*d^2);
L = 40*metri;
h[p_] := f*L*(v)^2/(2*10*d)Here, no concentrated losses are taken into account and so the following pressure drop is mainly due to distributed friction.
Plot[h[p], {p, 0, 1000}]
Plot[{h[p], metri 10^-3*10^((1/c)*Log10[p] + N[Log10[0.00079]])}, {p,
0, 1000}]
System curve fit
Otherwise, we can estimate the system curve from the tabulated date of the pressure drop for unit length for the copper pipe of the corresponding diameter. They came with log-log curves.
The relation could be expressed as:
f(x)=10^{c\log_{10}(x)+\log_{10}(y_0)}.Where the y0 values represent the y values when log(x) is zero and the slope is:
c={\log_{10}(y_2)-\log_{10}(y_1)\over \log_{10}(x_2)-\log_{10}(x_1)}.Eventually, rescale unit of measure and length:
c = (Log[100] - Log[2])/(Log[4] - Log[0.550]);
d = N[shift /.
Solve[10^(c*Log10[1] + Log10[shift]) == 25, shift][[1, 1]]];
mbarToMeters = 0.010199773339984;
meters = 40;
mQtoL = 1/1000;
impianto[x_] := mbarToMeters*meters*10^(c*Log10[x*mQtoL] + Log10[d]);
GraphicsRow[{Plot[impianto[x], {x, 0, 10000}],
LogLogPlot[impianto[x], {x, 0, 10000}]}]
In the end, we obtain the system and pump curves as well as their intersection point.
Pump dynamic range
A linear relation is considered between the minimal and maximal power.
\Phi ={1000\over 3600}\left[0.6+(0.85-0.6){P-30\over100-30}\right]\quad \left[\text l/\text s\right]Φ represents the current flow and P is the actual pump power reader by the Arduino. Such information is then integrated into the Home Assistant platform to compute energy usage.
Flow sensor definition
A web sensor is used to read the pump power value from an HTTP GET call through the REST integration.
- platform: rest
resource: http://192.168.1.xxx/ratioAV
name: Percentuale potenza pompa solare
unique_id: potenza_percentuale_pompa_solare
force_update: trueThen, it is converted from power percentage to flow in liters per second. basically, it is the previously mentioned linear interpolation with a conversion factor written in the yaml languages of Home Assistant.
- sensor:
- name: "Portata pompa solare litriSec"
state: >
{% set relPower=states('sensor.percentuale_potenza_pompa_solare')|float(2) %}
{% set conv=1/3600 %}
{% set mqToLt=1000 %}
{% set minFlow=0.6 %}
{% set maxFlow=0.85 %}
{% set minPower=30 %}
{% set maxPower=100 %}
{% if relPower>5 %}
{{mqToLt*conv*(minFlow+(maxFlow-minFlow)*((relPower-minPower)/(maxPower-minPower))) }}
{% else %}
{{0}}
{% endif %}
unit_of_measurement: l/s
Low pass filtering
In addition, a low-pass filter is included to smooth the fast oscillations of the pump modulation. By observing a timespan of more than 3 hours, an"averaged" signal may be more useful. So, two virtual sensors have been added. The first one is a simple hourly mean while the second one is a low pass filter tuned for a less efficient smoothing.
sensor:
- platform: statistics
name: "Percentuale potenza pompa clean"
entity_id: sensor.percentuale_potenza_pompa_solare
state_characteristic: average_linear
max_age:
hours: 1
- platform: filter
name: "Filtered percentuale potenza pompa solare"
entity_id: sensor.percentuale_potenza_pompa_solare
filters:
- filter: lowpass
time_constant: 10
precision: 2Thermal energy
The last step is energy computation. Starting from the current power, it is computed including the heat capacity of 10% glycol and water solution, as mentioned before.
- sensor:
- name: "Potenza termica solare"
state: >
{% set flow=states('sensor.portata_pompa_solare_litrisec')|float(2) %}
{% set c=4000 %}
{{ (flow*c)*(states('sensor.tasmota_ds18b20_3_temperature') |float(2) - states('sensor.tasmota_sala_termica_ds18b20_5_temperature') |float(2)) }}
unit_of_measurement: WThen, a time integral is added with a left time rule to integrate the power over time. In addition, a utility meter type entity is added to split the time integration into a 1-hour time block.
utility_meter:
- hourly_energy_calore_pannelli_solare:
source: sensor.calore_pannelli_solare
name: Hourly Energy Calore Pannelli Solare
cycle: hourly
Finally, it is possible to compute the energy consumption or the be precise the amount of energy that both the gas boiler and the solar thermal system give to the technical water inside the puffer.
CSST pipe
Another interesting topic is the corrugated stainless steel tube used for the solar system. The corrugated profile allows bending the conduct and easily modeling it through the entire path. However, with these geometries, it may seem to present a very elevated pressure drop due to friction. Besides this consideration, we should account also for the low range of velocity of the thermal fluid.
The thermal fluid, typically water with the addition of glycol, circulates with a very low flow because it has to be slowly heated by the solar panel where a small mesh of capillaries is present. So, with a good calibration of the system and pump the friction due to pipe geometry shouldn't be a problem. However, a CFD analysis may be interesting.
Computational Fluid Dynamics (CFD) is a branch of fluid mechanics that involves the use of numerical methods and algorithms to solve and analyze problems related to fluid flow, heat transfer, and related phenomena. It combines principles from physics, mathematics, and computer science to simulate and understand the behavior of fluids in various applications.
CFD
In traditional fluid dynamics, the governing equations that describe the motion of fluids are often complex and difficult to solve analytically. CFD provides a computational approach to solving these equations by discretizing the domain into small elements (grid cells or mesh) and approximating the fluid properties and equations over these elements. The fluid flow equations, typically the Navier-Stokes equations, are then solved numerically on this discretized domain.
The previous figure shows a segment of the CSST pipe and the velocity field in the bi-dimensional half of the rotational-symmetric tube. In this particular problem, we can simplify the formulation accounting only for the 2D slice and its symmetry around the pipe axis. A velocity field is applied at the inlet with a parabolic profile obtained from the pump flow while a pressure BC is applied at the outlet. The segment is considered to be just before the pressure gauge which gives us the pressure condition.
A simple comparison between a corrugated tube and a smooth one is shown. They present the same internal diameter. The velocity field in the CSST show also a horizontal component due to the corrugation profile which accounts for friction losses. The third figure shows the differences between the two different fields in terms of velocity magnitude.
Thermal sensor
The aforementioned temperature probes are mainly DS18B20 digital temperature sensors based on One Wire connection. However, the following analysis is relatively valid also for PT100/1000 and other sensors with a stainless steel 6mm shell that appear really similar from a thermo-mechanical point of view.
There are some probe holders well used to put the temperature sensor nearest to the fluid of interest. Generally, they suggest using a thermal compound to increment the thermal conductivity and obtained more accurate measurements. Furthermore, discussion with the technical staff they suggest I use thermal paste also on sensors placed over the tube.
So, let's analyze what happens when we put a temperature probe with a stainless steel shell in contact with a copper pipe filled with hot water.
Thermal compound
The thermal compound is a chemical compound with high thermal conductivity used to enhance thermal dissipation generally in electronic equipment. It helps fill microscopic gaps to avoid trapping air in and so it allows us to reduce the response time of the temperature variation and let's see if it may be useful or not for this application.
As for a general thermal system, it satisfies the law of thermodynamics, and the heat is transferred from the hotter end to the colder one. The ability to transfer heat is linked to the thermal conductivity k and the spontaneous thermal flow is directed opposite to the temperature gradient, as expressed by Fuorier's law:
\mathbf q=-k\nabla T
From the conservation of energy, the rate of change of heat accumulation is equal to the spatial derivative of the heat flow:
{\partial Q\over \partial t}+{\nabla\cdot \mathbf q}=0So, by considering the time rate of change of heat stored as a function of the temperature and eventually a heat source qv, it reads:
\rho c_p {\partial T\over \partial t}-\nabla\cdot (k\nabla T)=\dot q_vIn the case of an isotropic and homogenous medium in space, without external or internal sources, it leads to the heat equation:
{\partial T \over \partial t}=-{\alpha}{\nabla^2 T}With the thermal diffusivity function of the thermal conductivity k, density ρ and specific heat capacity cp:
\alpha={k\over c_p\rho}Finite Element Method analysis
With the physical phenomenon well described we can try to solve the heat equation over the aforementioned domain: a probe in contact with a copper pipe filled with hot water. In addition, the probe is considered made of a stainless steel shell filled with air and all the elements are covered with a thermal insulator.
To solve the equation in these domains we need a numerical solver and I decided to use the Wolfram Mathematica PDE solver. The problem is analyzed in 2D representing a slice of the probes and the process of solution is made of different steps:
Variable definition
We start with the definition of the thermal properties of the different materials:
Subscript[\[Rho], air] =
QuantityMagnitude[ThermodynamicData["Air", "Density"]];
Subscript[Cp, air] =
QuantityMagnitude[
ThermodynamicData["Air", "IsobaricHeatCapacity"]];
Subscript[k, air] =
QuantityMagnitude[ThermodynamicData["Air", "ThermalConductivity"]];
Subscript[\[Rho], w] =
QuantityMagnitude[ThermodynamicData["Water", "Density"]];
Subscript[Cp, w] =
QuantityMagnitude[
ThermodynamicData["Water", "IsobaricHeatCapacity"]];
Subscript[k, w] =
QuantityMagnitude[ThermodynamicData["Water", "ThermalConductivity"]];
Subscript[\[Rho], Cu] = 8.96; (* [(g/(cm^3))] *)
Subscript[Cp, Cu] = 0.39*10^3; (* [(J/(kg K))] *)
Subscript[k, Cu] = 401 ; (*[(W/(m K))]*)
Subscript[\[Rho], steel] = 7.5; (* [(g/(cm^3))] *)
Subscript[Cp, steel] = 0.5*10^3; (* [(J/(kg K))] *)
Subscript[k, steel] = 16.2 ; (*[(W/(m K))]*)
Subscript[\[Rho], ins] = 24*10^-3; (* [(g/(cm^3))] *)
Subscript[Cp, ins] = 0; (* [(J/(kg K))] *)
Subscript[k, ins] = 0.036; (*[(W/(m K))]*)
Subscript[\[Rho], paste] = 1.7*10^-3; (* [(g/(cm^3))] *)
Subscript[Cp, paste] = 1.1*10^3; (* [(J/(kg K))] *)
Subscript[k, paste] = 3;(*[(W/(m K))]*) The geometrical properties of the domains representing the pipe, insulator and probe shell:
Subscript[d, pipe] = 0.015;
Subscript[r, pipe] = Subscript[d, pipe]/2;
Subscript[t, pipe] = 0.0008;
Subscript[d, probe] = 0.006;
Subscript[r, probe] = Subscript[d, probe]/2;
Subscript[t, probe] = 0.0004;
Subscript[t, insulation] = 0.02;
domain = Disk[{0, 0}, Subscript[r, pipe] + Subscript[t, insulation]];
pipe = Disk[{0, 0}, Subscript[r, pipe]];
water = Disk[{0, 0}, Subscript[r, pipe] - Subscript[t, pipe]];
probeShell =
Disk[{0, Subscript[r, pipe] + Subscript[r, probe]}, Subscript[r,
probe]];
probe = Disk[{0, Subscript[r, pipe] + Subscript[r, probe]},
Subscript[r, probe] - Subscript[t, probe]];
Domain definition
From the aforementioned geometries, a domain boundary is defined:
Subscript[b, insulation] = ToBoundaryMesh[domain];
Subscript[b, pipe] = ToBoundaryMesh[pipe];
Subscript[b, shell] = ToBoundaryMesh[probeShell];
Subscript[b, probe] = ToBoundaryMesh[probe];
Subscript[b, water] = ToBoundaryMesh[water];
bmesh = FEMUtils`BoundaryElementMeshJoin[Subscript[b, insulation],
Subscript[b, pipe], Subscript[b, shell], Subscript[b, water],
Subscript[b, probe]]As well as the marker highlighting different regions to explicate different thermal properties:
pipeCoordinate = {Subscript[r, pipe] - Subscript[t, pipe]/2, 0};
waterCoordinate = {0, 0};
shellCoordinate = {0, Subscript[r, pipe] + Subscript[t, probe]/2};
probeCoordinate = {0, Subscript[r, pipe] + Subscript[r, probe]};
insulationCoordinate = {(
Subscript[d, pipe] + Subscript[t, insulation])/2, 0};
markerCoordinates = {{waterCoordinate}, {pipeCoordinate}, \
{probeCoordinate}, {shellCoordinate}, {insulationCoordinate}};
Show[{bmesh["Wireframe"],
Graphics[
MapThread[{PointSize[0.02], #1, Point /@ #2} &, {markerColors,
markerCoordinates}]]}]
And so we can mesh all the regions:
regions = <|"water" -> 1, "pipe" -> 2, "shell" -> 3, "probe" -> 4,
"insulation" -> 5|>;
markersSpecification = {{waterCoordinate,
regions["water"]}, {pipeCoordinate,
regions["pipe"]}, {shellCoordinate,
regions["shell"]}, {probeCoordinate,
regions["probe"]}, {insulationCoordinate, regions["insulation"]}};
mesh = ToElementMesh[bmesh, "RegionMarker" -> markersSpecification,
"RegionHoles" -> {waterCoordinate}, "MaxCellMeasure" -> 0.0000001,
"MeshQualityGoal" -> 0.5];
markerColors = {Blue, Brown, LightGray, LightBlue, Gray};
GraphicsRow[{Show[
mesh["Wireframe"[
"MeshElementStyle" ->
Map[Directive[FaceForm[#], EdgeForm[]] &,
markerColors[[2 ;; 5]]]]],
Graphics[{markerColors[[1]], water}]],
mesh["Wireframe"[
Sequence[
PlotRange -> {{-1.2 Subscript[d, pipe]/2,
1.2 Subscript[d, pipe]/2}, {0, Subscript[d, pipe]}},
"MeshElementStyle" ->
Map[FaceForm[#] &, markerColors[[2 ;; 5]]]]]]}]Material parsing
To define the material properties I defined an auxiliary function:
parameterRegion[elementMarkers_, material_, parameter_] :=
With[{mrk1 = elementMarkers["water"],
mrk2 = elementMarkers["pipe"],
mrk3 = elementMarkers["probe"],
mrk4 = elementMarkers["insulation"],
mrk5 = elementMarkers["shell"],
var1 = material["water"][parameter],
var2 = material["pipe"][parameter],
var3 = material["probe"][parameter],
var4 = material["insulation"][parameter],
var5 = material["shell"][parameter]},
Which[ElementMarker == mrk1, var1, ElementMarker == mrk2, var2,
ElementMarker == mrk3, var3, ElementMarker == mrk4, var4,
ElementMarker == mrk5, var5]]It allows me to automatically associate the thermal properties with the different markers of different regions:
materials = <||>;
materials["water"] = <|"k" -> Subscript[k, w],
"\[Rho]" -> Subscript[\[Rho], w], "Cp" -> Subscript[Cp, w]|>;
materials["pipe"] = <|"k" -> Subscript[k, Cu],
"\[Rho]" -> Subscript[\[Rho], Cu], "Cp" -> Subscript[Cp, Cu]|>;
materials["shell"] = <|"k" -> Subscript[k, steel],
"\[Rho]" -> Subscript[\[Rho], steel],
"Cp" -> Subscript[Cp, steel]|>;
materials["probe"] = <|"k" -> Subscript[k, air],
"\[Rho]" -> Subscript[\[Rho], air], "Cp" -> Subscript[Cp, air]|>;
materials["insulation"] = <|"k" -> Subscript[k, ins],
"\[Rho]" -> Subscript[\[Rho], ins], "Cp" -> Subscript[Cp, ins]|>;
pars = <||>;
pars["MassDensity"] = parameterRegion[regions, materials, "\[Rho]"];
pars["SpecificHeatCapacity"] =
parameterRegion[regions, materials, "Cp"];
pars["ThermalConductivity"] =
parameterRegion[regions, materials, "k"];
W need to define the boundary and initial conditions before solving the problem.
Variables
Here, the temperature is a function of time and space (e.g. bi-dimensional problem)
vars = {T[t, x, y], t, {x, y}}Boundary conditions
Boundary conditions are required as well as the initial ones. The temperature is considered to start from a cold environment of 10°C and it represents also the boundary of the insulator. It is a typical scenario of winter time with the pipe in an external environment.
\begin{cases} T(0,x,y)=T_{\operatorname{cold}}\quad \in \Omega\\ T(t,x,y)=T_{\operatorname{cold}} \quad \in \partial \Omega_{o}\end{cases}op = HeatTransferPDEComponent[vars, pars];
Subscript[T, cold] = 10;
Subscript[T, hot] = 60;
ic = {T[0, x, y] == Subscript[T, cold]};
Subscript[\[CapitalGamma], water] =
HeatTemperatureCondition[
x^2 + y^2 == (Subscript[r, pipe] - Subscript[t, pipe])^2, vars,
pars, <|"SurfaceTemperature" -> Subscript[T, hot]|>];
Subscript[\[CapitalGamma], cold] =
HeatTemperatureCondition[
x^2 + y^2 == (Subscript[r, pipe] + Subscript[t, insulation])^2,
vars, pars, <|"SurfaceTemperature" -> Subscript[T, cold]|>];
pde = {op == 0, ic, Subscript[\[CapitalGamma], water],
Subscript[\[CapitalGamma], cold]};And we can simply solve the numerical problem:
tend = 50;
Tfun = NDSolveValue[pde, T, {t, 0, tend}, {x, y} \[Element] mesh];Do we really need thermal paste?
Besides the animation, we can use this numerical formulation to investigate the need for thermal paste. Let's add a domain representing the thermal paste and observe how much the temperature varies.
Paste geometry
Thermal paste is added as a disk around the pipe and probe with a dimension of 3/4 of the pipe one.
pasteSize = 3/4 Subscript[r, pipe];
pasteDomain =
RegionDifference[
RegionDifference[Disk[{0, 3/4 Subscript[r, pipe]}, pasteSize],
probeShell], pipe];
regions["paste"] = 6;
pasteCoordinateLeft = {-(pasteSize/2), Subscript[r, pipe]};
pasteCoordinateRight = {pasteSize/2, Subscript[r, pipe]};
markerCoordinates = {{waterCoordinate}, {pipeCoordinate}, \
{probeCoordinate}, {shellCoordinate}, {insulationCoordinate}, \
{pasteCoordinateLeft}, {pasteCoordinateRight}};
markerColors = {Blue, Brown, LightGray, Gray, LightBlue, LightRed,
Yellow};The thermal compound domain is added by subtracting the pipe and the probe from the aforementioned disk.
pasteBoundary =
ImplicitRegion[
x^2 + (y - 3/4 Subscript[r, pipe])^2 == pasteSize^2 &&
x^2 + y^2 >= (0.99 Subscript[r, pipe])^2 &&
x^2 + (y - Subscript[r, pipe] - Subscript[r,
probe])^2 >= (0.99 Subscript[r, probe])^2, {x, y}];
Subscript[b, paste] =
ToBoundaryMesh[
pasteBoundary, {{-Subscript[r, pipe] - Subscript[t, insulation],
Subscript[r, pipe] + Subscript[t,
insulation]}, {-Subscript[r, pipe] - Subscript[t, insulation],
Subscript[r, pipe] + Subscript[t, insulation]}},
"MaxBoundaryCellMeasure" -> 0.001];
Subscript[b, insulation] =
ToBoundaryMesh[domain, "MaxBoundaryCellMeasure" -> 0.01];
Subscript[b, pipe] =
ToBoundaryMesh[pipe, "MaxBoundaryCellMeasure" -> 0.01];
Subscript[b, shell] =
ToBoundaryMesh[probeShell, "MaxBoundaryCellMeasure" -> 0.01];
Subscript[b, probe] =
ToBoundaryMesh[probe, "MaxBoundaryCellMeasure" -> 0.01];
Subscript[b, water] =
ToBoundaryMesh[water,
"MaxBoundaryCellMeasure" -> 0.01]; bmeshWithPaste =
FEMUtils`BoundaryElementMeshJoin[Subscript[b, insulation], Subscript[
b, probe], Subscript[b, shell], Subscript[b, pipe], Subscript[b,
water], Subscript[b, paste] ]
GraphicsRow[{Show[Subscript[b, paste]["Wireframe"]],
Show[{bmeshWithPaste["Wireframe"],
Graphics[
MapThread[{PointSize[0.02], #1, Point /@ #2} &, {markerColors,
markerCoordinates}]]}]}]
New mesh
Also, new marker and material properties need to be specified.
markersSpecification = {{waterCoordinate,
regions["water"]}, {shellCoordinate,
regions["shell"]}, {pipeCoordinate,
regions["pipe"]}, {probeCoordinate,
regions["probe"]}, {pasteCoordinateLeft,
regions["paste"]}, {insulationCoordinate,
regions["insulation"]}};
meshWithPaste =
ToElementMesh[bmeshWithPaste, "RegionMarker" -> markersSpecification,
"RegionHoles" -> {waterCoordinate}, "MaxCellMeasure" -> 0.0000001,
"BoundaryMeshGenerator" -> RegionPlot, "MeshQualityGoal" -> 0.5]And let's plot the new mesh with the corresponding domains.
GraphicsRow[{Show[
meshWithPaste[
"Wireframe"[
"MeshElementStyle" ->
Map[Directive[FaceForm[#], EdgeForm[]] &,
markerColors[[2 ;; 6]]]]],
Graphics[{markerColors[[1]], water}]],
meshWithPaste[
"Wireframe"[
Sequence[
PlotRange -> {{-1.2 Subscript[d, pipe]/2,
1.2 Subscript[d, pipe]/2}, {0, Subscript[d, pipe]}},
"MeshElementStyle" ->
Map[FaceForm[#] &, markerColors[[2 ;; 6]]]]]]}]Add material properties
The same structure of the auxiliary function is used.
parameterRegion[elementMarkers_, material_, parameter_] :=
With[{mrk1 = elementMarkers["water"],
mrk2 = elementMarkers["pipe"],
mrk3 = elementMarkers["probe"],
mrk4 = elementMarkers["insulation"],
mrk5 = elementMarkers["shell"],
mrk6 = elementMarkers["paste"],
var1 = material["water"][parameter],
var2 = material["pipe"][parameter],
var3 = material["probe"][parameter],
var4 = material["insulation"][parameter],
var5 = material["shell"][parameter],
var6 = material["paste"][parameter]},
Which[ElementMarker == mrk1, var1, ElementMarker == mrk2, var2,
ElementMarker == mrk3, var3, ElementMarker == mrk4, var4,
ElementMarker == mrk5, var5, ElementMarker == mrk6, var6]]Furthermore, new material properties are added accounting for thermal compound conduvitity.
materials["paste"] = <|"k" -> Subscript[k, paste],
"\[Rho]" -> Subscript[\[Rho], paste],
"Cp" -> Subscript[Cp, paste]|>;
pars["MassDensity"] = parameterRegion[regions, materials, "\[Rho]"];
pars["SpecificHeatCapacity"] =
parameterRegion[regions, materials, "Cp"];
pars["ThermalConductivity"] =
parameterRegion[regions, materials, "k"];The new domain is pretty similar to the older one with the addition of the thermal compound around the junction between the probe and the pipe.
Boundary conditions
op = HeatTransferPDEComponent[vars, pars];
Subscript[T, cold] = 10;
Subscript[T, hot] = 60;
ic = {T[0, x, y] == Subscript[T, cold]};
Subscript[\[CapitalGamma], water] =
HeatTemperatureCondition[
x^2 + y^2 == (Subscript[r, pipe] - Subscript[t, pipe])^2, vars,
pars, <|"SurfaceTemperature" -> Subscript[T, hot]|>];
Subscript[\[CapitalGamma], cold] =
HeatTemperatureCondition[
x^2 + y^2 == (Subscript[r, pipe] + Subscript[t, insulation])^2,
vars, pars, <|"SurfaceTemperature" -> Subscript[T, cold]|>];
pde = {op == 0, ic, Subscript[\[CapitalGamma], water],
Subscript[\[CapitalGamma], cold]};Solve
tend = 50;
TfunPaste =
NDSolveValue[pde, T, {t, 0, tend}, {x, y} \[Element] meshWithPaste];The new solution shows a larger contour domain of higher temperature near the probe sheel highlighting how the thermal compound increases the conductivity, as expected.
However, a numerical investigation could be interesting to detect quantitative differences. Without thermal paste, the temperature reached in 6 seconds is 0.5°C lower in both the shell and the internal region of the probe.
Compute average
The average temperature in the probe may me more relevant:
\bar T_{\operatorname{shell}}={\iint_{\Omega_{\operatorname{shell}}}T(t,x,y) \:dxdy\over \iint_{\Omega_{\operatorname{shell}}}1\:dxdy}And the numerical computation:
shellDomain = RegionDifference[probeShell, probe];
Region[shellDomain];
npoints = 30;
lastIntegrationTime = 7;
probeArea = NIntegrate[1, {x, y} \[Element] shellDomain];
meanTprobe =
Table[{ t,
1/probeArea NIntegrate[
Tfun[t, x, y], {x, y} \[Element] shellDomain]}, {t, 0,
lastIntegrationTime, lastIntegrationTime/npoints}];
meanTprobeWithPaste =
Table[{ t,
1/probeArea NIntegrate[
TfunPaste[t, x, y], {x, y} \[Element] shellDomain]}, {t, 0,
lastIntegrationTime, lastIntegrationTime/npoints}];Plot data
ListLinePlot[{meanTprobe,
meanTprobeWithPaste, {{0,
0.99 Subscript[T, hot]}, {lastIntegrationTime,
0.99 Subscript[T, hot]}}, {{0,
0.95 Subscript[T, hot]}, {lastIntegrationTime,
0.95 Subscript[T, hot]}}, {{0,
Max[meanTprobe]}, {lastIntegrationTime, Max[meanTprobe]}}, {{3.32,
0}, {3.32, 57}}, {{2.94, 0}, {2.94, 57}}, {{4.1, 0}, {4.1,
58.8}}}, PlotRange -> {{2, lastIntegrationTime}, {50, 60}},
"PlotLabels" -> {"\!\(\*SubscriptBox[\(T\), \(av\)]\) no paste",
"\!\(\*SubscriptBox[\(T\), \(av\)]\) with paste", "99%", "95%"},
PlotStyle -> {Thick, Thick, Directive[Gray, Dashed],
Directive[Gray, Dashed], Directive[Gray, Dashed],
Directive[Gray, Dashed], Directive[Gray, Dashed],
Directive[Gray, Dashed]}]
However, in both cases, the temperature reaches 95% of the boundary temperature with a time difference of 0.7s over the 3.5s of the faster one.
Thermal compound application
But what could happen if we put too much thermal compound?
Test different geometry
Add a different geometry for the paste
probeShell =
Disk[{0, 1.1 Subscript[r, pipe] + Subscript[r, probe]}, Subscript[r,
probe]];
probe = Disk[{0, 1.1 Subscript[r, pipe] + Subscript[r, probe]},
Subscript[r, probe] - Subscript[t, probe]];
Subscript[b, shell] =
ToBoundaryMesh[probeShell, "MaxBoundaryCellMeasure" -> 0.01];
Subscript[b, probe] =
ToBoundaryMesh[probe, "MaxBoundaryCellMeasure" -> 0.01];
Subscript[b, water] =
ToBoundaryMesh[water,
"MaxBoundaryCellMeasure" -> 0.01]; bmeshWithPaste2 =
FEMUtils`BoundaryElementMeshJoin[Subscript[b, insulation], Subscript[
b, probe], Subscript[b, shell], Subscript[b, pipe], Subscript[b,
water], Subscript[b, paste] ]
markerCoordinates = {{waterCoordinate}, {pipeCoordinate}, \
{probeCoordinate}, {{0,
1.1 Subscript[r, pipe] +
Subscript[t, probe]/
2}}, {insulationCoordinate}, {pasteCoordinateLeft}, \
{pasteCoordinateRight}};
markerColors = {Blue, Brown, LightGray, Gray, LightBlue, LightRed,
Yellow};
GraphicsRow[{Show[Subscript[b, paste]["Wireframe"]],
Show[{bmeshWithPaste2["Wireframe"],
Graphics[
MapThread[{PointSize[0.02], #1, Point /@ #2} &, {markerColors,
markerCoordinates}]]}]}]
markersSpecification = {{waterCoordinate,
regions["water"]}, {{0,
1.1 Subscript[r, pipe] + Subscript[t, probe]/2},
regions["shell"]}, {pipeCoordinate,
regions["pipe"]}, {probeCoordinate,
regions["probe"]}, {pasteCoordinateLeft,
regions["paste"]}, {insulationCoordinate,
regions["insulation"]}};
meshWithPaste2 =
ToElementMesh[bmeshWithPaste2,
"RegionMarker" -> markersSpecification,
"RegionHoles" -> {waterCoordinate}, "MaxCellMeasure" -> 0.00000008,
"BoundaryMeshGenerator" -> RegionPlot, "MeshQualityGoal" -> 0.5]
markerColors = {Blue, Brown, LightGray, LightBlue, Gray, LightRed};
GraphicsRow[{Show[
meshWithPaste2[
"Wireframe"[
"MeshElementStyle" ->
Map[Directive[FaceForm[#], EdgeForm[]] &,
markerColors[[2 ;; 6]]]]],
Graphics[{markerColors[[1]], water}]],
meshWithPaste2[
"Wireframe"[
Sequence[
PlotRange -> {{-1.2 Subscript[d, pipe]/2,
1.2 Subscript[d, pipe]/2}, {0, Subscript[d, pipe]}},
"MeshElementStyle" ->
Map[FaceForm[#] &, markerColors[[2 ;; 6]]]]]]}]
Solve
tend = 50;
TfunPaste2 =
NDSolveValue[pde,
T, {t, 0, tend}, {x, y} \[Element] meshWithPaste2];
Comparison
Clearly, the little differences highlighted are not relevant in these applications where the temperature is read in time intervals of minutes over differences of a few seconds.
As we can see from the visual comparison of the temperature behavior in time:
Probe and shell temperature
It is also highlighted that in both cases with the thermal compound, the equilibrium temperature of the shell is 0.4 °C higher than the standard case. However, as mentioned before the time differences of a few seconds are not relevant for this particular application.
