finite element – Modeling shear strain in comsol

I am trying to model coated circular inclusion subjected to remote shear strain conditions in plain strain settings. To do this, I have created a square inclusion of length 30 [m] which models the infinite domain, within which lies the coated inclusion of diameter 1[m].

The remote shear strain is constant. To model this, I fixed the bottom of the square (fixed constraint), and applied shear stress on the three sides, where the shear stress is computed as: 2*G2*\varepsilon_xy where \varepsilon_xy is the remote shear strain. On the top surface, the stress acts along the x-direction, right edge: along +y-direction, and along left edge: -ve y-direction. I have set the Geometric nonlinearity option to geometrically linear so that the solver solves only for the small strain conditions. Additionally, I have provided roller type BC to the bottom-right point.

Is my setup correct? I am new to modeling elastic problems in comsol, and I find that the analytical solution and the FEM solution obtained with this setup do not match. I am not sure where the issue is. Any help is appreciated!