Consider this Lagrangian density $$\mathcal{L}(x)=-\frac{1}{2}[a\partial_{\mu}A^{\nu}\partial^{\mu}A_{\nu}+b\partial_{\mu}A^{\nu}\partial_{\nu}A^{\mu}+c(\partial_{\mu}A^{\mu})^2+dA_{\mu}A^{\mu}]$$
The Euler-Lagrange e.q is: $$a\partial^{\nu}\partial_{\nu}A^{\mu}+(b+c)\partial^{\mu}\partial^{\nu}A_{\nu}-dA^{\mu}=0.$$
I want to know the detail reasoning process.the index operation make me confusion.
It's just a manner of being careful with the indices. Thake the Euler-Lagrange equation in covariant formulation
$$\partial_\rho \frac{\delta \mathcal{L}}{\delta \partial_\rho A_\sigma}-\frac{\delta\mathcal{L}}{\delta A_\sigma}=0$$ and let us start from the second bit. Clearly, the only part in the lagrangian which does not depend on the derivative of the field is the $d$-term which gives simply
$$\frac{\delta}{\delta A_\sigma}\frac{d}{2}A_\mu A^\mu = d\delta^{\sigma}_\mu A^\mu = dA^\sigma$$ where the factor $1/2$ is cancelled by the factor $2$ of the derivative. As you can see the free index on the right is $\sigma$ as well as on the left, which is what we want.
The other bit is the more involved one. Take it bit by bit
$$\frac{\delta}{\delta\partial_\rho A_\sigma}\frac{a}{2}\partial_\mu A_\nu \partial^\mu A^\nu = a\delta_\mu^\rho\delta_\nu^\sigma \partial^\mu A^\nu = a\partial^\rho A^\sigma$$ where one Kroneker delta comes from the derivative and the other from the field. Again the factor $1/2$ is cancelled by the derivative since it is a quadratic term. Then,
$$\frac{\delta}{\delta\partial_\rho A_\sigma} \frac{b}{2}\partial_\mu A_\nu\partial^\nu A^\mu = b\delta_\mu^\rho\delta_\nu^\sigma \partial^\nu A^\mu = b\partial^\sigma A^\rho$$
notice that the index are swapped with respect to the previous one. The last bit remaining is
$$\frac{\delta}{\delta\partial_\rho A_\sigma}\frac{c}{2}(\partial_\mu A^\mu)^2 = c\delta_\mu^\rho \delta^\mu_\sigma\partial_\mu A^\mu = c\partial^\sigma A^\rho.$$
Now it is just a manner of putting all together taking the other derivative from the first bit of the EL-equations
$$\partial_\rho [a\partial^\rho A^\sigma + (b+c)\partial^\sigma A^\rho] = a\partial^\rho\partial_\rho A^\sigma+ (b+c)\partial_\rho\partial^\sigma A^\rho.$$
Putting it all back together being wary of the signs and you get the desired result. You can clearly rename the indices however you want.
A second way to derive the result is through the use of the $\delta$-symbol, also called "variation". It usually is a bit easier with managing indices, as you do not have to manage Kronecker-deltas (and it also clarifies the relation between symmetries and equations of motion).
The basic properties needed are, that $\delta$ is linear $$ \delta(A^\mu + B^\mu) = \delta A^\mu + \delta B^\mu $$ it commutes with the derivative $$ \delta \partial_\mu A^\nu = \partial_\mu \delta A^\nu $$ and it is a derivation (Leibnitz rule) $$ \delta (A^\mu B^\nu) = (\delta A^\mu) B^\nu + A^\mu \delta B^\nu. $$ Applying these rules in the action and integrating by parts, so that no derivative acts on $\delta$, you can read off the equations of motion.
Writing this abstractly is not very helpful, so let's check out your example:
By linearity we calculate each term separately:
$$ \int d^Dx \delta( a \partial_\mu A^\nu \partial^\mu A_\nu) = \int d^Dx a \left[ (\partial_\mu \delta A^\nu) \partial^\mu A_\nu + \partial_\mu A^\nu (\partial^\mu \delta A_\nu) \right] \\ = \int d^Dx a \left[- (\partial_\mu \partial^\mu A_\nu) \delta A^\nu - (\partial^\mu\partial_\mu A^\nu) \delta A_\nu \right] \\ =\int d^Dx [(-2a) \partial^\mu\partial_\mu A^\nu]\delta A_\nu, $$ where in the second line we integrated by parts (discarding the boundary terms) and in the third, we collected terms and raised/lowered indices such that $\delta A_\nu$ is isolated.
You can already see that the term in brackets is exactly the term you get in the equations of motion, but to get there, you do not have to bother with as many Kronecker-deltas. Repeating the same steps with the other terms of the Lagrangian will get the desired equations of motion.
