In order to cope with the decarbonization challenge faced by many countries, fusion is one of the few alternatives to fossil fuels for the production of electricity. Two devices invented in the middle of the previous century have emerged as the most promising to confine a plasma at temperatures larger than $10^7\,$K for sufficiently long times and at high densities: the tokamak and the stellarator. In these machines, fusion reactions release energy that is envisaged to be used for electricity production. Due to the lack of toroidal symmetry, the stellarator has been plagued since its beginning by high levels of transport, preventing a satisfactory confinement of the plasma. However, with the discovery of quasi-symmetry and the advent of high-performance computers, stellarators have been optimized, and recent results of the optimized stellarator Wendelstein 7-X show that it is possible to reduce the stellarator transport level to that of tokamaks. This has lead to an increasing interest in this machine, something that is also reflected in the number of stellarator private companies that appeared in the last few years. While the interest in the stellarator as a viable device for the production of energy has increased, plasma turbulence in the boundary region of this machine is still poorly understood. The main goal of this thesis is to study the turbulent dynamics of the plasma in stellarators through global, fluid, flux-driven simulations. We use the GBS code that solves the drift-reduced Braginskii equations, valid in the limit of high collisionality, typical of the boundary of magnetic fusion devices. The code is extended to accommodate non-axisymmetric magnetic fields, and as a first application we perform a set of simulations of a diverted tokamak with applied 3D perturbations (RMPs). We show that increasing the perturbation amplitude reduces the peak heat flux towards the divertor targets, hence suggesting that RMPs could be used in the future to help mitigating the heat exhaust issue. The first global fluid simulation of a stellarator is then performed. We construct a vacuum magnetic field using the theory of the Dommaschk potentials, and we taylor the simulation domain to create an island divertor stellarator configuration. The simulations show that the plasma dynamics is different from typical tokamak simulations. In particular, a large, coherent mode with low poloidal mode number is responsible for most of the transport, contrasting with the broad-band, blobby turbulence often seen in tokamaks. Moreover, the radial extension of the mode is comparable to the poloidal one, in contrast to previous estimates in tokamak geometry that predict a larger extension in the radial direction. We find out that one of the key parameters leading to these differences --- coherent vs. broad-band; small-scale vs. large-scale structures --- is the magnetic shear. In fact, simulations of tokamaks with very low magnetic shear share many features with the simulated stellarators (which are shearless). We further explore different stellarator configurations to conclude that not only the shear but also the ellipticity and the torsion have an effect on turbulence. A validation of the GBS code is also presented, by comparing results of simulations of the TJ-K stellarator against experimental data. As in the experiments, transport in the simulations is dominated by a large, coherent mode, with mode numbers in good agreement.