I want to solve the plane stress problem of the following stress boundary:

I built a system of equations based on the stress balance equation and the deformation compatibility equation.

Needs["NDSolve`FEM`"]
Ω = 
  RegionDifference[Rectangle[{-1, -1}, {1, 1}], Disk[{0, 0}, 0.1]];
Region[Ω]


op = { D[σx[x, y], x] + D[τxy[x, y], y], 
  D[σy[x, y], y] + D[τxy[x, y], x], 
  2 D[τxy[x, y], x, y] + D[σx[x, y], x, x] + 
   D[σy[x, y], y, y]}
(*∂Subscript[σ, \
xx](x,y)/∂x+∂Subscript[τ, xy](x,y)/\
∂y\[Equal]0
∂Subscript[σ, \
yy](x,y)/∂y+∂Subscript[τ, xy](x,y)/\
∂x\[Equal]0;*)

Γ = {DirichletCondition[{σx[x, y] == 
      0., σy[x, y] == 0., τxy[x, y] == 0.}, 
    x^2 + y^2 == 0.1^2], 
   DirichletCondition[{σx[x, y] == 10., σy[x, y] == 
      0., τxy[x, y] == 0.}, x == 1 && -1 <= y <= 1], 
   DirichletCondition[{σx[x, y] == -10., σy[x, y] == 
      0., τxy[x, y] == 0.}, x == -1 && -1 <= y <= 1], 
   DirichletCondition[{σx[x, y] == 0., σy[x, y] == 
      10., τxy[x, y] == 0.}, y == 1 && -1 <= x <= 1], 
   DirichletCondition[{σx[x, y] == 
      0., σy[x, y] == -10., τxy[x, y] == 0.}, 
    y == -1 && -1 <= x <= 1]};


{ufun, vfun, wfun} = 
 NDSolveValue[{op == {0, 0, 
     0}, Γ}, {σx, σy, τxy}, {x, 
    y} ∈ Ω,  StartingStepSize -> 0.1, 
  MaxStepSize -> 0.01, WorkingPrecision -> 20]

ContourPlot[ufun[x, y], {x, y} ∈ Ω, 
 ColorFunction -> "Temperature", AspectRatio -> Automatic, 
 PlotPoints -> 30, WorkingPrecision -> 20, Contours -> Range[0, 5, 1],
  PlotRange -> Full]

But the solution result is obviously wrong:

How can I use Mathematica to solve this kind of plane stress problem?

There is a similar post here, but I would like to know if there is a general method to solve this type of problem that does not require additional processing skills.