z-test vs t-test under assumptions A traditional t-test to compare two means assumes
independent random variables, $X_1,...,X_n \sim  n(\mu_1,\sigma_1 ^2)$
and $Y_1, ...,Y_m \sim n(\mu_2,\sigma_2 ^2)$.
One usually assumes normal distributions with $\sigma_1$ = $\sigma_2$. Then the statistics
$$\frac{\bar{X}-\bar{Y}}{S\sqrt{1/n+1/m}}$$
is used and it's observed value is compared to $t_{n+m-2,1-\alpha/2}$ where $S^2$ is the pooled estimate of variance, valid under both the null and alternative,
given the assumption of equal variances.
If the variances are known one can use
$$\frac{\bar{X}-\bar{Y}}{\sqrt{\sigma_1^2/n+\sigma_2^2/m}}$$
and compare with a $z_{1-\alpha/2}$ test.
If we assume $n=m$ and for example $\sigma_1^2=1$ and $\sigma_2^2=2$ and investigate using simulation which test is better. How would we do that. It seems wise to check level and powers but how would one do that?
 A: Without loss of generality, you can assume that the mean of the first group, $\mu_1$ is zero, and the mean of the second group is $\mu_2 = \mu_1+\delta = \delta$.
You'll need to assume some particular sample sizes (specific m and n), or do this over a grid of such values.
Let's take level: here $\delta = 0$. 
So you know the parameters of both distributions. Basically:
Repeat nsim times:
    simulate sample 1
    simulate sample 2
    compute z statistic 
    if z rejects, increment the count of z rejections
    compute t statistic
    if t rejects, increment the count of t rejections

Compute rejection rates for t and z

The case for power is the same, expect you do it for some set of values of $\delta$.
A common thing to do at the end is plot the rejection rates against delta, so you can get an idea what the power curve looks like.

Edit: Okay, you work in R; that would have been handy to know.
In R you probably want to use replicate rather than a loop.
Since you would do this for a set of delta values, you'd probably wrap that in a call to sapply.
[and then, if you then did that over various values of n, you might even put that in a loop]

Here's an example of some R code for estimating level and power at a fixed n
nsim <- 1000
sigxsq <- 1
sigysq <- 2
n <- 10
m <- 10
sigx <- sqrt(sigxsq)
sigy <- sqrt(sigysq)

testasample <- function(n,m,delta=0,mux=0,muy=mux+delta,sigx=1,sigy=1,alpha=0.10,
                 zs=sqrt(sigx^2/n+sigy^2/m)){
  x <- rnorm(n,mux,sigx)
  y <- rnorm(m,muy,sigy)
  tt <- t.test(x,y,var.equal=TRUE)
  tp <- tt$p.value
  zp <- as.numeric(2*pnorm(-abs(diff(tt$estimate)/zs)))
  c(tp,zp) < alpha
}

del <- c(0,0.2,0.5,0.7,1,1.5,2)
res <- sapply(del,function(d){
         rowMeans(replicate(nsim,testasample(n,m,delta=d,sigx=sigx,sigy=sigy)))})

res <- cbind(del,t(res))
colnames(res)<-c("delta","trej","zrej")

It does a little unnecessary calculation, but that will hardly change the overall speed so I didn't bother with it.
Here's a plot that I generated after running the above code:
 plot(trej~delta,res,pch="t",col=4,ylim=c(0,1),ylab="rejection rate")
 points(zrej~delta,res,pch="z",col=2)
 abline(h=c(0,1,.1),col=8,lty=c(1,1,3))


Note that you can compute standard errors (and intervals) for these estimated rejection rates. 
A: n <- 1000
h <- 30
rt <- rep(0,h)
rz <- rep(0,h)
for (j in 2:h) {
    for (i in 1:n) {
        p <- j
        tstat <- 0
        zstat <- 0
        Xvec <- rnorm(p, mean = 0, sd = 1)
        Yvec <- rnorm(p,mean = 0, sd = sqrt(2))
        xsum <- 0
        ysum <- 0
    for (k in 1:p) { 
        xsum <- (Xvec[k]-mean(Xvec))^2+xsum 
        ysum <- (Yvec[k]-mean(Yvec))^2+ysum 
    }
        Ssq <- (xsum+ysum)/(2*p-2)
        tstat <- (mean(Xvec)-mean(Yvec))/(sqrt(Ssq)*sqrt(2/p))
        if (abs(tstat) > qt(0.95,2*p-2)) {
            rt[p] <- rt[p]+1
        }
        zstat <- (mean(Xvec)-mean(Yvec))/(sqrt(3/p))
        if (abs(zstat) > qnorm(0.95, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)) { 
            rz[p] <- rz[p]+1
        }
    }
}
rt <- rt/n
rz <- rz/n

