Let's assume three things. First, that the "appropriate boundary condition" is $u=0$ on $\partial \Omega$. Second that we write $L_k$ for the operator defined with $A_k,b_k,c_k, r_k$. Third, we need to avoid the possibility that $0$ is in the spectrum of the operator $L_k$ for any $k$ so that we can get estimates for solutions to $L_k u =f$ solely in terms of $f$. Let's assume this.
Now suppose that $L_k u_k = f \in H^{-1}$ for each $k \ge 1$. Then the usual $H^1(\Omega)$ a priori estimate for weak solutions shows that
$$
\| u_k \|_{H^1_0} \le C(\|A_k,b_k,c_k,r_k\|_{L^\infty},\Omega) \|f\|_{H^{-1}}
$$
where $C$ can be shown to increase with respect to the $L^\infty$ part. Consequently, if $(A_k,b_k,c_k,r_k) \to (A,b,c,r)$ in $L^\infty$ and $A$ remains uniformly elliptic on $\Omega$, then
$$
\sup_{k\ge 1}\| u_k \|_{H^1_0} \le \tilde{C}(\|A,b,c,r\|_{L^\infty},\Omega) \|f\|_{H^{-1}}
$$
for some other constant $\tilde{C}$.
We can now use this bound together with weak compactness and Rellich's theorem to extract a subsequence (still denoted $u_k$ for the sake of simplicity) such that $u_k \rightharpoonup u$ weakly in $H^1_0$ and $u_k \to u$ strongly in $L^2$. Then
$$
\int_\Omega A_k \nabla u_k \cdot \nabla v - u_k b_k \cdot \nabla v + c_k \cdot \nabla u_k v + r_k u_k v = <f,v>
$$
for all $v \in H^1_0(\Omega)$, and so we can send $k \to \infty$ and use the above convergence results to deduce that
$$
\int_\Omega A \nabla u \cdot \nabla v - u b \cdot \nabla v + c \cdot \nabla u v + r u v = <f,v>
$$
for all $v \in H^1_0$, and hence $u$ solves
$$
-\text{div}(A \nabla u) + \text{div}(bu) + c\cdot \nabla u + ru =f.
$$
Note: without the condition of avoiding $0$ we would have to include the $L^2$ norm of $u_k$ on the right side of our basic a priori estimate, and so we would have to have boundedness in $L^2$ in order to extract limits.
