Roche limits with a hint of internal bonds.

Hello everyone,
I have created a program that plots the orbit of a comet around our sun, as it orbits, a ratio between internal binding forces and tidal forces is calculated (this produces a Roche limit). I am using this equation to obtain the internal bonding force:
F = (G * M * Mu) / (r ^ 2) 'I have taken Mu = 1 kg.
Where r is the radius of the comet, M is the mass of the comet and the Mu is the mass of an object on the surface. I am wondering how I could include an extra component in this equation to make it more realistic. I want to take into account the extra binding forces caused by a shell of ice around the surface of the comet. I am not sure how to go about doing this. A shell of ice will result in a higher tidal force needed to break up the comet (Roche limit will be less). Any help on this will be very helpfull... lol.
Many thanks!!!!

(Note: I will also apply this extra ice shell component to an inner liquid center, I hope to end up with four Roche limit models,
1. Rubble-pile model (What I have already done)
2. Rubble-pile with ice surface
3. liquid sphere model (got data on this also)
4. liquid sphere with ice shell model)