Poisson PDE in a rectangular domain

Question

I found a solution for this problem, but this is in Scilab and I never use Scilab. Can anyone can help me to translate it in Mathematica? Here is the Link:http://imechanica.org/files/TorqueR.pdf

I solved this problem with NDSolve but now I want to implement it with FDM, like this.

Clear[a, b, n, m, Δx, Δy, G, θ, 
  sol, vars, ϕsol, eqns, tauyz, tauxz, ϕ, gammaxz, gammayz]
a = 0.04;
b = 0.04;
n = 5;
m = 5;
Δx = a/n;
Δy = b/m;
G = 13.125 10^6;
θ = 0.002;
ϕ[0, j_] = 0;
ϕ[n, j_] = 0;
ϕ[i_, 0] = 0;
ϕ[i_, m] = 0;
Ψ[i_, j_] = -2 G θ;
vars = Flatten[Table[ϕ[i, j], {i, 1, n - 1}, {j, 1, m - 1}]];
eqns = 
  Flatten[
    Table[
      (ϕ[i + 1, j] - 2 ϕ[i, j] + ϕ[i - 1, j])/Δx^2 + 
      (ϕ[i, j + 1] - 2 ϕ[i, j] + ϕ[i, j - 1])/Δy^2 == Ψ[i, j], 
      {i, 1, n - 1}, {j, 1, m - 1}]];
sol = Solve[eqns, vars][[1]];
ϕsol = 
  Interpolation[
    Flatten[Table[{i Δx, j Δy, ϕ[i, j]}, {i, 0, n}, {j, 0, m}] /. sol, 1]];
tauyz = -D[ϕsol[x, y], x];
tauxz = D[ϕsol[x, y], y];
gammayz = tauyz/G;
gammaxz = tauxz/G;

but this only works for a rectangle, but I need to solve for a rectangle with hole in the center, like this:

Answer 1

This is a very simple way to do it. I modified your code very little to add the condition, that is inside a hole, then set the solution to be the boundary condition at the edge of the hole inside, which I set to be zero.

Added one line:

  startRow = 4; endRow = 6; startCol = 4; endCol = 6;

Which tells the boundaries of the hole. And inside the Table added an If to check:

If[i >= startRow && i <= endRow && j >= startCol && j <= endCol,
 ϕ[i, j] == 0,
 .... same as before
 ],

With more mesh, the solution should become better.

Here is your code with the small modification

Clear[a, b, n, m, Δx, Δy, G, θ, sol, vars, ϕsol, eqns, tauyz, tauxz, ϕ, gammaxz, \
gammayz]
a = 0.04;
b = 0.04;
n = 100;
m = 100;
Δx = a/n;
Δy = b/m;
G = 13.125 10^6;
θ = 0.002;
ϕ[0, j_] = 0;
ϕ[n, j_] = 0;
ϕ[i_, 0] = 0;
ϕ[i_, m] = 0;
Ψ[i_, j_] = -2 G θ;
startRow = 40; endRow = 60; startCol = 40; endCol = 60;
vars = Flatten[Table[ϕ[i, j], {i, 1, n - 1}, {j, 1, m - 1}]];

eqns = Table[
   If[i >= startRow && i <= endRow && j >= startCol && j <= endCol,
    ϕ[i, j] == 0, (*assume B.C. in hole edge is zero *)
    (ϕ[i + 1, j] - 2 ϕ[i, j] + ϕ[i - 1, j])/Δx^2 + (ϕ[i, j + 1] - 
         2 ϕ[i, j] + ϕ[i, j - 1])/Δy^2 == Ψ[i, j]
    ],
   {i, 1, n - 1}, {j, 1, m - 1}
   ];

sol = Solve[Flatten@eqns, vars][[1]];

ϕsol = Interpolation[Flatten[Table[{i Δx, j Δy, ϕ[i, j]}, {i, 0, n}, {j, 0, m}] /. sol, 1]];
tauyz = -D[ϕsol[x, y], x];
tauxz = D[ϕsol[x, y], y];
gammayz = tauyz/G;
gammaxz = tauxz/G;