Here is an abstract nonsense proof:
The vector
$$\vec T(\vec{a},\vec{b},\vec{c}):=\vec{a} \times (\vec{b} \times \vec{c}) + \vec{b} \times (\vec{c} \times \vec{a}) + \vec{c} \times (\vec{a} \times \vec{b})$$
is easily seen to be a skew trilinear function of the three vector variables $\vec{a}$, $\vec{b}$, $\vec{c}$. It follows that its coordinates $T_i$ are three real-valued such functions, whence are multiples of the determinant function (so-called triple vector product) $[\vec{a},\vec{b},\vec{c}]$ with factors $\lambda_i$ independent of $\vec{a}$, $\vec{b}$, $\vec{c}$. Putting $\vec{p}:=(\lambda_1,\lambda_2,\lambda_3)$ we therefore have
$$\vec T(\vec{a},\vec{b},\vec{c})=[\vec{a},\vec{b},\vec{c}]\>\vec{p}$$
with a universal vector $\vec{p}$. This only makes sense when $\vec{p}=\vec{0}$.
In the cases when $A,B$ and $A+B$ are all invertible, the inverse of the equality is:
$$B^{-1}(A+B)A^{-1}=A^{-1}(A+B)B^{-1}$$
But simple calculation shows the left side is $B^{-1}+A^{-1},$ and the right side is $A^{-1}+B^{-1}.$
The general case can be shown by noticing that the set of pairs $(A,B)$ such that $A,B,A+B$ are invertible is dense in the set of pairs $A,B$ such that $A+B$ is invertible. Since the function is continous, that would finish it.
We can do this by taking an arbitrary $A,B$ and then replacing it with $A+\lambda I, B-\lambda I,$ where $\lambda$ is a positive value with smaller magnitude than any of the non-zero eigenvalues of $A,B.$
That continuity argument, of course, doesn't extend to matrices over discrete fields, and this equality is true in any ring. If $R$ is a ring (with identity $I$) and $a,c,d\in R$ so that $dc=cd=I,(*)$ then:
$$ad(c-a)=(c-a)da,\tag{1}$$
because he left side is $adc-ada=a-ada$ and the right side is $cda-ada=a-ada,$ so they are equal.
Now given $a,b\in R$ so that $a+b$ is invertible in $R,$ let $c=a+b,d=(a+b)^{-1}$. Then $b=c-a$ so (1) becomes $$a(a+b)^{-1}b=b(a+b)^{-1}a$$
(*) There are rings where $cd=I$ does not imply $dc=I,$ but in square matrices, $DC=I $ means $CD=I.$ So for the general ring, we need the condition $cd=I$ and $dc=I.$
Best Answer
No, it's not an accident. The cross product is orthogonal to each factor, so the vector has to be orthogonal to $b\times c$, hence in the plane spanned by $b$ and $c$. But it also has to be orthogonal to $a$. So, writing $$a\times(b\times c) = xb + yc$$ and dotting with $a$, you get $x(b\cdot a) + y(c\cdot a)=0$. So the answer must be some scalar multiple of the correct formula. Now you only have to check that that scalar is $1$ by substituting $a=b$ and $a=c$. Better yet, let $a$ be a unit vector in the plane spanned by $b$ and $c$ that is orthogonal to $b$.