In this thesis, we use lattice QCD to study a part of the QCD phase diagram, specifically the QCD phase transition at mu=0, where the QCD matter changes from hadron gas to quark-gluon plasma (QGP) with increasing temperature. This phase transition takes place as a crossover, but when theoretically changing the masses of the quarks, the order of the phase transition changes as well. We focus on the region of heavy quark masses with Nf=2 flavours, where we investigate the critical quark mass at the second order phase transition in the form of a Z2 point between the first-order and the crossover region. The first-order region is positioned at infinitely heavy quarks. As the quark masses decrease, the associated Z3 centre symmetry breaks explicitly, causing the first-order phase transition to weaken until it turns into the Z2 point and finally into a crossover. We study this Z2 point using simulations at Nf=2 and lattices of the sizes Nt = {6, 8, 10, 12}, partially building on previous work, in which the simulations for Nt = {6, 8, 10} were started. The simulations for Nt=12 are not finished yet though, but we were able to draw some preliminary conclusions. These simulations are run on GPUs and CPUs, using the codes Cl2QCD and open-QCD-FASTSUM, respectively. Afterwards, the data goes through a first analysis step in the form of the Python program PLASMA, preparing it for the two techniques we use to analyse the nature of the phase transition. As a first, reliable analysis method, we perform a finite size scaling analysis of the data to find the location of the Z2 point. Since we are using lattice QCD, performing a continuum extrapolation is necessary to reach the continuum result. In regard to this, the finite size scaling analysis method is hampered by the excessive amount of simulated data that is needed regarding statistics and the total number of simulations, which is why this thesis is only an intermediate step towards the continuum limit. This also leads to the second analysis technique we explore in this thesis. We start to design a Landau theory which describes the phase boundary for heavy masses at Nf=2 based on the simulated data. We develop a Landau functional for every Nt we have simulation data for. Albeit the results are not at the same precision as the ones from the finite size scaling analysis, we are able to reproduce the position of the Z2 point for every Nt. Even though we are not able to take a continuum extrapolation right now, after more development takes place in future works, this approach might, in the long run, lead to a continuum result that won't need as many simulations as the finite size scaling analysis.