Introduction: 2D Stormcloud, with lightning and lightning rod
In this notebook, we will model the path of a bolt of lightning when there is a lightning rod present. This simulation will create a grid of points and then solve Poisson's equation on that grid in order to find the electric potential throughout the 2D space. Based on that electric potential, it will predict where the lightning strike begins, then repeat the process to see how it advances.
After we've modelled the lightning strike we will see how a lightning rod setup will handle being struck by lightning, more specifically how it will handle the amount of energy in the strike and then evaluate whether the setups seem safe.
To calculate this we make simplified models of the rods and wires, calculate the resistances for different materials and gauges of wire, calculate for a range of currents in the strikes and finally compare the temperature rise in the system to the melting point of the PVC insulation around the wires.
First, we import the essential tools: numpy, matplotlib, pandas and numba to speed up calculations, and handle data.
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from numba import jit
import pandas as pd
Next, we define a function that will solve Poisson's equation on a rectangular grid, using a method called the "method of relaxation." This method relies on the fact that, as long as there are no charges in a region, the potential in that region will only change gradually. So, given a set of boundary conditions (places where the potential is specified or well defined) we can find the values of electric potential at certain points by averaging the potential of neighboring points. By doing this over and over, we eventually get a stable solution.
In the following function, the boundary conditions are given by the array b. In the points where b is defined (has a value), the boundary conditions are assumed to be this value. Where b is 'nan' (not defined—that is, we have yet to find it) we use the method of relaxation to find the value.
@jit(cache=True)
def solvepoisson(b,nrep):
# b = boundary conditions
# nrep = number of iterations
z = np.copy(b) # z = electric potential field
j = np.where(np.isnan(b)) #checks for where the points have no value, assigns them the value 0
z[j] = 0.0
znew = np.copy(z)
Lx = np.size(b,0) # determine the x range of the point grid
Ly = np.size(b,1) # determine the y range of the point grid
for n in range(nrep):
for ix in range(Lx):
for iy in range(Ly):
ncount = 0.0
pot = 0.0
if (np.isnan(b[ix,iy])): # check for cases in which the value is unspecified in the original grid
# Now, add up the potentials of all the the points around it
if (ix>0):
ncount = ncount + 1.0
pot = pot + z[ix-1,iy]
if (ix<Lx-1):
ncount = ncount + 1.0
pot = pot + z[ix+1,iy]
if (iy>0):
ncount = ncount + 1.0
pot = pot + z[ix,iy-1]
if (iy<Ly-1):
ncount = ncount + 1.0
pot = pot + z[ix,iy+1]
znew[ix,iy] = pot/ncount # Divide by the number of contributing
# surrounding points to find average potential
else:
znew[ix,iy]=z[ix,iy] # If the value is specified, keep it
tmp_z = znew # Swapping the field used for the calucaltions with the field from the previous iteration
znew = z # (to prepare for the next iteration)
z = tmp_z
return z
Now, we will use the poisson solver to simulate lightning. First, we will set up the boundary conditions, creating a grid of 50 by 50 points, setting all values to 'nan' (meaning they need to be solved for) except at the top (the cloud) where the potential is specified to be 2, and at the bottom (the ground) where it is specified to be zero.
We have to calculate the surface potential of the conductor first, this might seem daunting at first, we did a lot of back and forth to figure this out to be honest, untill we realised that we only need to deal with a constant potential. Since the lightning rod is connected to ground all charge within it and the conductor connecting it to ground will be evenly distributed across the surface of the rod, this means that the entire system will have a constant potential. We decided to set the potential to be $\frac{1}{10}$ of the top potential.
def set_conditions(Lx = 50, Ly = 50, L = 2, top_pot = 2.0, bot_pot = 0.0):
# First, we set up the boundary conditions
# and the necessary variables for our lightning rod
rod_pot = top_pot/10.0
rod_height = L/Ly
z = np.zeros((Lx,Ly),float)
b = np.copy(z)
c = np.copy(z)
b[:] = np.float('nan')
# Set the potential at the top of the grid to top_pot
b[:,0] = top_pot
# Set the potential at the bottom of the grid to bot_pot
b[:,Ly-1] = bot_pot
return z, b, c, Lx, Ly, rod_pot, rod_height
Now, to add in the lightning. To simulate the path of the lightning strike we will work backwards, starting with some charge on the ground and seeing where it moves to in order to get up to the cloud (this is simulating so-called "ground to cloud" lightning). This charge, because it is negative, prefers to stay at low potential, in this case V = 0. So, to find the path of the lightning, we will proceed as follows:
- Use the poisson solver to find the potential across the entire space, ground to cloud, starting with the ground
- Find where the charge is most likely to move to, based on the potential values multiplied by a random factor
- Set that location's potential equal to 0
- Update the neighboring positions to the lightning's path, making them possible locations for the lightning's next move
- We also set the potential where the lightning rod is placed, constant for each iteration of the loop.
In practice, this means we will be working with three arrays: The first holds the boundary values and the lightning's path (we call that one b; it was already defined, but will be updated based on the lightning's path). The second holds the possible places the lightning can move, stored as 'nan' values (we call that one zeroneighbor). The last holds the probabilistic values that lighting will move to particular positions (we call that one sprob).
# Create a function to calculate the bolt
def calc_bolt(z, b, c, Lx, Ly, rod_pot, rod_height):
# Create a copy of the boundary conditions matrix which will be used to check
#for possible locations for the lightning's path
zeroneighbor = np.copy(z)
zeroneighbor[:] = 0.0 #set all values in it equal to 0
#set the values next to the ground equal to 'nan'. This is where the lightning can start
zeroneighbor[:,Ly-2] = np.float('nan')
nrep = 3000 # Number of jacobi steps
eta = 1.0 #A factor that will be used in probability calculation
ymin = Ly-1 #The y value where we will stop (just above the ground)
ns = 0
i = 0
while (ymin>0):
# First find potential on the entire grid, based on the original boundary conditions
s = solvepoisson(b,nrep)
# Probability that lightning will move to a new position may depend on potential to power eta
sprob = s**eta
# We also multiply by a random number, uniform between 0 and 1, to introduce some randomness
# And we multiply with isnan(zeroneighbor) to ensure only 'nan' points can be chosen
sprob = sprob*np.random.uniform(0,1,(Lx,Ly))*np.isnan(zeroneighbor)
#now, find the point with max probability
[ix,iy] = np.unravel_index(np.argmax(sprob,axis=None),sprob.shape)
# Update the boundary condition array to set the potential where the lightning is to 0
b[ix,iy] = 0.0
# Update neighbor positions of the lightning path to 'nan'
#(making them possible choices for the next iteration)
if (ix>0):
zeroneighbor[ix-1,iy]=np.float('nan')
if (ix<Lx-1):
zeroneighbor[ix+1,iy]=np.float('nan')
if (iy>0):
zeroneighbor[ix,iy-1]=np.float('nan')
if (iy<Ly-1):
zeroneighbor[ix,iy+1]=np.float('nan')
x_pos = int(len(b[0])/2)
y_pos = len(b[1])
b[x_pos:x_pos+1,int(y_pos*(1-rod_height)):y_pos] = rod_pot
zeroneighbor[x_pos:x_pos+1,int(y_pos*(1-rod_height)):y_pos] = rod_pot
# sets the potential around the lightning rod in both b and zeroneighbour
# this is to make sure that the lightning goes past the rod.
ns = ns + 1
c[ix,iy] = ns #create an array of the lightning's path, scaled by the number of loops
if (iy<ymin): #iterate to the next set of y-values
ymin = iy
i+=1
if i > 10000:
break
return c, s, sprob
plt.figure(1,figsize=(20,15))
z, b, c, Lx, Ly, rod_pot, rod_height = set_conditions(Lx = 100, Ly = 100)
c1,s1,sprob1 = calc_bolt(z,b,c,Lx,Ly, rod_pot, rod_height)
z, b, c, Lx, Ly, rod_pot, rod_height = set_conditions(Lx = 100, Ly = 100)
c2,s2,sprob2 = calc_bolt(z,b,c,Lx,Ly, rod_pot, rod_height)
z, b, c, Lx, Ly, rod_pot, rod_height = set_conditions(Lx = 100, Ly = 100)
c3,s3,sprob3 = calc_bolt(z,b,c,Lx,Ly, rod_pot, rod_height)
z, b, c, Lx, Ly, rod_pot, rod_height = set_conditions(Lx = 100, Ly = 100)
c4,s4,sprob4 = calc_bolt(z,b,c,Lx,Ly, rod_pot, rod_height)
plt.subplot(4,2,1)
plt.imshow(c1.T, cmap="Purples_r") #create a plot of the lightning's path
plt.colorbar()
plt.subplot(4,2,2)
plt.imshow(s1.T, cmap="Purples_r") #create a plot of the final potential
plt.colorbar()
plt.subplot(4,2,3)
plt.imshow(c2.T, cmap="Purples_r") #create a plot of the lightning's path
plt.colorbar()
plt.subplot(4,2,4)
plt.imshow(s2.T, cmap="Purples_r") #create a plot of the final potential
plt.colorbar()
plt.subplot(4,2,5)
plt.imshow(c3.T, cmap="Purples_r") #create a plot of the lightning's path
plt.colorbar()
plt.subplot(4,2,6)
plt.imshow(s3.T, cmap="Purples_r") #create a plot of the final potential
plt.colorbar()
plt.subplot(4,2,7)
plt.imshow(c4.T, cmap="Purples_r") #create a plot of the lightning's path
plt.colorbar()
plt.subplot(4,2,8)
plt.imshow(s4.T, cmap="Purples_r") #create a plot of the final potential
plt.colorbar()
plt.tight_layout()
plt.show()
After running these simulations and plots dozens of times, we see that the majority of the lightning strikes will strike the lightning rod. Occasionally we get a lightning strike that does not strike the rod, and in these cases we noted that the strike struck quite a distance from the rod. This is due to the rod making a 45 degree "cone of protection" under it, as shown in the figure below.
This means that if we want to protect large structures, we have to set up several rods to protect it well enough. (see figure below with overhead ground wire)
We'll now look at these types of metals with conductivity $\sigma$, specific heat $H_{C}$ and density $\rho$.
Gold:
$\sigma = 4.1 \cdot 10^{7}\ S/m$
$H_{C} = 125.604\ J / (kg K)$
$\rho = 19300\ \text{kg/m}^{3}$
Copper:
$\sigma = 5.96 \cdot 10^{7}\ S/m$
$H_{C} = 376.812\ J / (kg K) $
$\rho = 8960\ \text{kg/m}^{3}$
Aliminium:
$\sigma = 3.77 \cdot 10^{7}\ S/m$
$H_{C} = 921.096\ J / (kg K)$
$\rho = 2600\ \text{kg/m}^{3}$
Silver:
$\sigma = 6.3 \cdot 10^{7}\ S/m$
$ H_{C} = 238.6476\ J / (kg K) $
$\rho = 10500\ \text{kg/m}^{3}$
Iron:
$\sigma = 1 \cdot 10^{7}\ S/m$
$ H_{C} = 460.548\ J / (kg K)$
$\rho = 7870\ \text{kg/m}^{3}$
We can calculate the resistance in a conductor with the formula: $R = \frac{L}{A\sigma}$ where $L$ is the length of the conductor, $A$ is the cross sectional area and $\sigma$ is the conductivity of the conductor.
We want to calculate the energy going through a conductor(the lightning rod and wire in this case), to do this we can use the formula $T = \frac{E}{H_{C}\ m}$, where E is the energy being passed through the system, $H_{C}$ is the specific heat of the conductor and $m$ is the mass of the conductor.
We know that $E = P t$ where $P = I^2R$, $t$ is the time the current is passed, inserting all this gives us:
We now have to find some baseline currents for our lightning strikes, through our research we found that lightning strikes can have very varying currents, from as low as 3kA to as high as 200kA. We decided to test the entire range of currents for a couple of different lengths of materials
Lightning rods are connected to ground by a lot of different means, some are connected to the plumbing of the building, some are connected to structural metal frames and some are simply connected to ground by a wire from the rod.
Estimating the amount of metal in a house's plumbing or structural frames is not that interesting in this case, since the amount of metal will easily be able to dissipate the energy from a lightning strike. We decided to look at some standard gauge wire, with PVC insulation, of different lengths.
We deem a setup safe if the temperature rise in the wire is below the melting point of the PVC insulation, which is 105 degrees celcius, if the conductor material heats up more than this we would deem it unsafe for use in this application.
# Here we define a dictionary with the values of the different properties
names = ["Gold", "Copper", "Aluminium", "Silver", "Iron"] # Conductor names
sigma = [(4.1* 10**7), (5.96* 10**7), (3.77* 10**7), (6.3* 10**7), (1* 10**7)] # S/m
spec_heat = [125.604, 376.812, 921.096, 238.6476, 460.548] # J / (kg K)
rho = [19300, 8960, 2600, 10500, 7870] # kg/m^3
materials = {} # setup an empty dictionary to populate
for i in range(len(names)): # populate the dictionary with the given values
values = {}
values["sigma"] = sigma[i]
values["spec_heat"] = spec_heat[i]
values["rho"] = rho[i]
materials[names[i]] = values
# Define a couple of functions for calculating Ohm's law and the temperature increase
# in a wire for a given wire
def resistance(L, A, sigma):
"""
Calculates the resistance of a given length of wire
"""
return (L / (A*sigma))
def temp(R, I, t, Hc, Vol, rho):
"""
Calculate the increase in temperature based
on the energy passed in the conductor
"""
return ((R * (I**2) * t) / (Hc * Vol * rho))
We found differing regulations about roof heights in Norway, from our research we decided that 4.5 meters and 8.0 meters were sensible values for single storey and two storey buildings respectively.
In our calculations we choose to assume that the cable connecting the rod to ground goes straight down to ground from the roof, to account for the length of wire that needs to be in the ground for a proper grounding to take place we assumed 0.5 metres would be sufficient.
# Setup the one storey and two storey house parameters.
one_storey_rod = 1.0 #m
two_storey_rod = 1.5 #m
rod_gauge = (0.01**2) * np.pi # m^2
one_storey_wire = 5.0 #m
two_storey_wire = 8.5 #m
wire_gauge = np.array([120, 95, 70, 55, 50, 35, 25, 16, 10, 6]) * (10**-6) #m^2
# These values are AWG 4/0, 3/0, 2/0, 1/0, 1, 2, 4, 6, 8 and 10 in m^2
for conductor in materials:
R_rod_os = resistance(one_storey_rod, rod_gauge, materials[conductor]["sigma"])# Calc resistance of one storey rod
R_rod_ts = resistance(two_storey_rod, rod_gauge, materials[conductor]["sigma"])# Calc resistance of two storey rod
R_os = resistance(one_storey_wire, wire_gauge, materials[conductor]["sigma"])# Calc resistance of one storey wire
R_ts = resistance(two_storey_wire, wire_gauge, materials[conductor]["sigma"])# Calc resistance of one storey wire
resistances = {}
resistances["os"] = (R_rod_os + R_os) # sums the rod and wire resistances
resistances["ts"] = (R_rod_ts + R_ts) # sums the rod and wire resistances
materials[conductor]["resistances"] = resistances # adds the resistances into our main dictionary
# Here we make a couple of nice looking tables using pandas
gauges = ["4/0","3/0","2/0","1/0","1","2","4","6","8","10"] # wire gauges
mat_list = [r"$Gold\ R(\Omega)$",r"$Copper\ R(\Omega)$",r"$Aluminium\ R(\Omega)$",
r"$Silver\ R(\Omega)$", r"$Iron\ R(\Omega)$"] # names with latex syntax for the plots and tables
resistances_os = []
resistances_ts = []
for conductor in materials:
resistances_os.append(materials[conductor]["resistances"]["os"])
resistances_ts.append(materials[conductor]["resistances"]["ts"])
res_table_os = pd.DataFrame(resistances_os,
columns=gauges,
index=mat_list)
res_table_ts = pd.DataFrame(resistances_ts,
columns=gauges,
index=mat_list)
def highlight_min(s):
'''
highlight the minimum in a Series.
'''
is_min = s == s.min()
return ['background-color: lightblue' if v else '' for v in is_min]
def highlight_max(s):
'''
highlight the max in a Series.
'''
is_max = s == s.max()
return ['background-color: red' if v else '' for v in is_max]
def color_val(val):
"""
Takes a scalar and returns a string with
the css property `'color: red'` for values over 1
strings, black otherwise.
"""
if val > 1:
color = 'pink'
elif 1 >= val > 0.75:
color = 'darkorange'
elif 0.75 >= val > 0.5:
color = 'orange'
elif 0.25 >= val > 0:
color = 'green'
else:
color = 'black'
return 'color: %s' % color
Here we can see a couple of tables which lists the resistances of our different situations, the first is the resistances of the different wire gauges, and materials, for one storey buildings, the second is for two storey buildings.
res_table_os.style.apply(highlight_min).apply(highlight_max)
res_table_ts.style.apply(highlight_min).apply(highlight_max)
# Here we make a quick graphical representation of the values shown above
plt.figure(figsize=(10,5))
plt.subplot(1,2,1)
plt.title("One Storey System")
for conductor in materials:
plt.plot(wire_gauge, materials[conductor]["resistances"]["os"], "-")
plt.plot(wire_gauge, materials[conductor]["resistances"]["os"], "o", label="%s"%(conductor))
plt.legend()
plt.xlabel(r"Cross Sectional Area of the wire $(m^{2})$")
plt.ylabel(r"$Resistance\ (\Omega)$")
plt.subplot(1,2,2)
plt.title("Two Storey System")
for conductor in materials:
plt.plot(wire_gauge, materials[conductor]["resistances"]["ts"], "-")
plt.plot(wire_gauge, materials[conductor]["resistances"]["ts"], "o", label="%s"%(conductor))
plt.legend()
plt.xlabel(r"Cross Sectional Area of the wire $(m^{2})$")
plt.ylabel(r"$Resistance\ (\Omega)$")
plt.tight_layout()
plt.show()
As we can see from these tables and plots Silver has the least resistance, followed by Copper, then Gold, then Aliminium and finally Iron. We can also take note that the resistance of the wires increases as the wire gauge becomes smaller. (4/0 is the largest wire, while 4 is the smallest)
# Calculates the volume of our conductors
for conductor in materials:
vol = {}
vol_os = (rod_gauge * one_storey_rod) + (wire_gauge * one_storey_wire)
vol_ts = (rod_gauge * two_storey_rod) + (wire_gauge * two_storey_wire)
vol["os"] = vol_os
vol["ts"] = vol_ts
materials[conductor]["vol"] = vol
volumes_os = [materials["Gold"]["vol"]["os"]]
volumes_ts = [materials["Gold"]["vol"]["ts"]]
vol_os = pd.DataFrame(volumes_os,
columns=gauges,
index=[r"Volume $m^{3}$"])
vol_ts = pd.DataFrame(volumes_ts,
columns=gauges,
index=[r"Volume $m^{3}$"])
vol_os
vol_ts
Here we can see the volumes of the different setups.
# Calculate the temperature rise
currents = np.linspace(3,200,101) * (10**3) # setting up the currents with a linspace from 3-200 kA
time_strike = 65 * (10**-6) # microseconds
for conductor in materials:
os = {} ; ts = {}
for i in range(len(materials[conductor]["resistances"]["os"])):
# Calculates the temperature rise in the one storey setup
temp_os = temp(R = materials[conductor]["resistances"]["os"][i],
I = currents,
t = time_strike,
Hc = materials[conductor]["spec_heat"],
Vol = materials[conductor]["vol"]["os"][i],
rho = materials[conductor]["rho"])
# Calculates the temperature rise in the Two storey setup
temp_ts = temp(R = materials[conductor]["resistances"]["ts"][i],
I = currents,
t = time_strike,
Hc = materials[conductor]["spec_heat"],
Vol = materials[conductor]["vol"]["ts"][i],
rho = materials[conductor]["rho"])
os["%s"%(i)] = temp_os
ts["%s"%(i)] = temp_ts
materials[conductor]["temps_os"] = os
materials[conductor]["temps_ts"] = ts
def plot(material,current):
plt.figure(figsize=(16,9))
print("%s: "%(material))
plt.subplot(1,2,1)
plt.title("One Storey System")
for i in range(len(materials[material]["temps_os"])):
plt.plot(current,materials[material]["temps_os"]["%s"%(i)], label="AWG: %s"%(gauges[i]))
plt.legend()
plt.xlabel(r"$Current\ (A)$")
plt.ylabel(r"$Temperature\ Rise\ (K)$")
plt.subplot(1,2,2)
plt.title("Two Storey System")
for i in range(len(materials[material]["temps_ts"])):
plt.plot(current,materials[material]["temps_ts"]["%s"%(i)], label="AWG: %s"%(gauges[i]))
plt.legend()
plt.xlabel(r"$Current\ (A)$")
plt.ylabel(r"$Temperature\ Rise\ (K)$")
plt.show()
plot("Gold", currents)
plot("Copper", currents)
plot("Aluminium", currents)
plot("Silver", currents)
plot("Iron", currents)
From these plots we can see that in all cases, the temperature goes up more the thinner the wire. This is quite intuitive, since there's less material to pass the energy through.
We can also note that copper seems to be the best conductor, it's temperatures are the lowest across all wire gauges, while iron is the worst due to it's temperatures being the highest across all wire gauges.
Now we will compare the maximum temperature increases to the melting point of the PVC insulation used on the wires.
# Gets the maximum temperature increase for each material and gauge
melt_point_pvc = 105 # degrees C
for conductor in materials:
max_temps_os = {}
max_temps_ts = {}
for i in range(len(materials[conductor]["temps_os"])):
# We add 20 degrees C to each max temp rise, since we assume we're working in an environment with
# a temperature of 20 degrees
max_temps_os["%s"%(i)] = np.amax(materials[conductor]["temps_os"]["%s"%(i)]) + 20 # add 20 degrees C
max_temps_ts["%s"%(i)] = np.amax(materials[conductor]["temps_ts"]["%s"%(i)]) + 20
materials[conductor]["max_temp_os"] = max_temps_os
materials[conductor]["max_temp_ts"] = max_temps_ts
melt_compare_os = []
melt_compare_ts = []
for conductor in materials:
temp_list_os = []
temp_list_ts = []
for i in range(len(materials[conductor]["max_temp_os"])):
temp_list_os.append((materials[conductor]["max_temp_os"]["%s"%(i)] / melt_point_pvc))
temp_list_ts.append((materials[conductor]["max_temp_ts"]["%s"%(i)] / melt_point_pvc))
melt_compare_os.append(temp_list_os)
melt_compare_ts.append(temp_list_ts)
melt_os = pd.DataFrame(melt_compare_os,
columns = gauges,
index = mat_list)
melt_ts = pd.DataFrame(melt_compare_ts,
columns = gauges,
index = mat_list)
melt_os.style.format("{:.2%}").apply(highlight_min).apply(highlight_max).applymap(color_val)
melt_ts.style.format("{:.2%}").apply(highlight_min).apply(highlight_max).applymap(color_val)
These tables show how close to the melting point of the PVC insulation that the cables got as a percentage of the melting point of the PVC insulation.
Conclusion
Lightning Rod Simulation
From our simulations we see that having a lightning rod present will in the majority of the cases cause the lightning to strike the rod. As we noted earlier the rod will make a 45 degree "cone of protection" under it, this means that one should take this into account when setting up lightning rods for practical applications.
Lightning Rod Materials
We chose to test a few metals for these calculations, Gold, Copper, Aluminium, Silver and Iron were our materials of choice.
From our calculations of resistance silver seemed to be the best candidate, due to it having the lowest resistance across all wire gauges. The difference in resistance between silver and the second best material(copper) was small, thus both are good conductors.
However when it came to the final result with the temperature increase, copper was clearly the best candidate, due to it having the lowest temperature increase across all wire gauges. The difference in temperature increase was not that large when it came to the bigger wire gauges, but as the wires became smaller the differences became significantly greater.
Most of the materials were well within safe temperature increase across the entire range of current, iron however was not for the thinnest gauges, aluminium and gold got dangerously close for the thinnest gauge.
For all applications copper seems to be the best alternative.