Even though hundreds of CubeSats have been launched, few have launched with a micropropulsion module on board. Propulsion would allow for an extended lifetime and better mission capabilities, thus greatly increasing the attractiveness of CubeSats. TU Delft is working on a type of thruster suitable for CubeSats, called the Vaporizing Liquid Microthruster (VLM). For micro-propulsion it is viable to generate sufficient thrust by using an external heat source to vaporize a liquid and expanding the vapour via a nozzle. The usage of such micropropulsion modules is still inhibited by issues in the nozzle and heating chamber. One issue is the multiphase flow occurring inside the micronozzle, which lowers the thrust and specific impulse performance of the thruster as well as the stability due to the explosive boiling phenomena. Previous research indicates that the multiphase flow in the nozzle most likely occurs inside the heating chamber and flows into the nozzle. Little is known on how exactly this multiphase flow looks like, ranging from annular flow to dispersed droplets, nor on how it is generated within the chamber. It is indicated that the boiling may lead to large pressure oscillations inside the chamber, which contributes to the instability of the thruster. These pressure oscillations occur due to multichannel interactions from nucleate boiling in microchannels, as well as explosive boiling. In this thesis, the lattice Boltzmann method (LBM) is used to simulate the complex phase change multiphase flow occurring inside the microchannels of the vaporizing liquid micro-thruster. This method has recently gained a lot of attraction in simulating microscale fluid problems. Research shows that LBM applied to multiphase flows can be an order of magnitude faster than conversional Navier-Stokes solvers. This thesis utilizes the Bhatnagar–Gross–Krook (BGK) collision operator in combination with the pseudopotential method to simulate the multiphase flow. A modified Guo forcing scheme is used to ensure thermodynamic consistency. The water is modelled using the non-ideal Peng-Robinson equation of state. The thermal problem is solved using the double distribution function method, in combination with a semi-hybrid source term solved via a multiple relaxation time (MRT) collision operator. The phase change term is derived from the local balance of entropy, using the equation of state to calculate the latent heat. The open-source lattice Boltzmann method solver OpenLB is extended to be capable of simulating phase change flow. Verification tests are performed to show that the code extension is correctly implemented. The code was released to the OpenLB, making it the first open-source lattice Boltzmann method framework capable of simulation phase change flow. In the analysis of the tool developed, the pool boiling simulations showed how nucleation, bubble formation, and departure occurs. Correct behaviour are shown in the thermodynamic consistency, Laplace pressure and wall wettability. The effect of the spurious currents on the thermal solution are analysed in more depth than found in literature. The procedure for the multichannel simulation verification and validation is given, including a novel analysis method. The novel analysis method can estimate the multiphase flow type occurring inside real thrusters without any additional equipment required in a test setup. The required experiment data are the specific impulse values at various massflow rates, and the massflow rate at which multiphase flow starts to occur. The simulations are strongly limited by the BGK collision operator, which is apparent by the large spurious currents at high density ratio multiphase flow. The spurious currents can cause simulation instabilities, but are mitigated by using a thicker diffuse interface. This mitigation results in a larger metastable phase change region, which inhibits the nucleation process. Thus, nucleation inside the microchannels did not occur. However, using a MRT collision operator should reduce the spurious currents without needing to increase the diffuse interface thickness, allowing for nucleation to take place in the microchannels.