After working this out and writing it down, I found that it was already done in 1969: Eric Langford, The probability that a random triangle is obtuse, Biometrika (1969) 56 (3): 689-690.
For this kind of problem, it's important to make good use of the symmetries and choose convenient coordinates. A significant simplification is achieved by noting that a triangle has at most one obtuse angle; thus the three events of one of the angles being obtuse are mutually exclusive, and thus the probability of the triangle being obtuse is just three times the probability $p$ of any given one of the angles being obtuse. Thus you don't have to deal with both the circle and the strip between the perpendiculars, since the probability of the third point being in the circle is $p$ and the probability of the third point being in the strip is $1-2p$ and you can obtain the desired probability $1-3p$ from either of them. In fact, you don't even need the strip, since the probability of the third point being beyond one of the perpendiculars is also $p$.
So there are three different approaches; we can either integrate either the area beyond one of the perpendiculars or the area of the strip or the area of the circle over all positions of two points. The advantage of using the perpendiculars or the strip is that you have to deal only with straight lines and not with circles; the disadvantage is that the integration region is complicated due to the corners. The advantage of the second approach is that you avoid the corners and get a relatively simple integration region; the disadvantage is that you get somewhat complicated integrands for the circle segments that extend beyond the square.
I'll work out the circle approach here. The easy part is integrating over the area of the circle without regard to whether it's inside the square. This is
$ \def\lint{\int\limits} \begin{eqnarray} && \lint_0^1\mathrm dx_1\lint_0^1\mathrm dx_2\lint_0^1\mathrm dy_1\lint_0^1\mathrm dy_2\,\pi\left(\left(\frac{x_2-x_1}2\right)^2+\left(\frac{y_2-y_1}2\right)^2\right) \\ &=& 2\cdot\frac\pi4\lint_0^1\mathrm dx_1\lint_0^1\mathrm dx_2\,(x_2-x_1)^2 \\ &=& \frac\pi2\left(\frac13-2\cdot\frac12\cdot\frac12+\frac13\right) \\ &=& \frac\pi{12}\;. \end{eqnarray} $
From this we have to subtract the integral over the circle segments that extend beyond the square. By symmetry, this is four times the integral over the segments that extend beyond the bottom edge of the square. These segments depend on $x_1$ and $x_2$ only through $x:=|x_2-x_1|$, so we can get rid of one integral by the transformation
$ \lint_0^1\mathrm dx_1\lint_0^1\mathrm dx_2\to2\lint_0^1\mathrm dx_1\lint_{x_1}^1\mathrm dx_2\to2\lint_0^1\mathrm dx\lint_x^1\mathrm dx_2\to2\lint_0^1\mathrm dx\,(1-x)\;. $
It will also be convenient to transform from $y_1$ and $y_2$ to $s:=y_1+y_2$ and $a:=y_2-y_1$. We can restrict to $a\gt0$ by symmetry, and the symmetry factor of $2$ cancels the Jacobi determinant $2$ of the transformation:
$\lint_0^1\mathrm dy_1\lint_0^1\mathrm dy_2\to\lint_0^1\mathrm da\lint_a^{2-a}\mathrm ds\;.$
The centre of the circle is at $y=s/2$, and its squared radius is $(x/2)^2+(a/2)^2$, so the angle of the segment extending beyond the lower edge is
$\alpha=2\arccos\frac s{\sqrt{x^2+a^2}}\;.$
The area of a circular segment with angle $\alpha$ and radius $r$ is $(\alpha-\sin\alpha)r^2/2$. Thus, with $\sin2\arccos u=2u\sqrt{1-u^2}$, the area is
$A(x,a,s)=\frac14\left(\left(x^2+a^2\right)\arccos\frac s{\sqrt{x^2+a^2}}-s\sqrt{x^2+a^2-s^2}\right)\;.$
The factor $1/4$ cancels with the symmetry factor $4$ for the four edges of the square. We need to integrate this over all values for which there is a circle segment extending beyond the lower edge, that is, for $s^2\lt x^2+a^2$. We also have $s\lt2-a$. These two bounds coincide for $x^2+a^2=(2-a)^2$, that is $a=1-x^2/4$, so we need two integrals for $a$ above and below that value:
$\lint_0^1\mathrm dx\,(1-x)\lint_0^{1-x^2/4}\mathrm da\lint_a^{\sqrt{x^2+a^2}}\mathrm ds\,A(x,a,s)+\lint_0^1\mathrm dx\,(1-x)\lint_{1-x^2/4}^1\mathrm da\lint_a^{2-a}\mathrm ds\,A(x,a,s)\;.$
The integrals over the second term in $A(x,a,s)$ are relatively easy to evaluate:
$\lint_0^1\mathrm dx\,(1-x)\lint_0^{1-x^2/4}\mathrm da\lint_a^{\sqrt{x^2+a^2}}\mathrm ds\,s\sqrt{x^2+a^2-s^2}=\frac{37}{2520}$
(computation) and
$\lint_0^1\mathrm dx\,(1-x)\lint_{1-x^2/4}^1\mathrm da\lint_a^{2-a}\mathrm ds\,s\sqrt{x^2+a^2-s^2}=\frac1{840}$
(computation), and the sum is $\displaystyle\frac{37}{2520}+\frac1{840}=\frac1{63}$.
The first term is a bit harder to integrate; Wolfram|Alpha fails to do it in one go. If $F(s)$ is an antiderivative for this term with respect to $s$ (where I'm suppressing $a$ and $x$ to avoid notational clutter), then the integral we need is
$\lint_0^1\mathrm dx\,(1-x)\left(\lint_0^{1-x^2/4}\mathrm da\,F\left(\sqrt{x^2+a^2}\right)+\lint_{1-x^2/4}^1\mathrm da\,F(2-a)-\lint_0^1\mathrm da\,F(a)\right)\;.$
With
$F(s)=\left(x^2+a^2\right)\left(s\arccos\frac s{\sqrt{x^2+a^2}}-\sqrt{x^2+a^2-s^2}\right)$
(computation) we have
$F\left(\sqrt{x^2+a^2}\right)=0\;,$
$F(2-a)=\left(x^2+a^2\right)\left((2-a)\arccos\frac{2-a}{\sqrt{x^2+a^2}}-\sqrt{x^2+4a-4}\right)\;,$
$F(a)=\left(x^2+a^2\right)\left(a\arccos\frac a{\sqrt{x^2+a^2}}-x\right)\;.$
Wolfram|Alpha evaluates the integral over $F(a)$ in one go if you give it extra time:
$\lint_0^1\mathrm dx\,(1-x)\lint_0^1\mathrm da\,F(a)=\frac{-52+21\pi-48\log2}{720}$
(computation). The hardest part is the integral over $F(2-a)$. The antiderivative with respect to $a$ looks rather daunting at first, until you realize that it's zero at the lower limit $a=1-x^2/4$. Thus we only need it for $a=1$. The result is still complicated but can be simplified to
$\frac1{12}\Bigg(\left(18x^2+5\right)\arccos\frac1{\sqrt{x^2+1}}-\frac1{15}x\left(6x^4+245x^2+75\right) \\ +x^3\left(8\log\left(x^2+1\right)-3x\arctan x\right)\Bigg).$
A final integration then yields
$\lint_0^1\mathrm dx\,(1-x)\lint_{1-x^2/4}^1\mathrm da\,F(2-a)=-\frac{517}{3150}+\frac{\pi-\log 2}{15}$
(computation).
Putting it all together, we obtain
$ \begin{eqnarray} p &=& \frac\pi{12}-2\left(-\frac{517}{3150}+\frac{\pi-\log 2}{15}-\frac{-52+21\pi-48\log2}{720}-\frac1{63}\right) \\ &=& \frac{97}{450}+\frac\pi{120}\;. \end{eqnarray} $
(computation). Thus, the probability $1-3p$ of the triangle being acute is
$\boxed{\displaystyle\frac{53}{150}-\frac\pi{40}\approx0.2748}$
in agreement with numerical simulations.
While we're at it, we can also calculate the probability for points uniformly distributed on the boundary of the square. We can assume the first point to be on one side, say the left one, and average the probability of the angle at the third point being obtuse as the second and third points range over all sides.
If the second point is on the same side as the first, the angle is obtuse only if the third point is also on the same side and between the other two; the probability for that is $1/3$ by symmetry.
If the second point is on the opposite side from the first, the angle is obtuse if the third point is also on one of these sides and between the other two, for a contribution of $2\cdot1/3$. It can also be obtuse if the third point is on one of the remaining sides. The width $w$ of the interval on the bottom edge for which this is the case is determined by the circle equation:
$ \left(\frac12\right)^2+\left(\frac{y_1-y_2}2\right)^2=\left(\frac w2\right)^2+\left(\frac{y_1+y_2}2\right)^2\;,\\ w^2=1-4y_1y_2 \;. $
The corresponding integrals are
$\lint_0^{1/4}\mathrm dy_1\lint_0^1\mathrm dy_2\,\sqrt{1-4y_1y_2}=\frac49-\frac{\log2}3$
(computation) and
$\lint_{1/4}^1\mathrm dy_1\lint_0^{1/(4y_1)}\mathrm dy_2\,\sqrt{1-4y_1y_2}=\frac{\log2}3$
(computation) for a total contribution of $8/9$ (multiplied by $2$ for the two remaining edges).
If the second point is on a side adjacent to the left one, say, the bottom one, the angle is obtuse if the third point is on one of these two sides and between the other points, which gives two contributions of $1/2$. The angle can also be obtuse if the third point is on one of the remaining sides, say, the right one. Again, the width of the interval for which this is the case is determined by the circle equation:
$ \left(\frac w2\right)^2+\left(1-\frac x2\right)^2=\left(\frac x2\right)^2+\left(\frac y2\right)^2\;,\\ w^2=y^2+4x-4 \;. $
The corresponding integral is
$\lint_0^1\mathrm dy\lint_{1-y^2/4}^1\mathrm dx\,\sqrt{y^2+4x-4}=\frac1{24}$
(computation).
Putting it all together, we have
$ \begin{eqnarray} p &=& \frac14\cdot\frac14\left(\frac13+\frac23+\frac89+2\left(2\cdot\frac12+2\cdot\frac1{24}\right)\right) \\ &=& \frac{73}{288}\;, \end{eqnarray} $
and the probability $1-3p$ of the triangle being acute is
$ \frac{23}{96}\approx0.2396 $
in agreement with numerical simulations.
Here's a table with the probabilities of three points forming an acute triangle for various distributions of the points in the plane:
$ \begin{array}{c|c|c} \text{circle}&\frac14&0.25\\ \hline \text{disk}&\frac4{\pi^2}-\frac18&0.2803\\ \hline \partial(\text{square})&\frac{23}{96}&0.2396\\ \hline \text{square}&\frac{53}{150}-\frac\pi{40}&0.2748\\ \hline \mathbb R^6\text{ symmetry}&\frac14&0.25 \end{array} $
References for these values are at MathWorld under Disk Triangle Picking, Square Triangle Picking and Gaussian Triangle Picking and in an article by Stephen Portnoy, A Lewis Carroll Pillow Problem: Probability of an Obtuse Triangle, Statistical Science, Volume 9, Number 2 (1994), 279-284 that contains a good discussion of what it might mean to choose triangles randomly and ends with the brilliant sentence "I would suggest that the value and enjoyment of mathematics is greatest when it is used to provide imperfect and incomplete solutions to ill-formulated problems".
The entry "$\mathbb R^6\text{ symmetry}$" refers to a part of Portnoy's article where he shows that any distribution invariant under rotations in $\mathbb R^6$ (treating all point coordinates as a single vector in $\mathbb R^6$) must lead to the same result and derives that result for a Gaussian distribution. Here's a simple argument for the result $1/4$ that I haven't found anywhere: Parametrize the coordinates $x_1,y_1,x_2,y_2$ of the first two points with $\phi,a_1,a_2,h$, where $\phi$ is the angle through which they have to be rotated to let the vector from one to the other point along the positive $x$ direction, $a_1$ and $a_2$ are the $x$ coordinates of the rotated points and $h$ is their common $y$ coordinate. The Jacobian of this transformation is $|a_1-a_2|$. Thus, to find the probability that one of the angles at the first two points is obtuse, we can average over all values of $\phi$, $a_1$, $a_2$, $h$ and the rotated coordinates $a_3,h_3$ of the third point, weighting each configuration by the distance between the first two points. Since $\phi$, $h$ and $h_3$ are independent of $a_1$, $a_2$ and $a_3$ and the obtuseness of the angles doesn't depend on them, we only have to average over $a_1$, $a_2$ and $a_3$. By symmetry, all six permutations of $a_1$, $a_2$ and $a_3$ are equally likely. They form three pairs of mirror images. The weights of the two pairs that have one of the first two points in the middle add up to to the weight of the pair that has the third point in the middle. Thus the probability for one of two angles to be obtuse is $1/2$, and hence the probability for one of three angles to be obtuse is $3/4$.