The bubble rise behavior in viscoelastic media is analyzed numerically with CFD. Three different constitutive models, Giesekus, linear and exponential Phan-Thien-Tanner, are used to evaluate three different biopolymer solutions. The terminal rise velocity over a range of bubble sizes is validated against experimental data. The local velocity fields are compared with respect to the shape and onset of the negative wake. Furthermore, the normal and shear components of the stress fields, transformed according to the local flow direction, are given. The simulations are performed with a volume of fluid solver in OpenFOAM.