Heat exchange between bodies

Hello to everyone, I'm having a problem modeling the thermal exchange betwwen different objects. I'm trying to use the PDE toolbox, and I've understood how to set up geometries and how to create meshes. My question is: given a cold rectangular mesh and two hot circular meshes, is it possible to place the circles on the rectangle and see the heat exchange effect? For clarity I'm attaching the code I wrote. As you can see, I can overlap the geometries and see the final effect, however my aim is to sum the models, that should have a well specified (constant) temperature.

rect1 = [3
    4
    -1
    1
    1
    -1
    0
    0
    -0.5
    -0.5];
C1 = [1
    0.7
    -0.25
    0.15];

C2 = [1
    -.7
    -0.25
    0.15];


C1 = [C1;zeros(length(rect1) - length(C1),1)];
C2 = [C2;zeros(length(rect1) - length(C2),1)];
gd = [rect1,C1,C2];
ns = char('rect1','C1','C2');
ns = ns'
sf = 'rect1+ C1 + C2';
[dl,bt] = decsg(gd,sf,ns);
figure(1); 
pdegplot(dl,"EdgeLabels","on","FaceLabels","on")
xlim([-1.5,1.5])

model2 = createpde()
[dl2,bt2] = csgdel(dl,bt);
figure(2)
pdegplot(dl2,"EdgeLabels","on","FaceLabels","on")
geometryFromEdges(model2,dl2);
specifyCoefficients(model2,m=0,d=0,c=1,a=0,f=1);
applyBoundaryCondition(model2,'dirichlet', ...
                             'Edge',1:model2.Geometry.NumEdges, ...
                             'u',10);
mesh = generateMesh(model2,'Hmax',0.25);
results = solvepde(model2)
figure(3)
axis equal
pdeplot(model2,'XYData',results.NodalSolution)

Answer 1

[Atithi]

Define the geometry of the objects using rectangles and circles.

Combine the geometries by summing them to create a single geometry that includes all the objects.

Create a PDE model, specify coefficients, apply boundary conditions, generate a mesh, solve the PDE, and visualize the results.

I have modified the code a bit

rect1 = [3; 4; -1; 1; 1; -1; 0; 0; -0.5; -0.5];
C1 = [1; 0.7; -0.25; 0.15];
C2 = [1; -0.7; -0.25; 0.15];
C1 = [C1; zeros(length(rect1) - length(C1), 1)];
C2 = [C2; zeros(length(rect1) - length(C2), 1)];
gd = [rect1, C1, C2];
ns = char('rect1', 'C1', 'C2');
ns = ns';
sf = 'rect1 + C1 + C2';
[dl, bt] = decsg(gd, sf, ns);

model = createpde();

% Step 4: Remove overlapping regions and create a valid geometry
[dl2, bt2] = csgdel(dl, bt);

% Step 5: Generate a mesh
geometryFromEdges(model, dl2);
mesh = generateMesh(model, 'Hmax', 0.25);

% Step 6: Specify coefficients
specifyCoefficients(model, 'm', 0, 'd', 0, 'c', 1, 'a', 0, 'f', 1);

% Step 7: Apply boundary conditions
applyBoundaryCondition(model, 'dirichlet', 'Edge', 1:model.Geometry.NumEdges, 'u', 10);

% Step 8: Solve the PDE
results = solvepde(model);

% Step 9: Visualize the results
figure;
pdeplot(model, 'XYData', results.NodalSolution);
axis equal;

Do let me know if it is working for you?