  We present a methodology to simulate the quantum thermodynamics of thermal
machines which are built from an interacting working medium in contact with
fermionic reservoirs at fixed temperature and chemical potential. Our method
works at finite temperature, beyond linear response and weak system-reservoir
coupling, and allows for non-quadratic interactions in the working medium. The
method uses mesoscopic reservoirs, continuously damped towards thermal
equilibrium, in order to represent continuum baths and a novel tensor network
algorithm to simulate the steady-state thermodynamics. Using the example of a
quantum-dot heat engine, we demonstrate that our technique replicates the well
known Landauer-B\"uttiker theory for efficiency and power. We then go beyond
the quadratic limit to demonstrate the capability of our method by simulating a
three-site machine with non-quadratic interactions. Remarkably, we find that
such interactions lead to power enhancement, without being detrimental to the
efficiency. Furthermore, we demonstrate the capability of our method to tackle
complex many-body systems by extracting the super-diffusive exponent for
high-temperature transport in the isotropic Heisenberg model. Finally, we
discuss transport in the gapless phase of the anisotropic Heisenberg model at
finite temperature and its connection to charge conjugation-parity, going
beyond the predictions of single-site boundary driving configurations.
