The volume penalization method (VP), a type of Fictitious Domain Method, is a widely used technique for solving partial differential equations (PDEs) in complex domains. Its appli- cations span various fields, from fluid-structure interactions like wave energy converters, bird and insect flight, fish swimming, and cardiovascular flows, to phase change applications such as glacier melting and additive manufacturing processes. This thesis presents robust and adap- tive VP techniques for simulating non-isothermal phase-changing flows, as well as isothermal multiphase fluid-structure interaction problems. Using the numerically constructed flux-forcing functions for arbitrarily complex boundaries, we extend the flux-based volume penalization (VP) method to handle more general boundary conditions, including spatially varying inhomogeneous Neumann and Robin boundary conditions. Several two- and three-dimensional test examples, including flux-driven thermal convection in a concentric annular domain, are considered to assess the spatial accuracy of the numerical solutions. In addition, we propose a projection method-based preconditioning strategy for solving VP incompressible and low-Mach Navier- Stokes equations. The solver converges faster as the penalty coefficient decreases, contrary to prior experience. The developed preconditioning strategy is used in a novel low Mach enthalpy method to solve solidification and melting problems with variable thermophysical properties, including density. The proposed method captures the density change-induced flow during phase change material (PCM) melting and solidification. A gas phase is also incorporated and coupled to the solid-liquid PCM region in this formulation. The new low Mach enthalpy method is validated against analytical solutions for a PCM undergoing a large density change during its phase transition. Furthermore, we propose a set of simple sanity checks to serve as benchmarks for evaluating computational fluid dynamics (CFD) algorithms that aim to capture the volume change effects of PCMs. Adaptive mesh refinement is employed to achieve fine grid resolution in domains requiring more accuracy, such as PCM-gas and liquid-solid interfaces.