We investigate the error of Euler-Maruyama approximate solution of reflecting　fractional Brownian motion with some boundary conditions on multi dimensional domains(so-called Skorohod SDE) using the penalty method. The fractional Brownian motion is applied in many fields including mathematical finance. Forthe purpose of observing numerical data by computer simulation of Skorohod SDE we define discretized approximate solution of　the SDE and obtain its error estimation. The penalty method is separated into two steps. The first step is to approximate reflecting fractional Brownian motions by usual fractional　Brownian motion with a drift term whose coefficients are gradient type in pathwise sense. The second step is to approximate the fractional Brownian motion with such drift term by Euler-Maruyama scheme.