Lets examine the quadrilateral as a linearly interpolated surface. If the vertices are
$$\bbox{
\vec{p}_1 = \left [ \begin{matrix} x_1 \\ y_1 \\ z_1 \end{matrix} \right ]
} , \quad \bbox {
\vec{p}_2 = \left [ \begin{matrix} x_2 \\ y_2 \\ z_2 \end{matrix} \right ]
} , \quad \bbox {
\vec{p}_3 = \left [ \begin{matrix} x_3 \\ y_3 \\ z_3 \end{matrix} \right ]
} , \quad \bbox {
\vec{p}_4 = \left [ \begin{matrix} x_4 \\ y_4 \\ z_4 \end{matrix} \right ]
}$$
and we parametrise the surface using a vector-valued function $\vec{p}(u, v)$, with $0 \le u \le 1$, $0 \le v \le 1$, then
$$\bbox{ \vec{p}(u, v) = (1 - v)\biggr( (1-u) \vec{p}_1 + u \vec{p}_2 \biggr) + v \biggr( (1-u) \vec{p}_3 + u \vec{p}_4 \biggr) }$$
i.e.
$$\bbox{ \vec{p}(u, v) = \vec{p}_1 + u \left( \vec{p}_2 - \vec{p}_1 \right) + v \left( \vec{p}_3 - \vec{p}_1 \right) + u v \left( \vec{p}_4 - \vec{p}_3 + \vec{p}_1 - \vec{p}_2 \right) }$$
If the four vertices are nonplanar, the surface is curved. In all cases, the surface passes through
$$\bbox{ \vec{p}(0.5, 0.5) = \frac{ \vec{p}_1 + \vec{p}_2 + \vec{p}_3 + \vec{p}_4 }{4} }$$
The surface tangent vectors are
$$\bbox{ \begin{aligned}
\vec{p}_u (u , v) &= \frac{\partial \vec{p}(u, v)}{\partial u} = \vec{p}_2 - \vec{p}_1 + v \left( \vec{p}_4 - \vec{p}_3 - \vec{p}_2 + \vec{p}_1 \right)
\\
\vec{p}_v (u, v) &= \frac{\partial \vec{p}(u, v)}{\partial v} = \vec{p}_3 - \vec{p}_1 + u \left( \vec{p}_4 - \vec{p}_3 - \vec{p}_2 + \vec{p}_1 \right) \\
\end{aligned} }$$
and the exact area of this surface is
$$\bbox{ A = \int_0^1 d u \int_0^1 d v \; \vec{p}_u (u, v) \cdot \vec{p}_v (u, v) = \frac{ \left( \vec{p}_1 + \vec{p}_2 - \vec{p}_3 - \vec{p}_4 \right)\cdot\left( \vec{p}_1 - \vec{p}_2 + \vec{p}_3 - \vec{p}_4 \right) }{4} }$$
The surface normal vector is $\vec{n}(u, v)$,
$$\bbox{ \vec{n}(u , v) = \frac{\partial \vec{p}(u, v)}{\partial u} \times \frac{\partial \vec{p}(u, v)}{\partial v} }$$
and the mean (expected value) of the normal vector is
$$\bbox{ \langle\vec{n}\rangle = \int_0^1 d u \int_0^1 d v \; \vec{n}(u ,v) = \left [ \begin{matrix}
\frac{(y_4 - y_1)(z_3 - z_2) - (z_4 - z_1)(y_3 - y_2)}{2} \\
\frac{(z_4 - z_1)(x_3 - x_2) - (x_4 - x_1)(z_3 - z_2)}{2} \\
\frac{(x_4 - x_1)(y_3 - y_2) - (y_4 - y_1)(x_3 - x_2)}{2} \\
\end{matrix} \right] = \vec{n}\left(\frac{1}{2} , \frac{1}{2}\right) }$$
If we consider a flow through the surface, we could assume the flow is perpendicular to the surface normal. Then, the effective surface is the area of the surface when projected to a plane perpendicular to the unit normal vector $\hat{n}$ used,
for example
$$\hat{n} = \frac{\langle\vec{n}\rangle}{\left\lVert\langle\vec{n}\rangle\right\rVert}$$
$$A_\perp = \int_0^1 d u \int_0^1 d v \; \Biggr( \vec{p}_u (u, v) - \hat{n} \biggr(\hat{n}\cdot\vec{p}_u (u,v) \biggr) \Biggr) . \Biggr( \vec{p}_v (u, v) - \hat{n} \biggr(\hat{n}\cdot\vec{p}_v (u,v) \biggr) \Biggr)$$
It turns out that for this choice of $\hat{n}$, $A_\perp = A$.
Therefore, choosing
$$\bbox[#ffffef, 1em]{ \vec{n} = \left [ \begin{matrix}
(y_4 - y_1)(z_3 - z_2) - (z_4 - z_1)(y_3 - y_2) \\
(z_4 - z_1)(x_3 - x_2) - (x_4 - x_1)(z_3 - z_2) \\
(x_4 - x_1)(y_3 - y_2) - (y_4 - y_1)(x_3 - x_2) \end{matrix} \right ] = \; \bigr( \vec{p}_4 - \vec{p}_1 \bigr) \times \bigr( \vec{p}_3 - \vec{p}_2 \bigr) }$$
$$\bbox{\hat{n} = \frac{\vec{n}}{\left\lVert\vec{n}\right\rVert} = \frac{\vec{n}}{\sqrt{\vec{n} \cdot \vec{n}}}}$$
$$\bbox[#ffffef, 1em]{ A = \frac{ \left( \vec{p}_1 + \vec{p}_2 - \vec{p}_3 - \vec{p}_4 \right)\cdot\left( \vec{p}_1 - \vec{p}_2 + \vec{p}_3 - \vec{p}_4 \right) }{4} }$$
is an obvious choice, in my opinion.
According to Bretschneider’s formula
$$K={\sqrt {(s-a)(s-b)(s-c)(s-d)-abcd\cdot \cos ^{2}\left({\frac {\alpha +\gamma }{2}}\right)}}$$
when $a,b,c,d$ are fixed, the area depends only on $\cos(\alpha + \gamma)$, where $\alpha$ and $\gamma$ are opposite included angles and $s$ is the semiperimeter. So, we have to find values of $\alpha$ and $\gamma$ that make $\alpha + \gamma$ as close as possible to $\pi/2$ and $\pi$.
The quadrilateral can be made cyclic, which gives us $\alpha+\gamma = \pi$. This minimizes $\cos(\alpha + \gamma)$ and hence maximizes the area.
The maximal possible area is therefore
$${\sqrt {(s-a)(s-b)(s-c)(s-d)}}.$$
Best Answer
Let $R$ be the radius of the inscribed circle.
Let $C$ be the tangent point of the inscribed circle with $AO_2$.
Let $O'$ be the center of the inscribed circle.
Since $\angle{O'CA}=90^\circ$, we have $$\text{Area}(O'AO_2)=\frac 12\times O'C\times AO_2=\frac{rR}{2}$$
Similarly, we have $$\text{Area}(O'AO_1)=\text{Area}(O'O_1B)=\text{Area}(O'BO_2)=\frac{rR}{2}$$
Therefore, the area of the quadrilateral $AO_1BO_2$ is $$4\times\frac{rR}{2}=\color{red}{2rR}$$
Added 1 :
I realized that I did not use the condition that $O_1 \in C_2, O_2 \in C_1$ from which we can say that $O_1O_2=r$.
So, we can say that $\triangle{AO_1O_2}$ is an equilateral triangle, so we have $$\text{Area}(AO_1O_2)=\frac 12\times r\times \frac{\sqrt 3}{2}r=\frac{\sqrt 3}{4}r^2$$
Therefore, the area of $AO_1BO_2$ is $$2\times\frac{\sqrt 3}{4}r^2=\color{red}{\frac{\sqrt 3}{2}r^2}$$
Added 2 :
The condition $O_1 \in C_2$ means that $O_1$ is on the circle $C_2$.
The condition $O_2 \in C_1$ means that $O_2$ is on the circle $C_1$.
It follows from these that $O_1O_2=r$.
Your drawing is not correct.
$O_1$ has to be on the circle $C_2$, and $O_2$ has to be on the circle $C_1$.