🔬
Solver Purpose & Physical Scope
Calculate linear and point thermal transmittances (ψ, χ), temperature factors (fRsi), and condensation risk. Features real-time 2D isothermal contours and FDM solver.
📂 Discipline: Conduction
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 383 times
📄 Source File:
thermal_bridge.f90
📁 calcul/Conduction /
thermal_bridge.f90
! ==============================================================================
! 2D Thermal Bridge Calculator (2D Pont Thermique)
! Steady-State Finite Difference Method (FDM) Solver
! References: ISO 10211, EN ISO 14683
! ==============================================================================
program thermal_bridge
implicit none
! Grid size
integer, parameter :: Nx = 50
integer, parameter :: Ny = 50
! Inputs
integer :: bridge_type ! 1=Balcony slab, 2=Wall corner, 3=Window frame, 4=Floor junction, 5=Column penetration
real(8) :: Ti, Te ! Inside and outside temperatures (C)
real(8) :: Rsi, Rse ! Surface resistances (m2-K/W)
real(8) :: bridge_material_k ! Conductivity of penetrating bridge material (W/m-K)
real(8) :: bridge_thick ! Bridge thickness (mm)
real(8) :: bridge_length ! Bridge length/extension (mm)
integer :: num_layers ! Number of wall layers
real(8) :: layer_thickness(10) ! Layer thickness (mm)
real(8) :: layer_k(10) ! Layer conductivity (W/m-K)
real(8) :: bridge_axis_length ! Length of bridge along wall axis (m)
! Grid and Solver parameters
real(8) :: H_domain ! Height of domain (m), fixed at 0.5 m
real(8) :: W_domain ! Width of domain (m), equal to total wall thickness
real(8) :: dx, dy
real(8) :: T(0:Nx+1, 0:Ny+1)
real(8) :: T_new(Nx, Ny)
real(8) :: k_grid(0:Nx+1, 0:Ny+1)
integer :: fluid_type(0:Nx+1, 0:Ny+1) ! 0=Solid, 1=Interior fluid, 2=Exterior fluid
! Derived parameters
real(8) :: layer_start(11)
real(8) :: wall_R_1D, wall_U_1D
real(8) :: q_2D, q_1D, psi, chi, Q_bridge, f_Rsi
real(8) :: T_si_min
! Loop & Solver variables
integer :: i, j, k, iter
real(8) :: x, y, error, max_diff
real(8) :: gW, gE, gN, gS, denom
H_domain = 0.5d0
! Read standard inputs
read(*, *) bridge_type
read(*, *) Ti
read(*, *) Te
read(*, *) Rsi
read(*, *) Rse
read(*, *) bridge_material_k
read(*, *) bridge_thick
read(*, *) bridge_length
read(*, *) num_layers
do i = 1, num_layers
read(*, *) layer_thickness(i)
read(*, *) layer_k(i)
enddo
read(*, *) bridge_axis_length
! Compute wall thickness and grid spacing
W_domain = 0.0d0
layer_start(1) = 0.0d0
do i = 1, num_layers
W_domain = W_domain + layer_thickness(i) / 1000.0d0 ! convert to meters
layer_start(i+1) = W_domain
enddo
dx = W_domain / dble(Nx)
dy = H_domain / dble(Ny)
! Setup fluid_type and default k_grid
fluid_type = 0
k_grid = 1.0d0
! Boundary default virtual neighbors
fluid_type(0, :) = 1 ! Left boundary virtual cells behave like Interior fluid (default)
fluid_type(Nx+1, :) = 2 ! Right boundary virtual cells behave like Exterior fluid (default)
! Map wall layers to k_grid
do i = 0, Nx+1
x = (dble(i) - 0.5d0) * dx
if (i == 0) x = 0.0d0
if (i == Nx+1) x = W_domain
! Find layer
k = 1
do j = 1, num_layers
if (x >= layer_start(j) .and. x <= layer_start(j+1)) then
k = j
exit
endif
enddo
k_grid(i, :) = layer_k(k)
enddo
! Apply Specific Thermal Bridge Geometries
if (bridge_type == 1) then
! 1. Balcony Slab: slab centered vertically, penetrates all layers
do j = 1, Ny
y = (dble(j) - 0.5d0) * dy
if (abs(y - H_domain/2.0d0) <= (bridge_thick/2000.0d0)) then
k_grid(1:Nx, j) = bridge_material_k
endif
enddo
else if (bridge_type == 2) then
! 2. Wall Corner: L-shaped wall
! Solid wall is for i <= 20 or j <= 20
! Interior room is for i > 20 and j > 20
do i = 1, Nx
do j = 1, Ny
if (i > 20 .and. j > 20) then
fluid_type(i, j) = 1 ! Interior fluid
endif
enddo
enddo
! For Corner, the outside boundaries at i=1 and j=1 are convective with Te
! The virtual cells at i=0 and j=0 behave like Exterior fluid
fluid_type(0, :) = 2
fluid_type(:, 0) = 2
! cutoff boundary at i=Nx (for j <= 20) and j=Ny (for i <= 20) are adiabatic
fluid_type(Nx+1, :) = 0 ! no fluid outside cutoff
fluid_type(:, Ny+1) = 0
else if (bridge_type == 3) then
! 3. Window Frame: wall meets window frame at j >= 25
! Frame is made of wood/PVC (low k) or metal (high k) with gas glazing
do j = 25, Ny
do i = 1, Nx
x = (dble(i) - 0.5d0) * dx
if (x < W_domain / 3.0d0) then
k_grid(i, j) = bridge_material_k ! frame outer/inner profiles
else if (x < 2.0d0 * W_domain / 3.0d0) then
k_grid(i, j) = 0.024d0 ! air/glazing cavity
else
k_grid(i, j) = bridge_material_k
endif
enddo
enddo
else if (bridge_type == 4) then
! 4. Floor Junction: concrete floor slab joins from interior (left)
! but does not fully penetrate the insulation layer.
do j = 1, Ny
y = (dble(j) - 0.5d0) * dy
if (abs(y - H_domain/2.0d0) <= (bridge_thick/2000.0d0)) then
do i = 1, Nx
x = (dble(i) - 0.5d0) * dx
if (num_layers >= 2) then
! Penetrates layer 1 (interior structural wall) but stops at layer 2 (insulation)
if (x < layer_start(2)) then
k_grid(i, j) = bridge_material_k
endif
else
! Penetrates halfway if single layer
if (x < W_domain / 2.0d0) then
k_grid(i, j) = bridge_material_k
endif
endif
enddo
endif
enddo
else if (bridge_type == 5) then
! 5. Column Penetration (Point thermal bridge)
! Represented in 2D as column cutting through center
do j = 1, Ny
y = (dble(j) - 0.5d0) * dy
if (abs(y - H_domain/2.0d0) <= (bridge_thick/2000.0d0)) then
k_grid(1:Nx, j) = bridge_material_k
endif
enddo
endif
! Initialize temperatures based on fluid_type
T = (Ti + Te) / 2.0d0
T_new = (Ti + Te) / 2.0d0
do i = 0, Nx+1
do j = 0, Ny+1
if (fluid_type(i, j) == 1) then
T(i, j) = Ti
else if (fluid_type(i, j) == 2) then
T(i, j) = Te
endif
enddo
enddo
! ==========================================================================
! Gauss-Seidel Solver Loop
! ==========================================================================
do iter = 1, 15000
max_diff = 0.0d0
do i = 1, Nx
do j = 1, Ny
! Skip fluid cells (their temperatures are fixed)
if (fluid_type(i, j) /= 0) cycle
! East neighbor link (i+1, j)
if (i == Nx .and. bridge_type == 2) then
! corner cutoff boundary is adiabatic
gE = 0.0d0
else
if (fluid_type(i+1, j) == 0) then
gE = (2.0d0 * k_grid(i, j) * k_grid(i+1, j)) / (k_grid(i, j) + k_grid(i+1, j)) / (dx**2)
else if (fluid_type(i+1, j) == 1) then
gE = 1.0d0 / (Rsi + dx / (2.0d0 * k_grid(i, j))) / dx
else
gE = 1.0d0 / (Rse + dx / (2.0d0 * k_grid(i, j))) / dx
endif
endif
! West neighbor link (i-1, j)
if (fluid_type(i-1, j) == 0) then
gW = (2.0d0 * k_grid(i, j) * k_grid(i-1, j)) / (k_grid(i, j) + k_grid(i-1, j)) / (dx**2)
else if (fluid_type(i-1, j) == 1) then
gW = 1.0d0 / (Rsi + dx / (2.0d0 * k_grid(i, j))) / dx
else
gW = 1.0d0 / (Rse + dx / (2.0d0 * k_grid(i, j))) / dx
endif
! North neighbor link (i, j+1)
if (j == Ny) then
! top boundary is adiabatic
gN = 0.0d0
else
if (fluid_type(i, j+1) == 0) then
gN = (2.0d0 * k_grid(i, j) * k_grid(i, j+1)) / (k_grid(i, j) + k_grid(i, j+1)) / (dy**2)
else if (fluid_type(i, j+1) == 1) then
gN = 1.0d0 / (Rsi + dy / (2.0d0 * k_grid(i, j))) / dy
else
gN = 1.0d0 / (Rse + dy / (2.0d0 * k_grid(i, j))) / dy
endif
endif
! South neighbor link (i, j-1)
if (j == 1) then
gS = 0.0d0
else
if (fluid_type(i, j-1) == 0) then
gS = (2.0d0 * k_grid(i, j) * k_grid(i, j-1)) / (k_grid(i, j) + k_grid(i, j-1)) / (dy**2)
else if (fluid_type(i, j-1) == 1) then
gS = 1.0d0 / (Rsi + dy / (2.0d0 * k_grid(i, j))) / dy
else
gS = 1.0d0 / (Rse + dy / (2.0d0 * k_grid(i, j))) / dy
endif
endif
denom = gE + gW + gN + gS
if (denom > 0.0d0) then
T_new(i, j) = (gE * T(i+1, j) + gW * T(i-1, j) + gN * T(i, j+1) + gS * T(i, j-1)) / denom
error = abs(T_new(i, j) - T(i, j))
if (error > max_diff) max_diff = error
endif
enddo
enddo
! Update T array
do i = 1, Nx
do j = 1, Ny
if (fluid_type(i, j) == 0) then
T(i, j) = T_new(i, j)
endif
enddo
enddo
! Check convergence
if (max_diff < 1.0d-7) exit
enddo
! ==========================================================================
! Heat flow integrations & 1D reference calculation
! ==========================================================================
! 1D Wall R-value
wall_R_1D = Rsi + Rse
do i = 1, num_layers
wall_R_1D = wall_R_1D + (layer_thickness(i) / 1000.0d0) / layer_k(i)
enddo
wall_U_1D = 1.0d0 / wall_R_1D
! Integrate interior heat flow q_2D (W/m)
q_2D = 0.0d0
T_si_min = Ti
if (bridge_type == 2) then
! L-shaped wall corner heat flow integration on the solid boundary cells
! adjacent to the room (i=20 for j > 20 and j=20 for i > 20)
do j = 21, Ny
q_2D = q_2D + (Ti - T(20, j)) / (Rsi + dx / (2.0d0 * k_grid(20, j))) * dy
if (T(20, j) < T_si_min) T_si_min = T(20, j)
enddo
do i = 21, Nx
q_2D = q_2D + (Ti - T(i, 20)) / (Rsi + dy / (2.0d0 * k_grid(i, 20))) * dx
if (T(i, 20) < T_si_min) T_si_min = T(i, 20)
enddo
! corner reference 1D flow based on internal lengths
q_1D = wall_U_1D * (dble(Nx - 20) * dx + dble(Ny - 20) * dy) * (Ti - Te)
else
! Rectangular domain heat flow integration at i=1
do j = 1, Ny
q_2D = q_2D + (Ti - T(1, j)) / (Rsi + dx / (2.0d0 * k_grid(1, j))) * dy
if (T(1, j) < T_si_min) T_si_min = T(1, j)
enddo
q_1D = wall_U_1D * H_domain * (Ti - Te)
endif
! Transmittance calculations
psi = 0.0d0
chi = 0.0d0
if (bridge_type == 5) then
! Point thermal bridge (Column penetration)
! Scale to W/K using the 2D plane area slice
chi = (q_2D * H_domain - q_1D) / (Ti - Te)
Q_bridge = chi * (Ti - Te)
else
! Linear thermal bridges (Balcony slab, Corner, Window joint, Floor junction)
psi = (q_2D - q_1D) / (Ti - Te)
Q_bridge = psi * bridge_axis_length * (Ti - Te)
endif
! Condensation factor f_Rsi
f_Rsi = (T_si_min - Te) / (Ti - Te)
! ==========================================================================
! Output report
! ==========================================================================
print *, "=========================================================================="
print *, " THERMAL BRIDGE 2D FDM ANALYSIS REPORT (ISO 10211) "
print *, "=========================================================================="
select case (bridge_type)
case (1)
print *, " Bridge Type: Balcony slab penetration (Linear)"
case (2)
print *, " Bridge Type: Wall Corner junction (Linear)"
case (3)
print *, " Bridge Type: Window frame joint (Linear)"
case (4)
print *, " Bridge Type: Floor slab junction (Linear)"
case (5)
print *, " Bridge Type: Column penetration (Point)"
end select
print "(A, F8.2, A)", " Inside Temperature (Ti): ", Ti, " C"
print "(A, F8.2, A)", " Outside Temperature (Te):", Te, " C"
print "(A, F8.4)", " Inside Resistance (Rsi): ", Rsi
print "(A, F8.4)", " Outside Resistance (Rse):", Rse
print *
print *, "--------------------------------------------------------------------------"
print *, " UNDISTURBED 1D WALL PERFORMANCE"
print *, "--------------------------------------------------------------------------"
print "(A, F10.3, A)", " Total Wall Thickness: ", W_domain * 1000.0d0, " mm"
print "(A, F10.4, A)", " 1D Thermal Resistance: ", wall_R_1D, " m2-K/W"
print "(A, F10.4, A)", " 1D Heat Transmittance U: ", wall_U_1D, " W/m2-K"
print *
print *, "--------------------------------------------------------------------------"
print *, " 2D NUMERICAL THERMAL BRIDGE RESULTS"
print *, "--------------------------------------------------------------------------"
print "(A, F10.3, A)", " Total 2D Heat Flow: ", q_2D, " W/m"
print "(A, F10.3, A)", " Undisturbed 1D Heat Flow:", q_1D, " W/m"
if (bridge_type == 5) then
print "(A, F10.5, A)", " Point Transmittance (chi):", chi, " W/K"
else
print "(A, F10.5, A)", " Linear Transmittance (psi):", psi, " W/(m-K)"
print "(A, F10.2, A)", " Thermal Bridge Length: ", bridge_axis_length, " m"
endif
print "(A, F10.2, A)", " Additional Heat Loss: ", Q_bridge, " W"
print "(A, F10.2, A)", " Min Inside Surface Temp: ", T_si_min, " C"
print "(A, F10.4)", " Temp Factor (f_Rsi): ", f_Rsi
if (f_Rsi < 0.75d0) then
print *, " WARNING: Condensation & mold growth risk detected (f_Rsi < 0.75)!"
else
print *, " Surface temperature factor is safe (f_Rsi >= 0.75)."
endif
print *, "=========================================================================="
end program thermal_bridge
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 thermal_bridge.f90 -o thermal_bridge_calc
2. Execution with input.txt redirection:
thermal_bridge_calc < input.txt
📄 Sample input.txt File Structure
Sample Data:
1 20.0 -10.0 0.13 0.04 1.7 200.0 300.0 2 150.0 1.5 50.0 0.03 1.0
Parameter Description:
Bridge Type (1=Concrete column bridge, 2=Linear slab bridge) Inner temp [°C] Outer temp [°C] Inner film resistance Rsi [m2-K/W] Outer film resistance Rse [m2-K/W] Bridge thermal conductivity k_bridge [W/m-K] Bridge thickness t_bridge [mm] Bridge width/length w_bridge [mm] Number of layers Layer 1 thickness [mm] Layer 1 thermal conductivity [W/m-K] Layer 2 thickness [mm] Layer 2 thermal conductivity [W/m-K] Length along bridge axis L [m]