implementation of shell elements in a topology optimization algorithm

I am working on developing a topology optimization solver based on the finite element method and I want to add a triangular shell element in it. I used the classic finite element method but I didn’t separate in-plane and off-plane strains. My strain-displacement matrix consists of the following:

$$
B_i=left(
begin{array}{cccccc}
frac{partial N_i}{partial x} & 0 & 0 & 0 & zfrac{partial N_i}{partial x} & 0
0 & frac{partial N_i}{partial y} & 0 & -zfrac{partial N_i}{partial y} & 0 & 0
0 & 0 & 0 & 0 & 0 & 0
frac{partial N_i}{partial y} & frac{partial N_i}{partial x} & 0 & -zfrac{partial N_i}{partial x} & zfrac{partial N_i}{partial y} & 0
0 & 0 & frac{partial N_i}{partial y} & 0 & kappa N_i & 0
0 & 0 & frac{partial N_i}{partial x} & -kappa N_i & 0 & 0
end{array}right)
$$

And the shape function is:

$$
N(xi,eta)=begin{cases}
N_1=1-xi-eta
N_2=xi
N_3=eta
end{cases}
$$

I tried to validate my implementation with an example a cantilever beam: fixed on the left side and a concentrated force applied in the lower right corner and here is my Von Mises stresses distribution:

And my topology optimization is (optimal criteria used):

The results I obtain from Abaqus with the same properties, for Von Mises stresses distribution are depicted as follows:

The results that I get from my code are different from the results of Abaqus and I am not sure why. But here are their solutions and modifications I tried without any tangible difference in the obtained results:

• I used a transformation matrix otherwise the calculation doesn’t run
• The transformation matrix applied only in the calculation of the stiffness matrix (its dimension is $18times 18$)
• I tried to formulate the shape function in $(x,y)$ basis instead of $(xi,eta)$ but there was no difference

What do you think might be the source of error in my implementation? If you need more information about the way I implemented my code please tell me to clarify.