understanding discrete-time convolution in LTI systems
You know that $h(n)$ is the system's response to the input $x(n)=\delta(n)$. Because of time-invariance, the response to $x(n)=\delta(n-k)$ is given by $h(n-k)$. If you write the input signal as

$$x(n)=\sum_kx(k)\delta(n-k)\tag{1}$$

then, due to linearity, you can write the output signal as the sum of the responses to each individual term $x(k)\delta(n-k)$ in (1), which is $x(k)h(n-k)$. So you get the convolution sum

$$y(n)=\sum_k x(k)h(n-k)\tag{2}$$

To compute (2) for a given value of $n$ you need to multiply $x(k)$ with $h(n-k)$ for all relevant values of $k$. Since you have $h(n\color{red}{-k})$ in (2) (i.e. with a minus sign before $k$) you need to flip the impulse response (shifted by the current value of $n$), multiply it with $x(k)$ and compute the sum.