I would like to ask the following question. Let's say we have two independent variables of x and y, and both are exponentially distributed with mean $\mu_1, \mu_2$. And let $S=x+y$. May I know whether the mean of S would be equal to $\mu_1+\mu_2$.
Thanks
This is generally true for any random variables you wish with finite means.
In other words, if $X,Y$ are random variables with $\mathbb{E}[X] < \infty$ and $\mathbb{E}[Y] < \infty$, then $\mathbb{E}[X+Y] = \mathbb{E}[X] + \mathbb{E}[Y]$.
In addition, $\mathbb{E}[cX] = c\mathbb{E}[X]$ for any constant $c \in \mathbb{R}$.
Together, these are known as the linearity of the expected value.
UPDATE
For your second question, to find the distribution of the sum, define $S = X+Y$.
Note that the event $\{X=x\}$ is a zero-probability event for any contintuous $X$ (but would be well defined for discrete $X$). To avoid this issue, define $E_x = \lim_{\epsilon \to 0} \{|X-x|<\epsilon\}$, which is a well-defined continuous analogy of $\{X=x\}$, with $\mathbb{P}[E_x] = f_X(x)$. Then $$ \begin{split} F_S(s) &= \mathbb{P}[S < s] \\ &= \mathbb{P}[X+Y < s] \\ &= \int_\mathbb{R} \mathbb{P}[X+Y < s|E_x]f_X(x) dx \\ &= \int_\mathbb{R} \mathbb{P}[Y < s-x]f_X(x) dx \\ &= \int_\mathbb{R} F_Y(s-x) f_X(x) dx, \end{split} $$ which you can now explicitly evaluate if you know the distributions of $X$ and $Y$. You get a nicer identity when differentiating that to get the pdf on the LHS: $$ f_S(s) = \int_\mathbb{R} f_Y(s-x) f_X(x) dx, $$ so the pdf of the sum is a convolution of the elementary pdfs.
It should be easy to tell that this is not likely to be of the same form as $X,Y$ even if $X,Y$ are iid. (One famous exception is the normal distribution.)
In your specific case, 2 exponentials add up to a di-gamma function...