  Using numerical methods, we examine the dynamics of nonlinear perturbations
in the expanding time direction, under a Gowdy symmetry assumption, of FLRW
fluid solutions to the Einstein-Euler equations with a positive cosmological
constant $\Lambda>0$ and a linear equation of state $p = K\rho$ for the
parameter values $1/3<K<1$. This paper builds upon the numerical work in
\cite{Marshalloliynyk:2022} in which the simpler case of a fluid on a fixed
FLRW background spacetime was studied. The numerical results presented here
confirm that the instabilities observed in \cite{Marshalloliynyk:2022} are also
present when coupling to gravity is included as was previously conjectured in
\cite{Rendall:2004,Speck:2013}. In particular, for the full parameter range
$1/3 < K <1$, we find that the density contrast of the nonlinear perturbations
develop steep gradients near a finite number of spatial points and becomes
unbounded there at future timelike infinity. This instability is of particular
interest since it is not consistent with the standard picture for late time
expansion in cosmology.
