You correctly wrote down the pooled estimator:
$$ \bar{U} = \frac{1}{m} \sum_{i=1}^m U_i$$
Where $U_i$ represents the analytic results from the $i$-th imputed dataset. Normally, analytic results have some normal approximating distribution from which we draw inference or create confidence bounds. This is mainly done using the mean value ($U_i$) and its standard error. T-tests, linear regressions, logistic regressions, and basically most analyses can be adequately summarized in terms of that value $U_i$ and its standard error $\text{se}(U_i)$.
Rubin's Rules uses the law of total variance to write down the variance as the sum of a between and within imputation variance:
$$\text{var}(\bar{U}) = E[\text{var}(\bar{U}|U_i)] + \mbox{var}\left(E[\bar{U}|U_i]\right)$$
The first term is the within-variance such that $E[\text{var}(\bar{U}|U_i) = \frac{1}{m}\sum_{i=1}^m V_i$ where $V_i$ is the variance of the analysis result from the $i$-th complete or imputed dataset. The latter term is the between-imputation variance: $ \mbox{var}\left(E[\bar{U}|U_i]\right) = \frac{M+1}{M-1} \sum_{i=1}^m\left(U_i - \bar{U}\right)^2$. I've never quite grasped the DF correction here, but this is basically the accepted approach.
Anyway, since the recommended number of imputations is small (Rubin suggests as few as 5), it is typically possible to compute this number by hand fitting each analysis. A by-hand example is listed below:
require(mice)
set.seed(123)
nhimp <- mice(nhanes)
sapply(1:5, function(i) {
fit <- lm(chl ~ bmi, data=complete(nhimp, i))
print(c('coef'=coef(fit)[2], 'var'=vcov(fit)[2, 2]))
})
Gives the following output:
coef.bmi var
2.123417 4.542842
3.295818 3.801829
2.866338 3.034773
1.994418 4.124130
3.153911 3.531536
So the within variance is the average of the imputation specific point estimate variances: 3.8 (average of second column). The between variance is 0.35 variance of the first column). Using the DF correction we get variance 4.23. This agrees with the pool command given in the mice package.
> fit <- with(data=nhimp,exp=lm(chl~bmi))
> summary(pool(fit))
est se t df Pr(>|t|) lo 95 hi 95 nmis fmi lambda
(Intercept) 119.03466 54.716451 2.175482 19.12944 0.04233303 4.564233 233.505080 NA 0.1580941 0.07444487
bmi 2.68678 2.057294 1.305978 18.21792 0.20781073 -1.631731 7.005291 9 0.1853028 0.10051760
which shows the SE = 2.057 for the model coefficient, (Variance = SE**2 = 4.23).
I fail to see how increasing the number of imputed datasets creates any particular issue. If you cannot supply an example of the error, I don't know how to be more helpful. But by-hand combination is certain to accommodate a variety of modeling strategies.
This paper discusses other ways that the law of total variance can derive other estimates of the variance of the pooled estimate. In particular, the authors point out (correctly) that the necessary assumption for Rubin's Rules is not normality of the point estimates but something called congeniality. WRT normality, most point estimates that come from regression models have rapid convergence under the central limit theorem, and the bootstrap can show you this.
