1. A method of operating a computer system to simulate the fluid and structural behavior of a volume of the earth near a wellbore, comprising the steps of: retrieving, from a memory resource in the computer system, parameters defining properties of at least one fluid phase at each of a plurality of mesh nodes in a grid representative of a sub-surface formation in the volume to be simulated, the parameters comprising a velocity parameter for fluid at each mesh node; retrieving, from a memory resource in the computer system, parameters defining properties for each of a plurality of particles of a frangible material within the volume, the parameters for each particle comprising a velocity parameter, and a damage parameter having possible states comprising undamaged and damaged; selecting initial pressure conditions at each of the plurality of mesh nodes; then, at each of a plurality of time steps over a simulation time period, solving a system of equations comprising a state equation for fluid flow at each of the plurality of mesh nodes and a state equation for particle behavior for each of the plurality of particles; and processing the results of the solving steps to estimate fluid flow and changes in the sub-surface formation over the simulation time period; wherein the solving step comprises: for each of a first plurality of time steps separated by a first time period: setting the velocity parameter to zero for each particle having a damage parameter in the undamaged state; setting an elastic stress term to a constant for each particle having a damage parameter in the undamaged state; then solving the state equation for fluid flow at each of the plurality of mesh nodes; and solving the state equation for particle behavior for each of the plurality of particles; and after the first plurality of time steps, for each of a second plurality of time steps separated by a second time period, the second time period being shorter than the first time period: setting the velocity parameter for fluid in each mesh node to a constant value corresponding to a recent value of the velocity parameter in that mesh node; then solving the state equation for particle behavior for each of a subset of the plurality of particles within a selected portion of the volume; and then, for each particle in the subset, evaluating a damage model responsive to one or more stress values determined in the solving step, to determine the state of the damage parameter for that particle.