thermal induced in simply supported beam

i want to write matlab code (FEM) for thermal induced vibration isotropic simply supported beam with uniform cross section

i write this code but its didnt work can anyone help me to detect the error, its give me :

Error in Untitled3 (line 30)

K(i:i+1,i:i+1) = K(i:i+1,i:i+1) + k;

the code is :

% Define beam properties

E = 30e6; % Young's modulus (Pa)

I = 1e-4; % Moment of inertia (m^4)

L = 1; % Length of beam (m)

alpha = 12e-6; % Coefficient of thermal expansion (1/K)

DeltaT = 10; % Temperature change (K)

% Discretize the beam into finite elements

n = 100; % Number of elements

dx = L/n; % Element size

x = 0:dx:L; % Discretized x-coordinates

% Define the stiffness matrix

k = (E*I/dx^3)*[12, 6*dx, -12, 6*dx; 6*dx, 4*dx^2, -6*dx, 2*dx^2; -12, -6*dx, 12, -6*dx; 6*dx, 2*dx^2, -6*dx, 4*dx^2];

% Define the thermal vector

T = alpha*DeltaT*[-1; -1; 1; 1];

% Define the mass matrix

m = (dx/420)*[156, 22*dx, 54, -13*dx; 22*dx, 4*dx^2, 13*dx, -3*dx^2; 54, 13*dx, 156, -22*dx; -13*dx, -3*dx^2, -22*dx, 4*dx^2];

% Define the boundary conditions

u0 = 0; % Displacement at first node

un = 0; % Displacement at last node

% Assemble the global stiffness and thermal vectors

K = zeros(n+1);

Tg = zeros(n+1,1);

for i = 1:n

K(i:i+1,i:i+1) = K(i:i+1,i:i+1) + k;

Tg(i:i+1) = Tg(i:i+1) + T;

end

% Modify the matrix to account for boundary conditions

K(1,1) = K(1,1) + 1;

K(n+1,n+1) = K(n+1,n+1) + 1;

% Solve for the nodal temperatures

Tn = K\Tg;

% Solve for the nodal displacements

[V,D] = eig(K,m);

[~,idx] = min(abs(diag(D)));

wn = sqrt(D(idx,idx));

fn = V(:,idx);

% Plot the results

figure;

plot(x,Tn);

xlabel('x (m)');

ylabel('Temperature (K)');

figure;

plot(x,fn);

xlabel('x (m)');

ylabel('Displacement (m)');

Answer 1

[Raag]

Hi mohanad,

As per my understanding, the issue you are facing arises from how the global stiffness matrix K is being assembled in your FEM code for thermal-induced vibration of a simply supported beam.

The error occurs because your element stiffness matrix k is a 4×4 matrix, which is appropriate for beam elements with 2 degrees of freedom (DOFs) per node (transverse displacement and rotation). However, your global matrix K is initialized as if each node only has 1 DOF, leading to a mismatch in dimensions during assembly.

To fix the issue, we need to update the size of the global matrices to account for 2 DOFs per node. We can do this in the following way:

ndof = 2 * (n + 1); % Total DOFs 

K = zeros(ndof, ndof); 

And we need to update the indexing in your assembly loop to match the 4 DOFs per element in the following manner:

for i = 1:n 

    idx = 2*(i-1) + (1:4); % DOF indices for the current element 

    K(idx, idx) = K(idx, idx) + k; 

end 

The output will look like:

This ensures that the element matrices are correctly mapped into the global system.

For more details, refer:

1. Finite element method basics: https://www.mathworks.com/help/pde/ug/basics-of-the-finite-element-method.html

2. structural mechanics: https://www.mathworks.com/help/pde/structural-mechanics.html