  Taylor bubbles are a feature of the slug flow regime in gas-liquid flows in
vertical pipes. Their dynamics exhibits a number of transitions such as
symmetry-breaking in the bubble shape and wake when rising in downward-flowing
and stagnant liquids, respectively, as well as breakup in sufficiently
turbulent environments. Motivated by the need to examine the stability of a
Taylor bubble in liquids, a systematic numerical study of a steadily-moving
Taylor bubble in stagnant and flowing liquids is carried out, characterised by
a dimensionless inverse viscosity ($Nf$), and E\"{o}tv\"{o}s ($Eo$), and Froude
($Fr$) numbers based on the centreline liquid velocity, using a Galerkin
finite-element method. A boundary-fitted domain is used to examine the
dependence of the steady bubble shape on a wide range of $Nf$ and $Eo$. Our
analysis of the bubble nose and bottom curvatures shows that the intervals $Eo
= \left[ 20,30 \right)$ and $Nf=\left[60,80 \right)$ are the limits below which
surface tension and viscosity, respectively, have a strong influence on the
bubble shape. In the interval $Eo = \left(60,100 \right]$, all bubble features
studied are weakly-dependent on surface tension. This is Part I of a two-part
publication in which its companion paper (Abubakar & Matar, 2021) reports the
results of a linear stability analysis of the steady-states discussed herein.
