How to derive the torque expression for a BLDC?
Brushless AC

Lets start with a brushless AC machine first &amp; Field Orientated Control.

So BLAC machine's have their stator windings sinusoidally distributed (higher concentration close to the tooth, lower for the outer turns)

Now to control such machine is relatively complex and you can use a field orientated control ( F.O.C. ).Using Park's transform (Clark + rotation) the 3phase sinus stator current's can be transformed 1st into a rotating 2phase representation 

Park Transforms - General

$$
I_{\alpha \beta 0 } = \frac{2}{3}\begin{bmatrix}
1 &amp; \frac{-1}{2} &amp; \frac{-1}{2} \\ 
0 &amp; \frac{\sqrt{3}}{2} &amp;\frac{-\sqrt{3}}{2} \\ 
\frac{1}{2} &amp; \frac{1}{2} &amp; \frac{1}{2}
\end{bmatrix}
$$

or simply 

\$I_\alpha = I_a \$

\$I_\beta = \frac{2I_b + Ia}{\sqrt{3}} \$

And these phasors can further be reduced to two DC quantities via a rotating frame of reference transform

$$
\begin{bmatrix}
I_d\\ 
I_q
\end{bmatrix} = \begin{bmatrix}
Cos(\Theta ) &amp; Sin(\Theta )\\ 
-Sin(\Theta ) &amp; Cos(\Theta )
\end{bmatrix} \cdot \begin{bmatrix}
I_\alpha\\ 
I_\beta
\end{bmatrix}
$$

These two DC terms (with Id usually controlled to 0, unless field weakening is desired) are simple inputs to two classic PI loops to control Id and Iq (the output of which is Vd &amp; Vq).
This Vd &amp; Vq, via inverse Clark&amp;PArk produce the 3phase voltage that will be applied to the stator (insert SVM or SPWM)

Great, it works 

Brushless DC &amp; general Permanent Magnet machine equations

Thing is BLDC machines produce higher torque than a sinusoidally wound stator (downto the higher concentration of winding around the teeth to produce the flattened line-line profile). The downside is the torque ripple. This is mostly due to simplifying the control downto which equally limits the effective bandwidth of the controller.

A park-like transform can be applied to overcome some of these shortcomings as the aim is to produce a stator stimulus closer to the airgap profile - tpyically a Quazi squarewave (30-120-30) is superimposed onto a trapezoidal backEMF (60-60-60). 

First of all you need the machine equation for a BLDC machine (which is the same for a BLAC machine).
$$
\begin{bmatrix}
V_a\\ 
V_b\\ 
V_c
\end{bmatrix} = R_s \begin{bmatrix}
i_a\\ 
i_b\\ 
i_c
\end{bmatrix} + L\frac{\mathrm{d} }{\mathrm{d} t}\begin{bmatrix}
i_a\\ 
i_b\\ 
i_c
\end{bmatrix} + \begin{bmatrix}
e_a\\ 
e_b\\ 
e_c
\end{bmatrix}$$

\$e_a, e_b, e_c\$ are the 3 backEMF's that are generated for a rotating rotor. Via Faraday's law
$$
\varepsilon = - \frac{\mathrm{d} \Phi _B}{\mathrm{d} t}
$$

BackEMF is the rate of change of flux

$$
\frac{\mathrm{d} \Phi }{\mathrm{d} t} = \frac{\mathrm{d} \Phi }{\mathrm{d} \Theta}\frac{\mathrm{d} \Theta }{\mathrm{d} t} = \frac{\mathrm{d} \Phi }{\mathrm{d} \Theta}\omega
$$

$$
e = \omega {\Phi}'
$$

Thus:
$$
\begin{bmatrix}
e_a\\ 
e_b\\ 
e_c
\end{bmatrix} = \omega_e \frac{\mathrm{d} }{\mathrm{d} \vartheta_e}\begin{bmatrix}
\Phi_a\\ 
\Phi_b\\ 
\Phi_c
\end{bmatrix} = \omega_e \begin{bmatrix}
{\Phi}'_a\\ 
{\Phi}'_b\\ 
{\Phi}'_c
\end{bmatrix}$$

As the total torque produced is the summation of the 3 phases producing torque (factoring in the pole-pair count): 

Torque equation

\$T_e = P({\Phi}'_a i_a + {\Phi}'_b i_b + {\Phi}'_b i_b)\$

Now, the Park transform can be applied to voltage, currents and flux, thus:

$$
T_e = \frac{3}{2}P({\Phi}'_\alpha i_\alpha + {\Phi}'_\beta i_\beta + {\Phi}'_0 i_0)  = \frac{3}{2}P({\Phi}'_d i_d + {\Phi}'_q i_q + {\Phi}'_0 i_0) = \frac{3}{2}P{\Phi}'_q i_q
$$ 

For a BLAC machine \${\Phi}'_q\$ is a simple sinus profile (and thus a lookup table, CORDIC... can be used) but for a BLDC an appropriate lookup table of flux vs angle is required specific for the machine being used.