Why is circular convolution used in DSP? Why not linear convolution?
Given a discrete time LTI system with impulse response $h[n]$, one can compute its response to any input $x[n]$ by a convolution sum: $$y[n] = x[n] \star h[n] =  \sum_{k=-\infty}^{\infty} {h[k]x[n-k]} \tag{1}$$

Without anything further stated, above definition is for the linear convolution (aperiodic convolution) between $h[n]$ and $x[n]$, which are aperiodic discrete-time sequences of possibly infinite length, unless otherwise stated.  This is different from a circular convolution which is between two periodic sequences of period $N$, and computed over a single period.

You can compute a linear convolution in time domain by Eq.1, or in frequency domain using the following DTFT (discrete-time Fourier transform) property: 
$$y[n] = x[n]\star h[n] \implies  Y(e^{j\omega}) = X(e^{j\omega}) H(e^{j\omega}) \tag{2}$$

DTFT is naturally related with the linear convolution, as it deals with theoretically existing aperiodic sequences that may extend from $-\infty$ to $\infty$ reflected in its limits of the defining sum:
$$X(e^{j \omega}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n} \tag{3}$$

When you want to make a calculation by hand, using symbolic expressions for signals, such as $x[n]=a^n u[n]$ and $h[n]=b^n u[n]$, you can calculate the results in time or frequency domains as outlined above.  

Also, when you want to compute the same result by using a computer, you can use the time domain approach based on an LCCDE recursion (for IIR systems) or direct finite convolution sum (for FIR systems), BUT the frequency domain approach won't work with DTFT; as it's mainly a tool used to develop the mathematics of signals &amp; systems theory, and it's not suitable for digital computer implementations, as its variable $\omega$ is a real continuous number. 

What is used instead is the DFT (discrete Fourier transform) defined as $$ X[k] = X(e^{j\omega})|_{\omega = \frac {2\pi k}{N} } \tag{4}$$

where $k=0,1,...,N-1$ and $N$ is the length of the DFT, which's then called as an N-point DFT of signal $x[n]$. 

Eq.4 implies that the DFT sequence $X[k]$ is obtained as the uniform samples of the DTFT $X(e^{j\omega})$, which is a periodic function, hence DFT $X[k]$ is also periodic but we only consider its first period from $k=0$ to $N-1$.

Since DFT sequences are inherently periodic, then their convolutions will also be periodic (circular). Therefore, whereas a linear convolution between aperiodic signals $x[n]$ and $y[n]$ is implied by the I-DTFT expression $$y[n] =\mathcal{I-DTFT} \{ X(e^{j\omega})H(e^{j\omega}) \}$$ instead a circular convolution between two periodic sequences is implied by the I-DFT expression $$\tilde{y}[n] = \mathcal{I-DFT} \{ X[k]H[k] \}$$

So, given that we want to compute a linear convolution between two aperiodic sequences $x[n]$ and $h[n]$ of lengths $L_x$ and $L_h$ respectively, using frequency domain by their $N$ point DFTs, $X[k]$ and $H[k]$, we actually have to compute a circular convolution between the periodic extensions of the signals $\tilde{x}[n]$ and $\tilde{h}[n]$ of periods $N$. 

The key is in chosing a proper length $N$ of the DFTs, which must be long enough to avoid any time domain aliasing of the sequence $\tilde{y}[n]$, returned by the IDFT computation: $$ \tilde{y}[n] = \sum_{r=-\infty}^{\infty} y[n-rN] \tag{5}$$ 

where $y[n]$ is the result of the linear convolution that would be returned by the theoretical inverse DTFT and $\tilde{y}[n]$ is the periodic result of the periodic (circular) convolution implied by the inverse DFT. 

Note that if any one of the signals are of infinite length, then it's NOT possible to compute their linear convolution using the DFT approach, as $N$ would go to infinity, practically impossible.
The implementation of a linear convolution via DFT then has the following steps:


Choose N according to the following criteria: $$N \geq L_x + L_h -1$$ which guarantees an alias-free reconstruction of the inverse signal $y[n]$ from its DFT $Y[k]$ of the computed circular convolution via $X[k]H[k]$. 
Compute N-point DFTs $X[k]$ and $H[k]$ of $x[n]$ and $h[n]$.
Compute $Y[k] = X[k]H[k]$
Compute N-point inverse DFT of $Y[k]$ to produce the output $y[n]$