FEM2D: Thermal convection problems using Triangular Elements

Consider the 2D equation model, a second order PDE equation with general coefficients a11, a12, a21, a22, a00, f :

Contents

Equations for Thermal problems

Assume the case and .

State the problem:

Consider an square domain [-1,1]x[-1,1] with a circular hole in the middle (CircleHoleMesh01.m)

As a BC, consider the circle points at 15ºC and bottom boundary at 50ºC. Moreover, assume top boundary has convection with beta=2, and Tinf=-5ºC.

Load Geometry: node coordinates and elements

clear all;
close all;
eval('CircleHolemesh01');
numNod=size(nodes,1)
numElem=size(elem,1)
numbering=0;
plotElements(nodes,elem,numbering);

numNod =

   260


numElem =

   427

Select Boundary points

indTop=find(nodes(:,2)> 0.99);
indBot=find(nodes(:,2)< -0.99);
indCircle=find(sqrt(nodes(:,1).^2+nodes(:,2).^2)<0.501);
hold on;
plot(nodes(indTop,1),nodes(indTop,2),'mo');
plot(nodes(indBot,1),nodes(indBot,2),'mo');
plot(nodes(indCircle,1),nodes(indCircle,2),'yd');
hold off;

Define Coefficients vector of the model equation

In this case we use the Poisson coefficients defined in the problem above

a11=1;
a12=0;
a21=a12;
a22=a11;
a00=0;
f=0;
coeff=[a11,a12,a21,a22,a00,f];

Compute the global stiff matrix

Use the linearTriangElement function to compute the system by element

K=zeros(numNod);   %global Stiff Matrix
F=zeros(numNod,1); %global internal forces vector
Q=zeros(numNod,1); %global secondary variables vector
for e=1:numElem
  % compute element system
    [Ke,Fe]=linearTriangElement(coeff,nodes,elem,e);
  %
  % Assemble the elements
  %
    rows=[elem(e,1); elem(e,2); elem(e,3)];
    colums= rows;
    K(rows,colums)=K(rows,colums)+Ke; %assembly
    if (coeff(6) ~= 0)
        F(rows)=F(rows)+Fe;
    end
end % end for elements
% we save a copy of the initial F array
% for the postprocess step
Kini=K;
Fini=F;

Boundary Conditions (BC)

Impose Boundary Conditions for this example. In this case only essential and natural BC are considered. We will do a thermal example for convection BC (mixed).

fixedNodes=[indBot', indCircle']; %Fixed Nodes (global num.)
% Obs: No cal fer servir unique perque les fronteres son disjuntes
freeNodes = setdiff(1:numNod,fixedNodes); %Complementary of fixed nodes

% ------------ Convection BC
convecNodes=indTop';
beta=2; %Convection Coefficient
Tinf=-5; %Bulk temperature
[K,Q]=applyConvTriang(convecNodes,beta,Tinf,K,Q,nodes,elem);

% ------------ Essential BC
u=zeros(numNod,1); %initialize u vector
u(indBot)=50; %all of them are zero
u(indCircle)=15; %all of them are zero
F=F-K(:,fixedNodes)*u(fixedNodes);%here u can be different from zero only for fixed nodes

% modify the linear system.
Km=K(freeNodes,freeNodes);
Im=F(freeNodes)+Q(freeNodes);

Compute the solution

solve the System

format short e; %just to a better view of the numbers
um=Km\Im;
u(freeNodes)=um;

PostProcess: Compute secondary variables and plot results

titol='Temperature Distribution';
colorScale='jet';
plotContourSolution(nodes,elem,u,titol,colorScale);

Exercise 1:

Compute the temperature for the point p=[0.5, 0.2].

Sol: pElem = 219, nodes=[217, 199, 210], tempP = 1.4846e+01

(c)Numerical Factory 2018

Published with MATLAB® R2018a