Finite Difference - heat generation in a square (square within a square)

Question

I am trying to develop a code that will allow me to solve the following problem. Lets say that I have a square block, inside that square block I have a smaller square box thats at a certain temperature. The outside of the outer square is also at a known temperature. What I am trying to find is the temperature distribution in the larger square due to heat transfer from the small square.

The following code allows me to find the temperature distribution in the larger square.

    diffFormula = (T[i + 1, k] - 2 T[i, k] + T[i - 1, k])/x^2 +(T[i, k + 1] - 2 T[i, k] + T[i, k - 1])/y^2 == 0

The boundary conditions are

NP=5; (*number of grid points*) (*so a 6x6 matrix*)
x=y=0.05;
Table[T[0, k] = 325, {k, 0, NP}];
Table[T[i, NP] = 325, {i, 0, NP}];
Table[T[i, 0] = 325, {i, 0, NP}];
Table[T[NP, k] = 325, {k, 0, NP}];

Now I setup the list of equations

eq = Flatten[Table[diffFormula, {i, 1, NP - 1}, {k, 1, NP - 1}]];

Followed by a list of variables

var = Flatten[Table[T[i, k], {i, 1, NP-1}, {k, 1, NP-1}]];

Followed by the solution

solution = N[Flatten[Solve[eq, var]]]
s = N[Table[T[i, k], {k, 0, NP}, {i, 0, NP}]] /. solution;

My problem is, lets say that I would like positions T[2,2],T[2,3],T[3,2],T[3,3] to represent the smaller square that is generating heat at 400. Simply setting these coordinates to 400 does not work as it passes this value to the variables list which causes interference with Solve. I've tried many ways of making the adjustment, however I can not find the right way to go about this.

I would really appreciate any help, Thank you.

Answer 1

Just remove the undesirable equations and variables:

Clear@"`*";
diffFormula = (T[i + 1, k] - 2 T[i, k] + T[i - 1, k])/x^2 + 
              (T[i, k + 1] - 2 T[i, k] + T[i, k - 1])/y^2 == 0;

NP = 5;(*number of grid points*)
x = y = 0.05;
Table[T[0, k] = 325., {k, 0, NP}];
Table[T[i, NP] = 325., {i, 0, NP}];
Table[T[i, 0] = 325., {i, 0, NP}];
Table[T[NP, k] = 325., {k, 0, NP}];

inner = {T[2, 2], T[2, 3], T[3, 2], T[3, 3]};

eq = Flatten[Table[diffFormula, {i, 1, NP - 1}, {k, 1, NP - 1}]];
selectedeq = Select[eq, Count[#, Alternatives @@ inner, Infinity] < 2 &];

var = Flatten[Table[T[i, k], {i, 1, NP - 1}, {k, 1, NP - 1}]];
selectedvar = Complement[var, inner];

(# = 400) & /@ inner;

solution = Flatten@Solve[selectedeq, selectedvar];
s = Table[T[i, k], {k, 0, NP}, {i, 0, NP}] /. solution;

ListPlot@s

BTW, though the region you're dealing with is luckily simple, handling irregular region with FDM can be really frustrating. This post is an example. If FDM isn't necessary for you, have a look at the FEM capabilities of NDSolve new added in v10.