Stray Voltage Calculator

Y.addline = function(Y, Zseries, from, to) {
    n = NROW(Zseries)
    fromNodes = from + 1:n - 1
    toNodes = to + 1:n - 1
    Yseries = solve(Zseries)
    Y[fromNodes,fromNodes] = Y[fromNodes,fromNodes] + Yseries # diagonal
    Y[toNodes,toNodes]     = Y[toNodes,  toNodes]   + Yseries
    Y[fromNodes,toNodes]   = Y[fromNodes,toNodes]   - Yseries # off diagonal
    Y[toNodes,fromNodes]   = Y[toNodes,  fromNodes] - Yseries
    Y
}

Y.addshunt = function(Y, Zshunt, busnum) {
    n = NROW(Zshunt)
    fromNodes = busnum:(busnum+n-1)
    Yshunt = solve(Zshunt)
    Y[fromNodes,fromNodes] = Y[fromNodes,fromNodes] + Yshunt # diagonal
    Y
}
load(RpadBaseFile("AA_conductors.RData"))
HTMLon()
cat("Neutral conductor: ")
HTMLselect("neutralConductor1", a$label, default=16)

numberOfSections = 50
subBus = 1
totalLength = totalLengthMi * 5.28 # ohms/kfeet
sectionLength = totalLength/numberOfSections # kfeet
Zgrnd = averageGround/(groundsPerKFeet * totalLength) *
        numberOfSections # ohms (average ground per section)
d.ab = 6 # feet (distance between the phase and neutral)

aa.n1 = a[a$label==neutralConductor1,]

r = 0.0001 # ohms per 1000 feet
gmr = 0.04 # feet
r.n = aa.n1$rac/5.28 # ohms per 1000 feet
gmr.n = aa.n1$gmr

Zcond = array(0i,c(2,2)) # ohms per 1000 feet
Zcond[1,1] = r + .01807*harmonic +
             harmonic*1i*.0529*log10(278.9*sqrt(rho/harmonic)/gmr)
Zcond[2,2] = r.n + .01807*harmonic +
             harmonic*1i*.0529*log10(278.9*sqrt(rho/harmonic)/gmr.n)
Zcond[1,2] = Zcond[2,1] = .01807*harmonic +
             harmonic*1i*.0529*log10(278.9*sqrt(rho/harmonic)/d.ab)

Zcond = Zcond * sectionLength # actual ohms
numberOfConductors = NROW(Zcond)

# preallocate the Ybus
n = (numberOfSections + 1)*numberOfConductors
Y = array(0i,c(n,n))

# make the Ybus from the impedances
for (i in seq(from=1,to=n-numberOfConductors,by=numberOfConductors))
    Y = Y.addline(Y,Zcond,from=i,to=i+numberOfConductors)

# add the shunt grounds
for (i in seq(from=2*numberOfConductors,to=n-1,by=numberOfConductors))
     Y = Y.addshunt(Y,Zgrnd,busnum=i) # the last one is the grounded neutral

# add the sub ground
Y = Y.addshunt(Y,Zsub,busnum=numberOfConductors*subBus) 

# ground the phase wire at the sub
Y = Y.addline(Y, .00001, from = 1+numberOfConductors*(subBus-1),
              to = numberOfConductors*subBus)

# make I
Isrc = array(0i,n)
# inject evenly distributed loads
if (loadType == "evenlyDistributed") {
    from = seq(from=1+numberOfConductors*subBus,to=n,by=numberOfConductors)
    # phase A
    to = seq(from=numberOfConductors*(1+subBus),to=n,by=numberOfConductors)
    # neutral
    Isrc[from] = UnbalanceI/(numberOfSections-subBus) # amperes
    Isrc[to] = -Isrc[from]
} else if (loadType == "lumpedAtEnd") {
    # inject isolated load at the end
    to = n - numberOfConductors + 1
    from = n
    Isrc[from] = Isrc[from] + UnbalanceI # amperes
    Isrc[to] = Isrc[to] - Isrc[from]
} else {
    Isrc[numberOfConductors*subBus] = Isrc[numberOfConductors*subBus] + 
                                      UnbalanceI # amperes
}
# Find the voltages:
V = solve(Y,Isrc)

Vn = V[seq(from=numberOfConductors,to=n,by=numberOfConductors)]

x = 0:numberOfSections * sectionLength / 5.28

idx.p = seq(from=1,to=n,by=numberOfConductors)
idx.n = seq(from=numberOfConductors,to=n,by=numberOfConductors)
V = rbind(V[idx.p] - V[idx.p + numberOfConductors],
          V[idx.n] - V[idx.n + numberOfConductors])
Ycond = solve(Zcond)
I = Ycond %*% V
Ip = I[1,]
In = I[numberOfConductors,]
Ig = Ip + In
cat("Substation Vn =",abs(Vn[subBus]),"V\n")

newgraph(width=6,height=3)
par(mfcol = c(1,2), mar = c(3, 3, 2, 1))

plot(x, abs(Vn),
     type = "l", ylim = c(0,2*abs(Vn[subBus])), 
     xlab = "", ylab = "")
mtext("Neutral-to-earth voltage, V",side=3,line=.7)
mtext("Distance from the sub, miles",side=1,line=-1,outer=TRUE)

plot(x, abs(Ig), 
     type = "l", col = "blue", ylim = c(0,UnbalanceI),
     xlab = "", ylab = "")
mtext("Currents, A",side=3,line=.7)
lines(x,abs(In),col="red")
lines(x,abs(Ip))
legend(.5, 5, c("phase","neutral","ground"), 
       col = c("black","red","blue"), 
       lty = 1, bty = "n", cex = .7)
HTMLon()
showgraph()