Heat losses in electrical devices are, usually, an unwanted effect. In order to maintain the temperature within normal operation limits, these devices are often cooled by advective fluid flow. On the other hand, many electrical devices are specifically intended to for heating up a fluid. In both cases, temperature distribution, electromagnetic fields and fluid flow are mutually dependent. Mathematically, the coupling between these three phenomena represents a complex, non-linear problem which is accessible only to numerical simulations. In this work, a methodology for the simulation of thermally coupled electromagnetic field and fluid flow is presented. The large difference in the time scales between fluid- and thermodynamics on the one hand and electromagnetic fields on the other suggests the application of a weak coupling approach. In order to obtain fluid flow and electromagnetic field solutions for complex, three dimensional devices, specialized simulation codes have to be employed. The codes used in the present work are MAFIA for the electromagnetic field and FASTEST for the fluid dynamics simulations, respectively. One difficulty related with this approach is that the simulation data reside on different grids. In addition, the fluid flow simulations are based on a Finite Volume type discretization, whereas the electromagnetic field solution is obtained in the FIT formulation. Thus, data transfer and interpolation between the two simulation codes is necessary. In particular, the heat sources generated in the electromagnetic simulation have to be interpolated into the fluid flow simulation grid. Vice versa, the temperature solution must be brought to the electromagnetic grid, in order to compute the temperature dependent material properties. These tasks have been addressed to by two different implementations of the coupling interface between the fluiddynamics and the electromagnetic field solver. The newly developed coupling tool, 3Dint, generates the geometrical relationship between the simulation grids, implements the interpolation of the coupling quantities and synchronizes program execution in the transient simulation. In addition, the consistency of the transient solution is imposed by means of a fixed-point iteration which is applied in every simulation time step. For purposes of comparison, also a coupling interface based on the commercially available library, MpCCI, is implemented. The numerical performance of the proposed algorithm is investigated for the simple case of a lightning arc propagating between two parallel electrodes. In the context of this example, the numerical convergence and the stability of the coupling approach is shown. The accuracy of the algorithm with respect to the spatial and temporal discretization is investigated. In particular, the convergence orders for the different interpolation schemes used are numerically established. The method is applied in the simulation of three large technical examples of practical importance. The first example represents a water-cooled, high-current cable which is used linear accelerators. The second is a heating module with applications in the food industry. The last example is the numerical modelling of the lightning arc as it occurs in low voltage power circuit breakers.