! ============================================================================== ! 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