The nonlinear, resistive, 3D magnetohydrodynamic equations are solved numerically to demonstrate the possibility of forming and sustaining a spheromak by forcing tangential flows at the plasma boundary. The method can by explained in terms of helicity injection and differs from other helicity injection methods employed in the past. Several features which were also observed in previous dc helicity injection experiments are identified and discussed.Spheromak plasmas have been formed using several different techniques [1,2]. The existence of multiple formation methods clearly shows that the spheromak is a preferred (lowest energy) state toward which a magnetohydrodynamic (MHD) system naturally evolves when the appropriate boundary conditions are imposed. The physical process responsible for the formation is magnetic relaxation: on the time scale of MHD instabilities the plasma relaxes to the minimum energy state compatible with its magnetic helicity content (which remains approximately constant) [3]. Once formed, the spheromak will decay in a resistive time scale because resistivity does dissipate magnetic helicity. For this reason, the sustainment of the configuration during times longer than the resistive decay time requires some helicity injection method. Since relaxation operates on a shorter time scale, the spheromak configuration is maintained regardless of the details of the specific helicity injection mechanism. Some particular examples are the coaxial helicity injection method (CHI) [4][5][6], the merging of helicitycarrying filaments (MHF) [7] and the helicity injected torus with steady inductive helicity injection (HIT-SI) [8].In this Letter we report the first evidence coming from nonlinear, resistive, 3D MHD numerical simulations that demonstrate the possibility of forming and sustaining a spheromak by forcing tangential flows at the plasma boundary. The method can by explained in terms of helicity injection and differs from other helicity injection methods employed in the past (CHI, MHF and HIT-SI).An enhanced helicity injection mode was recently reported in spheromak experiments with large plasma rotation [9]. Althought not analyzed in terms of boundary flows, this observation could support the feasibility of the mechanism proposed and studied in this Letter.We model the plasma using the resistive MHD equations with finite viscosity and zero β. The evolution equations for u and B are:where E = −u × B + ηJ and Π = (∇u + ∇u T ) − (2/3) × (∇·u) (see Ref.[10] for further details). These equations are normalized with a (chamber radius), ψ G (imposed flux) and c A (Alfvèn speed). In addition, B and J are scaled with √ µ 0 . Time is expressed in units of the Alfvèn time τ A = a/c A . The normalized resistivity η and the kinematic viscosity ν are set to 5 × 10 −5 . With these values the resistive time is τ r ∼ 800 (defined as in Ref.[11]). The resulting system is solved with the Versatile Advection Code [12].The domain is a cylinder of elongation h/a = 1, with perfectly conducting wall conditions (B · n = 0 ...