The flow of polymer solutions in the high Elasticity number, El, regime in complex geometries may lead to strong viscoelastic behaviour and eventually become unstable as the Weissenberg number, Wi, is increased beyond a critical level. So far, the success of numerical simulations in predicting the highly non-linear behaviour of polymer solutions in complex flows has been limited. In this thesis, selected constitutive models are evaluated under the high El flow regime in the cross-slot and contraction benchmark flows using a numerical technique based on the finite volume method. The numerical technique is implemented within the OpenFOAM framework and thoroughly validated in the benchmark flow. A modification to the FENE dumbbell model based on the non-affine deformation of polymer solutions is proposed, which enabled the prediction of some non-linear material functions and also enhanced numerical stability, allowing a higher Wi to be attained. Asymmetric flow instability in the cross-slot flow has been studied. Time-dependent stability diagrams were constructed based on Wi and the strain, ε, both of which govern the stretching of a polymer chain. In the contraction flow, elastic instability is simulated for the first time in this geometry. Substantial time-dependent asymmetric flow patterns were predicted as seen in experiments. The effect of the contraction ratio is investigated through a stability diagram. Three-dimensional finite element simulations were also carried out to study the effect of the aspect ratio in the contraction flow of a Phan-Thien-Tanner fluid. The simulations suggest that a lip vortex mechanism is a signature for the onset of strong viscoelastic behaviour.