# Heat Transfer for engineers



## Stichedup (Aug 21, 2008)

Here is a Heat Transfer macro that I created in 2007 excel for an engineering class that I had at the University of Delaware.  It finds the 2D temperature distribution of a plate after 90 seconds, 360 seconds, and at steady state if the following conditions are applied to the plate:

1) two of the sides are at constant temperature
2) the other two sides have a constant heat flux (user input)
3) the fourier number is within limit (time step is not to large)
4) It is in celcius
5) the plate is 0.3 meters per side --> this can be changed in the macro

The user chosen inputs have to be put into these cells:
Cell A3 = Grid width, B3 = grid heigth --> This is the number of nodes across the plate.
cell B4 = heat transfer coeficient (237 for aluminum)
Cell B5 = specific heat (903 for aluminum)
Cell B6 = density (2702 for aluminum)
Cell B7 = time step *** you cannot make this to large or it will never converge.  0.5 [seconds] or lower
Cell B8 = heat generation
Cell B9 = Heat flex on the left side
Cell B10 = Heat flex on the bottom side
Cell B11 = Initial temperature of the plate
Cell B12 = steady temperature of the top side of the plate
Cell B13 = steady temperature of the right side of the plate
Cell B14 = Delta X =0.3/(A3-1)
Cell B15 = Delta Y =0.3/(B3-1)
Cell B16 = Alpha =B4/(B5*B6)


```
Sub HeatTransfer()

Application.ScreenUpdating = True               'This will increase the calculation time, but it will also let you know that the macro is working.
Application.Calculation = xlCalculationManual   'This is so excel does not do calculations in between every loop code --> Lowers calculation time
Application.EditDirectlyInCell = True           'Allows the worksheet's cell's format to be edited by the macro
ActiveSheet.Unprotect                           'Allow the macro to acces data from the worksheet

'PLEASE NOTE THAT "dx" must equal "dy"
'Defining constants
M = Cells(3, 1) 'Number of nodes in the x-direction
N = Cells(3, 2) 'Number of nodes in the y-direction
Const X = 0.3  'Length in the x-direction
Const Y = 0.3  'Length in the y-direction
Const ACC = 0.0001 'This is the converging accuracy without the cooling blocks
DeltaT = Cells(7, 2) 'This is the unit change in time
G = Cells(8, 2) 'Constant internal heat generation
Qleft = Cells(9, 2) 'Constant heat flux along the left side of the plate
Qbottom = Cells(10, 2) 'Constant heat flux along the bottom of the plate
Tinitial = Cells(11, 2) 'Initial temperature
Ttop = Cells(12, 2) 'The constant temperature along the top of the plate
Tright = Cells(13, 2) 'The constant temperature along the right side of the plate
'Aluminum Properties found for Pure Aluminum in appendix Table A-1 on page 840
Kth = Cells(4, 2) 'Thermal conductivity of aluminum in W/mºK
Cp = Cells(5, 2) 'Specific heat of aluminum in J/kgºK
Rho = Cells(6, 2) 'Density of Aluminum in kg/m^3
Alpha = Cells(16, 2) 'Is a property of aluminum equal to "Kth / (Rho * Cp)"
'Maxtrix
ReDim Temp(Cells(3, 1), Cells(3, 2)) As Single
ReDim T(Cells(3, 1), Cells(3, 2)) As Single
Dim isConverged As Boolean
Dim continue As Boolean
Dim dx As Single
Dim dy As Single
Dim dt As Single
Dim k As Single
Dim time As Single
'Variables
dx = X / (M - 1)
dy = Y / (N - 1)
dt = DeltaT
k = Kth
Fo = Alpha * dt / (dx * dy) 'For my own conveinience -> it is basically the Fourier number
time = 0
Dim Time90 As Single
Dim Time360 As Single
Dim Converge As Single
Dim dTime As Single
dTime = timer
Application.StatusBar = timer - dTime

'entering initial temperatures
For i = 1 To M
    For j = 1 To N
        T(i, j) = Tinitial
        Temp(i, j) = Tinitial
    Next j
Next i

'This is the first iterations used for the 90 second loop
MyCell = 1
Do
    Application.StatusBar = timer - dTime
    'Corner node at the bottom left
    i = 1
    j = N
    Temp(i, j) = 2 * Fo * (Qleft * dy / k + Qbottom * dx / k + T(i + 1, j) + T(i, j - 1) - 2 * T(i, j)) + 0.25 * G * Alpha * dt / k + T(i, j)
    'Top boundary condition -> constant temperature at Ttop
    For i = 1 To M
        j = 1
        Temp(i, j) = Ttop
    Next i
    'Right boundary condition -> constant temperature at Tright
    For j = 2 To N
        i = M
        Temp(i, j) = Tright
    Next j
    'Left boundary -> constant heat flux = 1000000 W/m^2
    For j = 2 To N - 1
        i = 1
        Temp(i, j) = Fo * (2 * Qleft * dy / k + T(i, j + 1) + T(i, j - 1) + 2 * T(i + 1, j) - 4 * T(i, j)) + 0.5 * G * Alpha * dt / k + T(i, j)
    Next j
    'Bottom boundary -> constant heat flux = 500000 W/m^2
    For i = 2 To M - 1
        j = N
        Temp(i, j) = Fo * (2 * Qbottom * dx / k + T(i + 1, j) + T(i - 1, j) + 2 * T(i, j - 1) - 4 * T(i, j)) + 0.5 * G * Alpha * dt / k + T(i, j)
    Next i
    'Center nodes -> constant thermal conductivity, density, and specific heat
    For j = 2 To N - 1
        For i = 2 To M - 1
            Temp(i, j) = Fo * (T(i + 1, j) + T(i - 1, j) + T(i, j + 1) + T(i, j - 1) - 4 * T(i, j)) + G * Alpha * dx * dy * dt / k + T(i, j)
        Next i
    Next j
    'Make the new temperatures the old temperatures
    For i = 1 To M
        For j = 1 To N
            T(i, j) = Temp(i, j)
        Next j
    Next i
    'displaying the time
    time = time + dt
    Cells(N + 2, 4) = time
Loop While (time < 90)
'display results in cells
For i = 1 To M
    For j = 1 To N
        Cells(j + 1, i + 3) = Temp(i, j)
    Next j
Next i
'timer
Time90 = timer - dTime
Cells(MyCell + N + 1, 18) = Time90
Time90 = timer

'Now the iterations will continue for the 360 sec loop
MyCell = N + 3
Do
    Application.StatusBar = timer - dTime
    'Corner node at the bottom left
    i = 1
    j = N
    Temp(i, j) = 2 * Fo * (Qleft * dy / k + Qbottom * dx / k + T(i + 1, j) + T(i, j - 1) - 2 * T(i, j)) + 0.25 * G * Alpha * dt / k + T(i, j)
    'Top boundary condition -> constant temperature at Ttop
    j = 1
    For i = 1 To M
        Temp(i, j) = Ttop
    Next i
    'Right boundary condition -> constant temperature at Tright
    i = M
    For j = 2 To N
        Temp(i, j) = Tright
    Next j
    'Left boundary -> constant heat flux = 100 W/m^2
    i = 1
    For j = 2 To N - 1
        Temp(i, j) = Fo * (2 * Qleft * dy / k + T(i, j + 1) + T(i, j - 1) + 2 * T(i + 1, j) - 4 * T(i, j)) + 0.5 * G * Alpha * dt / k + T(i, j)
    Next j
    'Bottom boundary -> constant heat flux = 50 W/m^2
    j = N
    For i = 2 To M - 1
        Temp(i, j) = Fo * (2 * Qbottom * dx / k + T(i + 1, j) + T(i - 1, j) + 2 * T(i, j - 1) - 4 * T(i, j)) + 0.5 * G * Alpha * dt / k + T(i, j)
    Next i
    'Center nodes -> constant thermal conductivity, density, and specific heat
    For j = 2 To N - 1
        For i = 2 To M - 1
            Temp(i, j) = Fo * (T(i + 1, j) + T(i - 1, j) + T(i, j + 1) + T(i, j - 1) - 4 * T(i, j)) + G * Alpha * dx * dy * dt / k + T(i, j)
        Next i
    Next j
    'Make the new temperatures the old temperatures
    For i = 1 To M
        For j = 1 To N
            T(i, j) = Temp(i, j)
        Next j
    Next i
    'displaying the time
    time = time + dt
    Cells(2 * N + 4, 4) = time
Loop While (time < 360)
'display results in cells
For i = 1 To M
    For j = 1 To N
        Cells(N + j + 3, i + 3) = Temp(i, j)
    Next j
Next i
'timer
Time360 = timer - Time90
Cells(MyCell + N + 1, 18) = Time360
Time360 = timer

'These iterations check for convergence (steady state answer)
MyCell = 2 * N + 5
Do
    Application.StatusBar = timer - dTime
    'Corner node at the bottom left
    i = 1
    j = N
    Temp(i, j) = 2 * Fo * (Qleft * dy / k + Qbottom * dx / k + T(i + 1, j) + T(i, j - 1) - 2 * T(i, j)) + 0.25 * G * Alpha * dt / k + T(i, j)
    'Top boundary condition -> constant temperature at Ttop
    j = 1
    For i = 1 To M
        Temp(i, j) = Ttop
    Next i
    'Right boundary condition -> constant temperature at Tright
    i = M
    For j = 2 To N
        Temp(i, j) = Tright
    Next j
    'Left boundary -> constant heat flux = 100 W/m^2
    i = 1
    For j = 2 To N - 1
        Temp(i, j) = Fo * (2 * Qleft * dy / k + T(i, j + 1) + T(i, j - 1) + 2 * T(i + 1, j) - 4 * T(i, j)) + 0.5 * G * Alpha * dt / k + T(i, j)
    Next j
    'Bottom boundary -> constant heat flux = 50 W/m^2
    j = N
    For i = 2 To M - 1
        Temp(i, j) = Fo * (2 * Qbottom * dx / k + T(i + 1, j) + T(i - 1, j) + 2 * T(i, j - 1) - 4 * T(i, j)) + 0.5 * G * Alpha * dt / k + T(i, j)
    Next i
    'Center nodes -> constant thermal conductivity, density, and specific heat
    For j = 2 To N - 1
        For i = 2 To M - 1
            Temp(i, j) = Fo * (T(i + 1, j) + T(i - 1, j) + T(i, j + 1) + T(i, j - 1) - 4 * T(i, j)) + G * Alpha * dx * dy * dt / k + T(i, j)
        Next i
    Next j
    'Check for convergence accuracy
    isConverged = True
    For i = 1 To M
        For j = 1 To N
            If Abs(T(i, j) - Temp(i, j)) > ACC Then isConverged = False
        Next j
    Next i
    'Make the new temperatures the old temperatures
    For i = 1 To M
        For j = 1 To N
            T(i, j) = Temp(i, j)
        Next j
    Next i
    'displaying the time
    time = time + dt
    Cells(3 * N + 6, 4) = time
Loop While (isConverged = False)
'display results in cells
For i = 1 To M
    For j = 1 To N
        Cells(2 * N + j + 5, i + 3) = Temp(i, j)
    Next j
Next i
'timer
Converge = timer - Time360
Cells(MyCell + N + 1, 18) = Converge
Totaltime = timer - dTime
Cells(MyCell + N + 2, 7) = Totaltime

MyCell = 2
For k = 1 To 3
    'Enter the energy balence formula into the top left cell
    Cells(MyCell, 3 + M + 2).Select
    ActiveCell.Value = "=$B$9*(0.5*$B$15)+$B$4*(0.5*$B$14)*(D3-D2)/$B$15+$B$4*(0.5*$B$15)*(E2-D2)/$B$14+$B$8*$B$14*$B$15/4"
    'Entering the energy balence formula into the bottom nodes (cells) at constant heat flux.  Because I cannot predict what grid size you are going to use, I am going to enter it in the third row and fill the equation down until it reaches the bottom.
    Cells(MyCell + 1, 3 + M + 3).Select
    ActiveCell.Value = "=$B$10*$B$14+$B$4*$B$14*(E2-E3)/$B$15+$B$4*(0.5*$B$15)*(D3+F3-2*E3)/$B$14+$B$8*$B$14*$B$15/2"
    Range(Cells(MyCell + 1, 3 + M + 3), Cells(MyCell + N - 1, 3 + M + 3)).FillDown
    Range(Cells(MyCell + N - 1, 3 + M + 3), Cells(MyCell + N - 1, 3 + 2 * M)).FillRight
    'Enter the energy balence formula into the bottom left corner.  Because I cannot predict what grid size you are going to use, I am going to enter it in the third row and fill the equation down until it reaches the bottom.
    Cells(MyCell + 1, 3 + M + 2).Select
    ActiveCell.Value = "=$B$9*(0.5*$B$15)+$B$10*(0.5*$B$14)+$B$4*(0.5*$B$14)*(D2-D3)/$B$15+$B$4*(0.5*$B$15)*(E3-D3)/$B$14+$B$8*$B$14*$B$15/4"
    Range(Cells(MyCell + 1, 3 + M + 2), Cells(MyCell + N - 1, 3 + M + 2)).FillDown
    'Enter the energy balence formula into the nodes (cells) along the left edge at constant heat flux
    Cells(MyCell + 1, 3 + M + 2).Select
    ActiveCell.Value = "=$B$9*$B$15+$B$4*(0.5*$B$14)*(D4+D2-2*D3)/$B$15+$B$4*$B$15*(E3-D3)/$B$14+$B$8*$B$14*$B$15/2"
    Range(Cells(MyCell + 1, 3 + M + 2), Cells(MyCell + N - 2, 3 + M + 2)).FillDown
    'Enter the energy balence formula into the top right corner.  Because I cannot predict the grid size that you are going to use, I am giong to enter it into a top node (cell) and fill the equation right until it reaches the right side.
    Cells(MyCell, 3 + M + 3).Select
    ActiveCell.Value = "=$B$4*(0.5*$B$14)*(E3-E2)/$B$15+$B$4*(0.5*$B$15)*(D2-E2)/$B$14+$B$8*$B$14*$B$15/4"
    Range(Cells(MyCell, 3 + M + 3), Cells(MyCell, 3 + 2 * M + 1)).FillRight
    'Enter the energy balence formula into the nodes (cells) along the top edge at a constant temperature
    Cells(MyCell, 3 + M + 3).Select
    ActiveCell.Value = "=$B$4*(0.5*$B$15)*(D2+F2-2*E2)/$B$14+$B$4*$B$14*(E3-E2)/$B$15+$B$8*$B$14*$B$15/2"
    Range(Cells(MyCell, 3 + M + 3), Cells(MyCell, 3 + 2 * M)).FillRight
    'Enter the energy balence formula into the bottom right corner.  Because I cannot predict the grid size that you are going to use, I am giong to enter it into a center node (cell) and fill the equation right until it reaches the right side.
    Cells(MyCell + 1, 3 + M + 3).Select
    ActiveCell.Value = "=$B$10*(0.5*$B$14)+$B$4*(0.5*$B$14)*(E2-E3)/$B$15+$B$4*(0.5*$B$15)*(D3-E3)/$B$14+$B$8*$B$14*$B$15/4"
    Range(Cells(MyCell + 1, 3 + M + 3), Cells(MyCell + 1, 3 + 2 * M + 1)).FillRight
    Range(Cells(MyCell + 1, 3 + 2 * M + 1), Cells(MyCell + N - 1, 3 + 2 * M + 1)).FillDown
    'Enter the energy balence formula into the nodes (cells) along the right edge at a constant temperature.  Because I cannot predict the grid size that you are going to use, I am giong to enter it into a center node (cell) and fill the equation right until it reaches the right side.
    Cells(3, 3 + M + 3).Select
    ActiveCell.Value = "=$B$4*(0.5*$B$14)*(E2+E4-2*E3)/$B$15+$B$4*$B$15*(D3-E3)/$B$14+$B$8*$B$14*$B$15/2"
    Range(Cells(MyCell + 1, 3 + M + 3), Cells(MyCell + 1, 3 + 2 * M + 1)).FillRight
    Range(Cells(MyCell + 1, 3 + 2 * M + 1), Cells(MyCell + N - 2, 3 + 2 * M + 1)).FillDown
    'Enter the energy balence formula into the center nodes (cells)
    Cells(3, 3 + M + 3).Select
    ActiveCell.Value = "=$B$4*$B$14*(E2+E4-2*E3)/$B$15+$B$4*$B$15*(D3+F3-2*E3)/$B$14+$B$8*$B$14*$B$15"
    Range(Cells(MyCell + 1, 3 + M + 3), Cells(MyCell + N - 2, 3 + 2 * M)).FillDown
    Range(Cells(MyCell + 1, 3 + M + 3), Cells(MyCell + N - 2, 3 + 2 * M)).FillRight
    MyCell = MyCell + N + 2
Next k

MyCell = 1
SpecialCell = 3
For p = 1 To 2
    MyCell = 1
    For k = 1 To 3
        'Formating the horizontal axis next to the data
        Range(Cells(MyCell, SpecialCell + 1), Cells(MyCell, M + SpecialCell)).Interior.Color = RGB(255, 255, 0)
        Range(Cells(MyCell, SpecialCell + 1), Cells(MyCell, M + SpecialCell)).Font.Color = Black
        Range(Cells(MyCell, SpecialCell + 1), Cells(MyCell, M + SpecialCell)).Font.Bold = True
        Range(Cells(MyCell, SpecialCell + 1), Cells(MyCell, M + SpecialCell)).BorderAround _
            Weight:=xlMedium
        Range(Cells(MyCell, SpecialCell + 1), Cells(MyCell, M + SpecialCell)).Borders(xlInsideVertical).Weight = xlMedium
        For i = 1 To M
            Cells(MyCell, i + SpecialCell) = 0 + 30 * (i - 1) / (M - 1)
        Next i
        'Formating the Vertical axis (row = C) next to the data
        Range(Cells(MyCell + 1, SpecialCell), Cells(N + MyCell, SpecialCell)).Interior.Color = RGB(255, 255, 0)
        Range(Cells(MyCell + 1, SpecialCell), Cells(N + MyCell, SpecialCell)).Font.Color = Black
        Range(Cells(MyCell + 1, SpecialCell), Cells(N + MyCell, SpecialCell)).Font.Bold = True
        Range(Cells(MyCell + 1, SpecialCell), Cells(N + MyCell, SpecialCell)).BorderAround _
            Weight:=xlMedium
        Range(Cells(MyCell + 1, SpecialCell), Cells(N + MyCell, SpecialCell)).Borders(xlInsideHorizontal).Weight = xlMedium
        For j = 1 To N
            Cells(j + MyCell, SpecialCell) = 30 - 30 * (j - 1) / (N - 1)
        Next j
        'Formating for the Cells under the data
        Range(Cells(MyCell + N + 1, SpecialCell), Cells(MyCell + N + 1, M + SpecialCell)).Interior.Color = Black
        Range(Cells(MyCell + N + 1, SpecialCell), Cells(MyCell + N + 1, M + SpecialCell)).Font.Size = 14
        Range(Cells(MyCell + N + 1, SpecialCell), Cells(MyCell + N + 1, M + SpecialCell)).Font.Color = RGB(255, 255, 255)
        Range(Cells(MyCell + N + 1, SpecialCell), Cells(MyCell + N + 1, M + SpecialCell)).Font.Bold = True
        Range(Cells(MyCell + N + 1, SpecialCell), Cells(MyCell + N + 1, M + SpecialCell)).Borders.Color = RGB(0, 255, 0)
        If SpecialCell > 4 Then GoTo line399
        'Formatting for the time cells
        Cells(MyCell + N + 1, SpecialCell) = "Time :"
        Cells(MyCell + N + 1, SpecialCell + 1).HorizontalAlignment = xlLeft
        'Formatting for the Grid Size cells
        Range(Cells(MyCell + N + 1, SpecialCell + 2), Cells(MyCell + N + 1, SpecialCell + 3)).MergeCells = True
        Cells(MyCell + N + 1, SpecialCell + 2) = "Grid Size :"
        Cells(MyCell + N + 1, SpecialCell + 2).HorizontalAlignment = xlRight
        Cells(MyCell + N + 1, SpecialCell + 4) = M & " X " & N
        Cells(MyCell + N + 1, SpecialCell + 4).HorizontalAlignment = xlLeft
        'Formatting for the Time Step cells
        Range(Cells(MyCell + N + 1, SpecialCell + 6), Cells(MyCell + N + 1, SpecialCell + 7)).MergeCells = True
        Cells(MyCell + N + 1, SpecialCell + 6) = "Time Step :"
        Cells(MyCell + N + 1, SpecialCell + 6).HorizontalAlignment = xlRight
        Cells(MyCell + N + 1, SpecialCell + 8) = dt
        Cells(MyCell + N + 1, SpecialCell + 8).HorizontalAlignment = xlLeft
        'Formatting for the Heat generation cells
        Range(Cells(MyCell + N + 1, SpecialCell + 9), Cells(MyCell + N + 1, SpecialCell + 10)).MergeCells = True
        Cells(MyCell + N + 1, SpecialCell + 9) = "Heat Gen:"
        Cells(MyCell + N + 1, SpecialCell + 9).HorizontalAlignment = xlRight
        Cells(MyCell + N + 1, SpecialCell + 11) = G
        Cells(MyCell + N + 1, SpecialCell + 11).HorizontalAlignment = xlLeft
        'Formatting for the Calculation Time cells
        Range(Cells(MyCell + N + 1, SpecialCell + 12), Cells(MyCell + N + 1, SpecialCell + 14)).MergeCells = True
        Cells(MyCell + N + 1, SpecialCell + 12) = "Calculation Time :"
line399:
        
        MyCell = MyCell + N + 2
    Next k
    SpecialCell = SpecialCell + M + 1
Next p
MyCell = 1
SpecialCell = 3

'Formatting for the total calculation timer
MyCell = 2 * N + 5
Range(Cells(MyCell + N + 2, 3), Cells(MyCell + N + 2, 6)).MergeCells = True
Range(Cells(MyCell + N + 2, 7), Cells(MyCell + N + 2, 8)).MergeCells = True
Range(Cells(MyCell + N + 2, 3), Cells(MyCell + N + 2, 8)).Interior.Color = Black
Range(Cells(MyCell + N + 2, 3), Cells(MyCell + N + 2, 8)).Font.Size = 14
Range(Cells(MyCell + N + 2, 3), Cells(MyCell + N + 2, 8)).Font.Color = RGB(255, 0, 0)
Range(Cells(MyCell + N + 2, 3), Cells(MyCell + N + 2, 8)).Font.Bold = True
Range(Cells(MyCell + N + 2, 3), Cells(MyCell + N + 2, 8)).Borders(xlInsideVertical).Color = RGB(255, 0, 0)
Range(Cells(MyCell + N + 2, 3), Cells(MyCell + N + 2, 8)).BorderAround _
    Weight:=xlMedium, Color:=RGB(255, 0, 0)
Cells(MyCell + N + 2, 3).HorizontalAlignment = xlRight
'Total calculation timer
Cells(MyCell + N + 2, 3) = "Total Calculation Time :"

'Conditional Formatting
Dim cfColorScale As ColorScale
MyCell = 1
For k = 1 To 3
    Range(Cells(MyCell + 1, 4), Cells(MyCell + N, 3 + M)).Select
    'Setting the selected cells to a conditional formatting of a three-color scale
    Set cfColorScale = Selection.FormatConditions.AddColorScale(ColorScaleType:=3)
    'Set the minimum threshold to red and maximum threshold to blue
    cfColorScale.ColorScaleCriteria(1).FormatColor.Color = RGB(190, 0, 0)
    cfColorScale.ColorScaleCriteria(2).FormatColor.Color = RGB(255, 255, 0)
    cfColorScale.ColorScaleCriteria(3).FormatColor.Color = RGB(0, 0, 255)
    MyCell = MyCell + N + 2
Next k
Cells(1, 3).Select

'resetting the original settings
Application.ScreenUpdating = True
Application.Calculation = xlCalculationAutomatic

End Sub
```


----------



## Stichedup (Aug 22, 2008)

Does anyone have any comments or concerns?  Has anyone else run the macro?  Is there anything that could/should be changed to improve calculation time?  Is the visual results ok?

Background:  The plate is a 30cm by 30cm aluminum plate that is 2.5 cm thick and is cooled on two sides with a 10ºC liquid water solution connected to a cooler.  The constant heat flux is applied by electric heating strips that are glued to the two other sides of the plate.  They have a digital power consumption meter.  The underside of the plate is resting on an insolating pad and the top is exposed to room temperature air.

Results:  It seams that the temperature convergence occurs to slowly in the macro. There is a 3-5% error in temperature if the time is scaled lower to appropriately match convergence time realized from the expiriment when a time step of 0.1 seconds is used on a 21x21 nodal grid.  I think this is pretty good, but requires to much math.

The other numbers to the right of the temperatures are an energy balence.


----------

