Abstract: The physical properties of gas hydrate-bearing coal play a crucial role in effectively evaluating the effectiveness of gas hydrate curing technology in preventing coal and gas outburst. To address the problems of inaccurate reproduction of the real deformation characteristics of the specimen and low computational efficiency in the numerical simulation of the physical properties of gas hydrate-bearing coal, the coupling method of FLAC^3D and PFC^3D is considered to be applied to the establishment of boundary conditions. Therefore, a three-axis discrete element numerical model is established by taking the indoor three-axis test of a gas hydrate-bearing coal with a hydrate saturation of 31% as an example. Considering the deformation characteristics of the rubber sleeve wrapping the coal samples, the FLAC^3D-PFC^3D numerical coupling method is used to create a lateral flexible boundary, i.e., FLAC^3D is used to simulate the macroscopic deformation characteristics of the rubber sleeve, which ensures the full evolution of the shear zone of the numerical model, and PFC^3D is used to simulate the micromechanical characteristics of the coal. The coupling is established through the wall control by the boundary at the interface of the contact between the discrete and the continuous ones. In order to calculate the volumetric strain of the model for the analysis of the model shear expansion characteristics, a three-dimensional flexible boundary volumetric strain and radial strain calculation method is proposed. In order to meet the three-dimensional numerical simulation working condition verification, a triaxial test platform of gas hydrate-bearing coal was constructed, and by comparing the strength and damage characteristics of the gas hydrate-bearing coal during the triaxial compression process, the detailed parameters of the numerical model were determined, and the reliability of the model was verified. The results show that the proposed model can not only accurately reflect the strength and damage characteristics of the gas hydrate-bearing coal during triaxial compression, but also achieve the prediction of its shear expansion characteristics and shear zone development at the fine-scale. The numerical model can be used as a necessary supplement to the triaxial test of gas hydrate-bearing coal. This study presents a novel approach that enables realistic and computationally efficient triaxial simulation of gas hydrate-bearing coal, thereby providing theoretical support for its mechanical enhancement and the development of a mesomechanics-informed constitutive model.