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, but I couldn’t solve it.

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], 
   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]

How can I use MMA 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.