The cryogenic buffer gas beam (CBGB) is an important tool in the study of cold and ultracold molecules. While there are known techniques to enhance desired beam properties of a CBGB, such as high flux, low velocity, or reduced divergence, they have generally not undergone detailed numerical optimization. Numerical simulation of buffer gas beams is challenging, because the relevant dynamics occur in regions where the density varies by orders of magnitude, rendering typical numerical methods unreliable or intractable. Here we simulate CBGBs with a hybrid approach that combines gas dynamics methods with particle tracing. The simulations capture important properties, such as velocitiy and divergence, across an assortment of designs, including two-stage slowing cells and de Laval nozzles. This approach should therefore be a useful tool for optimizing CBGB designs across a wide range of applications.