Large bubbles always undergo asymmetric paths due to the wake instability, which was rarely considered in the simulation of bubbly flow in complex flow systems such as bubble columns. Here, we propose a bubble path oscillation model considering the zigzag lateral motion caused by the periodic vortex shedding, which is added into the rectilinear bubble motion model in the Lagrange frame. This model is further adopted to simulate the flow of 4 mm bubbles with 680 < Re < 1000 in bubble columns, with the liquid flow predicted by the partially averaged Navier-Stokes model in the Euler frame. A parameter analysis is conducted by comparing the simulation results with the existing experimental data to determine the optimal oscillation amplitude of the path. The effects of oscillation amplitude on the bubble distribution, gas holdup, and flow field are further discussed. Our findings may contribute to the fundamental understanding and accurate simulation of bubbly flow.