How can you simulate the effect of a solvent on a polymer in an MD simulation? I'm doing a molecular dynamics simulation of polymer melts and I want to include the effect of a good/bad solvent on the melt. How can this be done?

My simulation models polymers using a bead-spring model and the interactions forces are Lennard-Jones for pair interactions and FENE for bonding forces.

Edit:

FENE stands for Finite Extensible Nonlinear Elastic. It's a log based elastic force.
The potential function used, Lennard Jones + FENE is described here (http://lammps.sandia.gov/doc/bond_fene.html)
You have generally two choices a) implicit and b) explicit solvent. Presumably, you don't have good interaction parameters for explicit solvent, since you are likely doing a generic polymer melt.

There are a variety of solvation models; LAMMPS, for instance, implements a dielectric pair style for which you can attenuate the dielectric constant of your solvent (and a few other things). Other popular approximations are generalized Born models. 

This paper uses GROMACS to simulate poly(styrene sulfonate) and poly(allylamine hydrochloride) in explicit and implicit solvent, as an example from the literature.