To reduce the climate impact of thermal engines such as heat pumps or refrigeration machines, refrigerants with a low global warming potential need to be paired with fitting lubricants. As the contamination of those liquid components influences the efficiency and lifetime of these machines, knowledge about their solubility behavior is of great interest. Molecular simulations offer mighty tools to investigate these solubilities while giving structural insight into the systems. Here the solubility behavior of CO2, R‐32, R‐1233zd(E), and R‐1234yf in PEB8, PEC4, and PEC8 is compared through the solvation free energy ∆GSolv obtained by molecular dynamics simulations. To derive ∆GSolv at low computational cost, an iterative method is used to find an optimal number and distribution of intermediate states. The resulting distributions are investigated with regard to different parameters of the employed softcore‐potential. ∆GSolv values for the different refrigerant‐lubricant pairings at different temperatures are provided, followed by a structural analysis.