The rheological behavior of associative polymer networks has been studied by direct Monte Carlo simulation. To describe the associations, an algorithm is used which for dilute monatomic systems reduces to the well-known mass action law. To describe the polymers, we adhere to the simplest possible model, namely hard beads connected by harmonic springs. On applying an external field the gel is stretched, and generally it starts to flow. When the how velocity is plotted against the applied stress we find Newtonian (weak gel) flow behavior if the association constant is small, and a gradual cross-over to nonlinear (strong gel) flow behavior when the association constant is larger. Furthermore, the applied flow is observed to induce gelation, hence for weak association constants we find shear thickening, whereas for strong association constants stress overshoot, thixotropic behavior and shear thinning are found. Thus all phenomena commonly observed in experiments seem to be reproduced qualitatively by the simple model. Furthermore, asa complement to existing analytic theories, the simulations are used to evaluate the chain-breaking rate function, and to measure the bond-length distribution under conditions of high shear.