Summary A novel smoothing particle hydrodynamics (SPH)‐like Lagrangian meshfree method, named as Lagrangian gradient smoothing method (L‐GSM), has been proposed to avoid the “tensile instability” issue in SPH simulation by replacing the SPH particle‐summation gradient approximation technique with a local grid‐based GSM gradient smoothing operator. The L‐GSM model has been proven effective and efficient when applied to a wide range of large deformation problems for fluids and flowing solids in two‐dimensional case. In this study, a three‐dimensional (3D) L‐GSM numerical framework is proposed for simulating large deformation problems with the existence of free surfaces through developing a widely adaptable 3D gradient smoothing domain (GSD) constructing algorithm. It includes three key novel ingredients: (i) the localized GSD based on an efficient distance‐oriented particle‐searching algorithm enabling both easy implementation and efficient computation; (ii) a novel algorithm for constructing 3D GSD to guarantee the effectiveness of the 3D GSM gradient operator adaptable to any extreme cases; (iii) a robust normalized 3D GSM gradient operator formulation that can restore the accuracy of gradient approximation even on boundary interface. The effectiveness of the proposed 3D GSD‐constructing algorithm is first verified under various distribution conditions of particles. The accuracy of the proposed adaptable 3D GSM gradient algorithm is then examined through conducting a series of numerical experiments with different spacing ratios. Finally, the 3D L‐GSM numerical framework is applied to solve a practical problem of free surface flows with large deformation: collapse of a soil column. The results reveal that the present adaptable 3D L‐GSM numerical framework can effectively handle the large deformation problems, like flowing solids, with a constantly changing arbitrary free surface profile.