CODEWITHSUNDEEP
matlabfinite-element-analysis
Sorade
Sorade

Reputation: 937

Setting Boundary Conditions in Finite Element

EDIT: See solution in the full code, following @bgb2's comment.

I'm currently trying to code a Finite Element Analysis to solve a 2D heat conduction problem. For now I'm looking at the steady state system:

enter image description here

where k is the thermal conductivity of the material.

Here is my code. I have verified that the stiffness matrix is correct, so I really think I'm doing something wrong when assigning the boundary conditions before solving for temperature.

Here is my commented code:

%% Input Geometry
% Local node geometry
% 4____3
% |    |
% |    |
% 1____2
% Global geometry
% 2____3____6(T=10)
% |    |    |
% |    |    |
% 1____4____5
% (T=3.5)
% Node 1 is at 3.5 degrees, and node 6 at 10 degrees.

% node global coordinates
node = [ 0   0  ;
         0   0.05;
         0.05 0.05;
         0.05 0  ;
         .1   0  ;
         .1   0.05];
% element connections, row=element, col=nodes
conn = [1, 4, 3, 2;
        4, 5, 6, 3];

%% Material Inputs
k_r = 1.4;

%% Boundary Conditions
T = zeros(6,1);
T(1) = 3.5; % node 1 fixed at 3.5 degrees
T(6) = 10; % node 6 fixed at 10 degrees

isol = [ 2, 3, 4, 5]; % unconstrained dofs,

%% Integration inputs
% Gauss Quadrature Points and weights for 2D quadrilateral elements
% local coordinates
GQ_coord = [-sqrt(1/3) -sqrt(1/3);
             sqrt(1/3) -sqrt(1/3);
            -sqrt(1/3)  sqrt(1/3);
             sqrt(1/3)  sqrt(1/3)];
% GQ weights
GQ_w = [1 1 1 1]; 

%% SOLVER
nn = size(node,1); % number of nodes
ndof = nn; % number of dofs in the problem
ne = size(conn,1); % number of elements

K = zeros(ndof, ndof); %global stiffness matrix
f = zeros(ndof,1); % Heat flux vector

for e = 1:ne
    n1 = conn(e, 1); % node id for first node in element e
    n2 = conn(e, 2); % node id for second node in element e
    n3 = conn(e, 3); % node id for first node in element e
    n4 = conn(e, 4); % node id for second node in element e

    x1 = node(n1,1); y1 = node(n1,2); % x and y coordinates for the 1st node
    x2 = node(n2,1); y2 = node(n2,2); % x and y coordinates for the 2nd node
    x3 = node(n3,1); y3 = node(n3,2); % x and y coordinates for the 3rd node
    x4 = node(n4,1); y4 = node(n4,2); % x and y coordinates for the 4th node

    global_coord = [x1, y1; x2, y2; x3 y3; x4 y4];

    %Compute Stiffness of element e
    ke = Ke(global_coord, GQ_coord, GQ_w, Vx, Vy, k_r)
    
    % locations where ke is to scatter to in the global stiffness matrix
    sctr = [ n1, n2 , n3, n4];
    
    %Add ke into global K
    K(sctr,sctr) = K(sctr,sctr) + ke;    
end

% Accounting for boundary conditions
EDIT following @bg2b's comment
f(2) = -K(1,2) * T(1);
f(3) = -K(6,3) * T(6) + -K(1,3) * T(1);
f(4) = -K(1,4) * T(1) + -K(6,4) * T(6);
f(5) = -K(6,5) * T(6);

% computing temperature
T(isol) = K(isol,isol)\f(isol);

%% FUNCTIONS
% Shape function
% takes local coord
function s = S(ni, coord)
    switch ni
        case 1
            s = 1/4 * (1-coord(1)) * (1-coord(2));
        case 4
            s = 1/4 * (1+coord(1)) * (1-coord(2));
        case 3
            s = 1/4 * (1+coord(1)) * (1+coord(2));
        case 2
            s = 1/4 * (1-coord(1)) * (1+coord(2));
    end
end

% Shape function derivative
% takes local coord as an [eta nu] pair
function ds = dS(ni, dim, coord)
    switch ni
        case 1
            if dim == 1
            ds = 1/4 * (coord(2)-1);
            else
            ds = 1/4 * (coord(1)-1);
            end
        case 4
            if dim == 1
            ds = 1/4 * (-coord(2)-1);
            else
            ds = 1/4 * (1-coord(1));
            end
        case 3
            if dim == 1
            ds = 1/4 * (coord(2) + 1);
            else
            ds = 1/4 * (coord(1) + 1);
            end
        case 2
            if dim == 1
            ds = 1/4 * (1-coord(2));
            else
            ds = 1/4 * (-coord(1)-1);
            end
    end
end

% Computation of the element's Jacobian
function J = Je(global_coord, GQ_coord)
    S = [dS(1,1, GQ_coord(1,:)) dS(2,1, GQ_coord(1,:)) dS(3,1, GQ_coord(1,:)) dS(4,1, GQ_coord(1,:));
         dS(1,2, GQ_coord(1,:)) dS(2,2, GQ_coord(1,:)) dS(3,2, GQ_coord(1,:)) dS(4,2, GQ_coord(1,:))];

    J = S * global_coord;
end

% Computation of the element stifness matrix
function ke = Ke(global_coord, GQ_coord, GQ_w, Vx, Vy, k_r)
    J = Je(global_coord, GQ_coord)
    ke = zeros(4,4);
    Jinv = inv(J);
    for i = 1:4
        for j = 1:4
            ke(i,j) = 0;
            for n = 1:4                    
                A = (Jinv(1,1) * dS(i,1, GQ_coord(n,:)) + Jinv(1,2) * dS(i,2, GQ_coord(n,:)));                  
                B = (Jinv(1,1) * dS(j,1, GQ_coord(n,:)) + Jinv(1,2) * dS(j,2, GQ_coord(n,:)));                    
                C = (Jinv(2,1) * dS(i,1, GQ_coord(n,:)) + Jinv(2,2) * dS(i,2, GQ_coord(n,:)));                    
                D = (Jinv(2,1) * dS(j,1, GQ_coord(n,:)) + Jinv(2,2) * dS(j,2, GQ_coord(n,:)));
          
                ke(i,j) = ke(i,j) + k_r * (A*B + C*D) * det(J)*GQ_w(i)*GQ_w(j);           
            end            
        end        
    end    
end

To save having to put the code for the stiffness matrix here it is:

K =[ 0.9333   -0.2333   -0.4667   -0.2333         0         0;
   -0.2333    0.9333   -0.2333   -0.4667         0         0;
   -0.4667   -0.2333    1.8667   -0.4667   -0.4667   -0.2333;
   -0.2333   -0.4667   -0.4667    1.8667   -0.2333   -0.4667;
         0         0   -0.4667   -0.2333    0.9333   -0.2333;
         0         0   -0.2333   -0.4667   -0.2333    0.9333]

So here is the cut down test code for you:

%% Input Geometry
% Local node geometry
% 4____3
% |    |
% |    |
% 1____2
% Global geometry
% 2____3____6(T=10)
% |    |    |
% |    |    |
% 1____4____5
% (T=3.5)
% Node 1 is at 3.5 degrees, and node 6 at 10 degrees.
   
%% Boundary Conditions
T = zeros(6,1);
T(1) = 3.5; % node 1 fixed at 3.5 degrees
T(6) = 10; % node 6 fixed at 10 degrees

isol = [ 2, 3, 4, 5]; % unconstrained dofs

f = zeros(ndof,1); % Heat flux vector
%global stiffness matrix
K =[ 0.9333   -0.2333   -0.4667   -0.2333         0         0;
       -0.2333    0.9333   -0.2333   -0.4667         0         0;
       -0.4667   -0.2333    1.8667   -0.4667   -0.4667   -0.2333;
       -0.2333   -0.4667   -0.4667    1.8667   -0.2333   -0.4667;
             0         0   -0.4667   -0.2333    0.9333   -0.2333;
             0         0   -0.2333   -0.4667   -0.2333    0.9333]

% Accounting for boundary conditions
% I THINK THAT'S WHERE I'M GOING WRONG
f(2) = K(2,1) * T(1); % Heat outflow
f(5) = -K(5,6) * T(6); % Heat inflow

% computing temperature
T(isol) = K(isol,isol)\f(isol);

Here is what I obtain, x is the global node index, y is the temperature. Because the only place for heat to leave the system is node 1 (as far as I can tell) I wouldn't expect to see temperature values below 3.5 degrees. enter image description here

PS: I know Matlab has a FEM toolbox, but this is really a learning exercise I've set myself.

PS 2: Feel free to leave answers in Python as I am familiar with that language as well

Upvotes: 0

Views: 228

Answers (0)

Related Questions