A couple things:

BN applies an affine transform. So you want to set the affine scaling weight to 1 before getting Y, i.e. B2.weight.data.fill_(1).
Use biased version of variance.
BN normalizes for data within each channel. So you should calculate Y's variance for each channel, rather than for all data in Y.

Altogether, this gives correct result:
X3 = torch.randn(150,20,100) * 2 + 4
X3 = Variable(X3)
B2 = nn.BatchNorm1d(20)
B2.weight.data.fill_(1)
Y = B2(X3)
Y = Y.transpose(0, 1).contiguous().view(20, -1)  # put data for each channel in the second dimension
print(Y.var(dim=-1, unbiased=False))