NPN voltage regulator

The voltage difference (9.1 V - 8.42 V) is right around 0.7 V, one silicon PN junction drop. This is because you're reading the zener voltage through Q5's base-emitter junction. With no load there's very little current so the voltage drop is low, with a load the PN junction starts to drop more voltage.

This is a common problem with pass transistors using a zener reference, you get one diode drop below the zener voltage or two diode drops if you use a Darlington configuration.

one way to avoid this is to use a lower voltage zener and a comparator circuit to compare the output taken through a voltage divider. The divider can be made variable to get different voltages.

how do I calculate the voltage at the base, emitter and collector of each transistor?

Look at the series potential divider formed by R4, R8, R9, R10 and, R11. Each one will have approximately 10 volts across it and, this plus the circa 10 volts on Q5's base adds up to 60 volts i.e. the input voltage. An emitter will be about 0.7 volts lower than the base on the same transistor.

LTSpice was showing me 49.4 V.

Yup, pretty much what I'd expect to see with a 60 volt supply.

Each transistor emitter lowers the collector voltage for the next lower transistor by 10 volts and, for a given output load current, all transistors will dissipate roughly the same power. And, because each NPN is rated at 40 volts, the peak voltage input could be safely up to 180 volts or a tad more.

It's a bit eccentric these days but workable.

I'm trying to understand this voltage regulation circuit.

Then redraw it:

^{simulate this circuit – Schematic created using CircuitLab}

Ideally, the base of each transistor doesn't require any recombination current. So the two green-colored currents shown above would be all there is to worry about.

Of course, each transistor does require some recombination current at their bases and the amount will vary with the load current. In addition, though this is only important for \$Q_5\$ in the above circuit, the \$V_{\small{BE}}\$ of each transistor will also vary with the load current.

In addition, the zener diode's voltage drop will vary with its own current. You can see this on its datasheet:

Assuming no base recombination currents for the moment, we'd expect to see about \$\frac{60\:\text{V}-9.1\:\text{V}}{5\,\cdot\,2.2\:\text{k}\Omega}\approx 4.63\:\text{mA}\$. This is very close to the test current shown and for which the nominal voltage is shown in the table. It appears that the designer targeted the right value in the schematic, so the first sniff test is passed.

The selected transistor, the 2N3904, has a worst-case \$\beta\approx 100\$ at a collector current of \$10\:\text{mA}\$, memory serving. (I'll leave it to you to verify.) It's typical is much higher. But that's all you can count on. And at different collector currents the worst-case \$\beta\$ is lower still.

The simulator, though, will use the nominal value of about \$\beta\approx 300\$ for the device and will use the exact same behavior for all 5 transistors.

Reality will be different. The transistors won't all behave the same and, more important to this circuit, their required base recombination currents won't be the same.

how do I calculate the voltage at the base, emitter and collector of each transistor?

This question begs another. Do you want to get values similar to what a simulator will give you? Or values that you'd get if you build one? You only talk about the simulator in your question, so that will guide my response.

You can just apply KCL to all the nodes. I'll label each node's index value as the same as the index value for the transistor that it connects to. For simplicity, I'll also just use \$R_1\$ for all the resistors, since they are the same value.

Technically, the load current as seen by the emitter of \$Q_5\$ will gradually split off into each of the base currents, such that the collector current of \$Q_1\$ will be slightly less. But let's make the assumption that all the base currents are the same and just call it \$I_{\small{B}}\$. For simulating reasonable loads in this circuit, that should be close enough. (But if you wanted to make this fit simulation for excessive loads, then the following development would need modification.)

So:

values = {vcc:60, r1:2200, zzt:15, vzt:9.1-15*5e-3, beta:300, vt:0.0259, isat:1e-14}
e1 = Eq( v1/r1 + v1/r1 + ib, vcc/r1 + v2/r1 )      # KCL V1
e2 = Eq( v2/r1 + v2/r1 + ib, v1/r1 + v3/r1 )       # KCL V2
e3 = Eq( v3/r1 + v3/r1 + ib, v2/r1 + v4/r1 )       # KCL V3
e4 = Eq( v4/r1 + v4/r1 + ib, v3/r1 + v5/r1 )       # KCL V4
e5 = Eq( v5/r1 + v5/zzt + ib, v4/r1 + vzt/zzt )    # KCL V5
eb = Eq( ib*(beta+1), iload )                      # Ie = Ib*(beta+1)
eshockley = Eq( vout, v5 - vt*ln(1+iload/isat) )   # Q5's emitter -- Shockley
ans=solve( [e1, e2, e3, e4, e5, eb, eshockley], [v1, v2, v3, v4, v5, ib, vout] )

There are a few simplifying assumptions that a simulator will not assume. But this should get close.

Let's assume what I see in your schematic, a \$20\:\text{mA}\$ load:

ans[vout].subs( parts | {iload:20e-3} ).n()
8.35783478415784

Now, this is in no small part due to the \$V_{\small{BE}}\$ of \$Q_5\$. There has to be a diode drop of sorts, following the Shockley equation. And there is. (Without any load, then the \$V_{\small{BE}}\$ of \$Q_5\$ will be very tiny and so the emitter voltage will be very close to the zener voltage at its base. So that explains why there may be so much difference when there's no load, at all.)

Let's look at the zener voltage itself with this same load:

ans[v5].subs( parts | {iload:20e-3} ).n()
9.09143074303690

And that's about right, keeping in mind that the current into the zener is a little lower than the datasheet specification. So far, nothing too surprising.

Another good question to wonder about is what the output impedance is. This is where we assume that there is a series resistance of some kind at the emitter of \$Q_5\$ so that if the load increases, then we'd expect to see a small voltage drop at the load itself.

Let's increase the load and see:

ans[vout].subs( parts | {iload:100e-3} ).n()
8.30420649639793

So it dropped from about \$8.36\:\text{V}\$ to about \$8.30\:\text{V}\$. And that is something you should expect.

There are two key parts that account for this. One part is something called \$r_e^{\:'}=\frac{V_T}{I_{\small{LOAD}}}\$. That says that the resistance varies with the load. And it does. The other factor is due to \$Z_{\small{ZT}}\$, as seen by \$Q_5\$, and reflected over to the emitter. It's value should be about 100 times smaller than \$Z_{\small{ZT}}=15\:\Omega\$, or about \$150\:\text{m}\Omega\$.

To test this idea, I'll first work out the load voltage, given a specified load current. Then I'll add a very tiny additional load current and recompute the new load voltage. Taking the difference between these and dividing by the tiny load current change will tell me the effective resistance at the emitter of \$Q_5\$ as seen by the load.

Don't forget! This will show both effects I mentioned. Not just one. Here's how I'll set this up:

dR = (ans[vout]-ans[vout].subs(iload,iload+di))/di      # sum of both effects
re = vt/iload                                           # little-re effect
remainder = dR - re                                     # should be 150 mOhms

Let's do some tests:

dR.subs( parts | {iload:20e-3,di:1e-6} )
1.44426569892511
dR.subs( parts | {iload:100e-3,di:1e-6} )
0.408296777870558

And it appears a heavier load does experience a lower apparent series resistance at \$Q_5\$'s emitter, as expected.

Let's see about the remainder, though. Will it be different for both loads? Or the same?

remainder.subs( parts | {iload:20e-3,di:1e-6} )
0.149265698925110
remainder.subs( parts | {iload:100e-3,di:1e-6} )
0.149296777870558

Not surprisingly, they are almost exactly the same. And that was predicted.

So there is an invariant piece of the emitter resistance seen by the load due to the zener diode and there is a load-dependent piece due to \$Q_5\$'s behavior.

I just want to remind you that I've neglected factors that the simulator doesn't neglect. So while the calculations here should roughly mimic simulator results, they will not match them exactly. There are, for example, bulk resistances at both the emitter and base of \$Q_5\$ that will also have an impact with the simulator results and that wasn't accounted for here.

But to answer your question, yes it is possible to compute something similar to what a simulator would report.

Just to confirm using LTspice:

And that appears to connect well with the above approach.

An alternative that "drops in" costs less and is probably more reliable.

Remove ALL resistors and transistors and replace with a single resistor.
Zener is rated at 100 mA max.
IF Vin is fixed at 60V Rmin = V/I = (60-5)/0.1 55 Ohms for 100 mA.
Rpower = 5.5W so a 10W aircooled resistor would work well.

This would be cheaper and almost certainly more reliable and the resistor would "solder in" to the existing circuit.

Rmin is controlled by Vin max.
I max at Vin_min = Imax x Vinmin/Vinmax.
For Vin variable R needs adjusting to achieve max zener dissipation at Vin max so R will be larger, so at Vinmin Iout will be lower.

Still probably worthwhile.
Doubling the zener power rating (if heatsinking adequate) doubles Ioutmax and needs a higher dissipation resistor.