This probability does tend to $\alpha:=\prod_{i=1}^\infty(1-2^{-i})$. In fact, this probability tends to $\alpha$ much more quickly than the probability that a single random $n\times n$ matrix is invertible. (I've written up a proof of this, which involves a bit more work than is needed to just show that the probability tends to $\alpha$.)
For an integer $n$, let $\alpha_n=\prod_{i=1}^n(1-2^{-i})$ be the probability that a random matrix in $\mathbb F_2^{n\times n}$ is nonsingular, and let $\alpha=\lim_{n\to\infty}\alpha_n\approx 0.2888$. Following Jyrki Lahtonen's comment, the probability that two random invertible matrices in $\mathbb F_2^{n\times n}$ have invertible sum is
$$p_n:=\frac1{\alpha_n}\Pr[\det A=\det(I+A)=1],$$
where $A\in\mathbb F_2^{n\times n}$ is chosen uniformly at random. We will find an explicit expression for $p_n$.
Given positive integers $n$ and $k$, let $F(n,k)$ denote the number of $k$-dimensional subspaces of $\mathbb F_2^n$. We black-box the fact
$$\sum_{k=0}^nF(n,k)(-1)^k2^{\binom k2}=\mathbf 1_{n=0}=\begin{cases}1&\text{if }n=0\\0&\text{if }n>0,\end{cases}\tag{$\star$}$$
which will be extremely useful. Our first use of this is the following.
Lemma. (Inclusion-exclusion in the subspace lattice) For a matrix $M$, we have
$$\mathbf1_{\det M=1}=\sum_{U\subset\mathbb F_2^n}(-1)^{\dim U}2^{\binom{\dim U}2}\mathbf1_{MU=0}.$$
Proof sketch. This follows quickly from applying ($\star$) within the kernel of $M$. $\square$
We can apply this lemma to find a formula for $\alpha_n$ as a sum, which will be useful later:
$$\alpha_n=\sum_{k=0}^nF(n,k)(-1)^k2^{\binom k2}2^{-nk}.$$
We now attack the problem. We first use our lemma to express $p_n$ as a sum:
$$p_n:=\frac1{\alpha_n}\sum_{U,V\subset\mathbb F_2^n}(-1)^{\dim U+\dim V}2^{\binom{\dim U}2+\binom{\dim V}2}\Pr[AU=(A+I)V=0].$$
The probability inside the sum is zero if $U$ and $V$ intersect nontrivially. If $U\cap V=\{0\}$, we claim that the events $[AU=0]$ and $[(A+I)V=0]$ are independent. It suffices to show this row-wise, i.e. that for a random vector $x\in\mathbb F_2^n$, $x\in U^\perp$ and $x+e_i\in V^\perp$ are independent, where $e_i$ is an elementary vector. This happens because the projections of $x$ onto $\mathbb F_2^n/U^\perp$ and $\mathbb F_2^n/V^\perp$ are independent because $U\cap V=\{0\}$; we omit the details.
We can now (more or less) explicitly compute $p_n$. From the above, we expand
\begin{align*}
p_n\alpha_n
&=\sum_{\substack{U,V\subset\mathbb F_2^n\\U\cap V=\emptyset}}(-1)^{\dim U+\dim V}2^{\binom{\dim U}2+\binom{\dim V}2}2^{-n(\dim U+\dim V)}\\
&=\sum_{U,V\subset\mathbb F_2^n}(-1)^{\dim U+\dim V}2^{\binom{\dim U}2+\binom{\dim V}2}2^{-n(\dim U+\dim V)}\left(\sum_{W\subset U,V}(-1)^{\dim W}2^{\binom{\dim W}2}\right)\\
&=\sum_{W\subset\mathbb F_2^n}(-1)^{\dim W}2^{\binom{\dim W}2}\sum_{W\subset U,V\subset\mathbb F_2^n}(-1)^{\dim U+\dim V}2^{\binom{\dim U}2+\binom{\dim V}2}2^{-n(\dim U+\dim V)}\\
&=\sum_{W\subset\mathbb F_2^n}(-1)^{\dim W}2^{\binom{\dim W}2}\left(\sum_{W\subset U\subset\mathbb F_2^n}(-1)^{\dim U}2^{\binom{\dim U}2}2^{-n\dim U}\right)^2\\
&=\sum_{W\subset\mathbb F_2^n}(-1)^{\dim W}2^{\binom{\dim W}2}\left(\sum_{\tilde U\in\mathbb F_2^n/W}(-1)^{\dim\tilde U+\dim W}2^{\binom{\dim\tilde U+\dim W}2}2^{-n(\dim\tilde U+\dim W)}\right)^2\\
&=\sum_{k=0}^nF(n,k)(-1)^k2^{\binom k2}\left(\sum_{\ell=0}^{n-k}F(n-k,\ell)(-1)^{\ell+k}2^{\binom{\ell+k}2}2^{-n(\ell+k)}\right)^2\\
&=\sum_{k=0}^nF(n,k)(-1)^k2^{\binom k2}\left((-1)^k2^{-nk}\sum_{\ell=0}^{n-k}F(n-k,\ell)(-1)^{\ell}2^{\binom\ell2+k\ell+\binom k2}2^{-n\ell}\right)^2\\
&=\sum_{k=0}^nF(n,k)(-1)^k2^{3\binom k2}2^{-2nk}\left(\sum_{\ell=0}^{n-k}F(n-k,\ell)(-1)^{\ell}2^{\binom\ell2}2^{-(n-k)\ell}\right)^2\\
&=\sum_{k=0}^nF(n,k)(-1)^k2^{3\binom k2}2^{-2nk}\alpha_{n-k}^2\\
&=\alpha_n^2-\frac{2^n-1}{2^{2n}}\alpha_{n-1}^2+\frac{(2^n-1)(2^{n-1}-1)}{3\cdot 2^{4n-3}}\alpha_{n-2}^2+\sum_{k=3}^nF(n,k)(-1)^k2^{3\binom k2}2^{-2nk}\alpha_{n-k}^2.
\end{align*}
The $k\geq 3$ terms in the above sum are negligible. We can compute
$$F(n,k)=\frac{(2^n-1)\cdots(2^n-2^{k-1})}{(2^k-1)\cdots(2^k-2^{k-1})}\leq \frac1\alpha2^{k(n-k)},$$
so the $k\geq 3$ terms are in total at most
$$\sum_{k=3}^\infty \frac1\alpha2^{-nk+\frac{k^2-3k}2}\leq 10\cdot 2^{-3n}$$
for large $n$. Doing the same to the $k=1$ and $k=2$ terms is enough to see that $p_n\to\alpha$.
To figure out how the error $|p_n-\alpha|$ behaves, we need to understand how $\alpha_n$ and $\alpha$ differ. We have
\begin{align*}
\log\alpha_n-\log\alpha=\sum_{m>n}-\log(1-2^{-m})
&=\sum_{m>n}\sum_{j=1}^\infty \frac{2^{-mj}}j\\
&=\sum_{j=1}^\infty\frac{2^{-nj}}{j(2^j-1)}=2^{-n}+\frac162^{-2n}+O(2^{-3n}),
\end{align*}
so
$$\alpha_n=\alpha\exp\left(2^{-n}+\frac162^{-2n}+O(2^{-3n})\right)=\alpha\left(1+2^{-n}+\frac232^{-2n}+O(2^{-3n})\right).$$
Now, we get
\begin{align*}
p_n
&=\alpha_n-\frac{2^n-1}{2^{2n}}\frac{\alpha_{n-1}^2}{\alpha_n}+\frac{(2^n-1)(2^{n-1}-1)}{3\cdot 2^{4n-3}}\frac{\alpha_{n-2}^2}{\alpha_n}+O(2^{-3n})\\
&=\alpha_n\left(1-\frac{1}{2^n-1}+\frac2{3(2^n-1)(2^{n-1}-1)}\right)+O(2^{-3n})\\
&=\alpha\left(1+2^{-n}+\frac232^{-2n}+O(2^{-3n})\right)\left(1-2^{-n}+\frac132^{-2n}+O(2^{-3n})\right)+O(2^{-3n})\\
&=\alpha\left(1+2^{-2n}+O(2^{-3n})\right).
\end{align*}
So, $p_n$ tends to $\alpha$ with error on the order of $4^{-n}$, while $\alpha_n$ tends to $\alpha$ with error on the order of $2^{-n}$.
