Abstract Laser assisted milling (LAMill) of ceramics shows some complicated characteristics such as discontinuous chips, crack formation, propagation and coalescence. In this paper, numerical simulation is conducted to explore the machining mechanism of LAMill. The distinct element method (or discrete element method, DEM) is applied to model the microstructure of a β-type silicon nitride ceramic. Clusters are used to simulate the real grain shape of the silicon nitride ceramic and parallel bonds are employed to represent the connections between intergranular glass phase and grains. Numerical tests (compression, bending and fracture toughness tests) are performed to evaluate the macroproperties of the synthetic material, thus matching the corresponding physical properties of the real silicon nitride. Moreover, a temperature-dependent synthetic DEM specimen is created and then used in simulations of LAMill. The DEM model is validated through comparing the simulation results with the experimental ones in terms of the cutting forces and subsurface damage under different depths of cut. It is shown that the model can successfully predict the subsurface damage in LAMill.