In order to avoid nasty equations, and square roots in particular, I'll be doing some Lie circle geometry. Since there are various conventions for Lie geometry, I'll start by introducing mine. I'll represent a circle with center $(x,y)$ and radius $r$ as a vector
$$\begin{pmatrix}x\\y\\\pm r\\x^2+y^2-r^2\\1\end{pmatrix}$$
or a multiple thereof (homogeneous coordinates). A line $ax + by + c = 0$ will be represented as
$$\begin{pmatrix}a\\b\\\pm \sqrt{a^2+b^2}\\-2c\\0\end{pmatrix}$$
or a multiple thereof. Two general circles (i.e. circles or lines) touch if $\langle c_1,c_2\rangle:=c_1^T\cdot L\cdot c_2=0$ where
$$L = \begin{pmatrix}-2\\&-2\\&&2\\&&&0&1\\&&&1&0\end{pmatrix}$$
represents the Lie quadric. Every circle satisfies the quadratic form $\langle c,c\rangle=0$.
Now some coordinates. I'll pick $B=(0,0)$, $A$ on the $x$ axis and $C$ on the $y$ axis. I want to avoid square roots, and since these appear in the representation of a line, I have to make sure that there is a square under that root. The tangent half-angle substitution can help here. Without loss of generality, I can write the line $AC$ as
$$AC: (1-s^2)x + (2s)y + t = 0$$
for suitably chosen $s, t\in\mathbb R$. (The special case of $s=\infty$ would correspond to a vertical line, so $AC$ and $BC$ would coincide. We can safely ignore this case.) Likewise, the line through $D$ can be written as
$$BD: (1-u^2)x + (2u)y = 0$$
If you want, you can compute the coordinates of $A,C,D$ from these lines.
\begin{align*}
A &= \left(0, -\frac{t}{2s}\right) \\
C &= \left(\frac{t}{s^2+1}, 0\right) \\
D &= \frac{1}{2(s-u)(su+1)}\Bigl(2tu, t(u^2-1)\Bigr)
\end{align*}
Now one can turn the lines into Lie circle vectors as above:
\begin{align*}
AB &: (1, 0, 1, 0, 0)^T \\
BC &: (0, 1, 1, 0, 0)^T \\
AC &: (1-s^2, 2s, 1+s^2, -2t, 0)^T \\
BD &: (1-u^2, 2u, \pm(1+u^2), 0, 0)^T
\end{align*}
The sign of $BD$ has to be different for the two incircles touching it from different sides. I'll be using the positive sign for $\triangle ABD$ and the negative for the incircle of $\triangle BCD$.
Now to find the incircles, one can solve a special case of the problem of Apollonius: Find a circle which touches three general circles (which in this case are lines). That means multiplying the vectors of the corresponding three lines by $L$, then finding the kernel of the corresponding $3\times5$ matrix. That kernel has dimension $2$, and one solution will always be the point at infinity $P_\infty=(0,0,0,1,0)^T$. So if the other element of the basis is a vector $v$, then you're looking for a linear combination $c$ of $P_\infty$ and $v$ which satisfies $c^T\cdot L\cdot c=0$ but which is not equal to $P_\infty$ itself. This combination can be given explicitely as
$$ c = \langle v, v\rangle P_\infty - 2\langle P_\infty, v\rangle v
= (v^T\cdot L\cdot v)P_\infty - 2(P_\infty^T\cdot L\cdot v)v $$
Using this technique, one can find the relevant incircles as
$$
c_{ABC} = \begin{pmatrix}
2st(s-1) \\
2st(s-1) \\
2st(s-1) \\
t^2 \\
4s^2(s-1)^2
\end{pmatrix}
\qquad
c_{ABD} = \begin{pmatrix}
2st(s-u) \\
2stu(s-u) \\
2st(s-u) \\
t^2u^2 \\
4s^2(s-u)^2
\end{pmatrix}
\qquad
c_{BCD} = \begin{pmatrix}
2t(su+1)(s-1)(u+1) \\
2t(su+1)(s-1)(u-1) \\
2t(su+1)(s-1)(u-1) \\
t^2(u+1)^2 \\
4(su+1)^2(s-1)^2
\end{pmatrix}
$$
The radius in each case is the third coordinate divided by the fifth. To have equal radius for the latter two incircles, you solve
$$\bigl(2st(s-u)\bigr)\bigl(4(su+1)^2(s-1)^2\bigr) =
\bigl(2t(su+1)(s-1)(u-1)\bigr)\bigl(4s^2(s-u)^2\bigr)$$
Bringing everything to one side and factoring the result, you get
$$
t \cdot s \cdot (s - 1) \cdot (s-u) \cdot (su + 1) \cdot
(su^2 + s^2 - 2su + s - 1) = 0
$$
$t=0$ means $B$ lies on $AC$. $s=0$ means $AC$ is vertical, so again $B$ would lie on that line if $A$ does. $s=1$ means $AC$ is horizontal so $B$ would lie on that line if $C$ does. $u=s$ as well as $u=-1/s$ makes $AC$ parallel to $BD$. All of these are degenerate situations which we can discount here. So the only relevant factor is the last one.
Which means that given a triangle $ABC$, one can read $s,t$ from its line $AC$, then compute $u$ from that in order to find the appropriate $D$.
$$ u = \frac{s \pm \sqrt{s\left(1-s^2\right)}}s =
1\pm\sqrt{s^{-1}-s} $$
(If anyone can give this relation between $s$ and $u$ a more geometric interpretation, I'd be glad to hear about it.)
The equation from the question statement can be written as
\begin{align*}
\frac{1}{r'} - \frac{1}{r} &= \frac{1}{BD}
\\
\frac{4s^2(s-u)^2}{2st(s-u)} -
\frac{4s^2(s-1)^2}{2st(s-1)} &=
\frac{2(s-u)(su+1)}{\sqrt{\bigl(2tu\bigr)^2 + \bigl(t(u^2-1)\bigr)^2}}
\\
\frac{2s(s-u)}{t} -
\frac{2s(s-1)}{t} &=
\frac{2(s-u)(su+1)}{\lvert t(u^2+1)\rvert}
\\
\frac{t(s-u)(su+1)}{s(1-u)} &=
\lvert t(u^2+1)\rvert
\end{align*}
You can do a case distinction here.
\begin{align*}
\frac{t(s-u)(su+1)}{s(1-u)} &=
t(u^2+1)
&
\frac{t(s-u)(su+1)}{s(1-u)} &=
-t(u^2+1)
\\
tu(su^2 + s^2 - 2su + s - 1) &= 0
&
t(su^3 - s^2u + su - 2s + u) &= 0
\end{align*}
As you can see, the left equation is satisfied whenever the condition $(su^2 + s^2 - 2su + s - 1)=0$ (which we found above to characterize equal radius $r'$) holds. If $A$ has a positive $y$ coordinate and $C$ a positive $x$ coordinate, then we can choose $0<s<1$ to represent the line. Since the normal vectors of $AB$ and $BC$ point towards the inside of the triangle, it would be better to use $-\infty<s<-1$ instead. In that case, $1-s^2$ and $2s$ are both negative (i.e. point towards the inside as well), and $t$ must be positive. And as $D$ must lie in the first quadrant, we have $-1<u<0$ or $1<u<\infty$. But again, with the normal vector pointing into $\triangle ABD$ (and the flipped normal vector into $\triangle BCD$ we have $1<u<\infty$ as the sole possibility. With $t>0$ and $1<u<\infty$ we know that $t(u^2+1)>0$ so the left equation is indeed the relevant one. The $1<u$ also tells you that in the above formula, $u = 1 \pm \sqrt{s^{-1} - s}$ you have to choose the positive solution. The other would lead to excircles. This only holds for positive coordinates of $A$ and $C$, though.
This answer became considerably longer than anticipated. (At some point, I had a mistake somewhere which made the results look a lot easier than they actually were.) I guess there may be other approaches which will do this shorter. The key benefit in the above answer is its potential to evaluate things over polynomial rings, with little to no need to take square roots. As you saw, the formula to compute $u$ given $s$ was never actually used, since we found the underlying polynomial factor instead. Extensive computations over polynomials can be done by hand, but that's quite prone to errors. I did the above using a computer algebra system (Sage in my case), and would agree that not doing so makes the approach impractical.