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)
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)
Code:
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
