thermal induced in simply supported beam

4 次查看（过去 30 天）
mohanad al-tahrawi
mohanad al-tahrawi 2023-1-19
回答： Raag 2025-7-4
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)');

请先登录，再回答此问题。

回答（1 个）

Raag
Raag 2025-7-4
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

类别

在 Help Center 和 File Exchange 中查找有关 General Physics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by