!!!At the very bottom, there is a much simpler version of this question!!!
I am trying to model a beam using FEM in Matlab. I understand how to use ode23 for single dimentional problems however, my problem is 10 dimentional and second order. It looks like [M]*(xdotdot)+[K]*x=Q where M and K are 10×10 matrices. I think I have my code close to being right but I'm not quite sure what is wrong. When I debug the code it stops when It gets to the ode function at the bottom. Note that it is calling function sys.m
Here is my code:
clc; clearformat shortEt=1:.01:10; % time
% -------------------------------------------------------
%{
% User input
l = 3; % length in inchesEI1 = 5*10^6; % EIEI2 = 2*EI1; % EIEI3 = EI1; % EIEI4 = EI2; % EIk1 = 200; % spring constant 1k2 = k1; % spring constant 2m = 268.2; % mass per unit lengthpo = 100; % distributed loadd1 = 2; % distance to Pd2 = 1; % distance to spring on Element 3d3 = 1.5; % distance to spring on Element 4P = 75; % point load%}
% -------------------------------------------------------% -------------------------------------------------------% User input% Validation code
l = 3; % length in inches
EI1 = 5*10^6; % EI
EI2 = 2*EI1; % EIEI3 = EI1; % EIEI4 = EI2; % EIk1 = 0; % spring constant 1
k2 = k1; % spring constant 2
m = 268.2; % mass per unit length
po = 100; % distributed load
d1 = 2; % distance to P
d2 = 1; % distance to spring on Element 3
d3 = 1.5; % distance to spring on Element 4
P = [75]'; % point load
zeta = .5;% -------------------------------------------------------% -------------------------------------------------------% Phi function handles
phi1 = @(zeta) 1-3*zeta^2+2*zeta^3;phi2 = @(zeta) l*zeta-2*l*zeta^2+l*zeta^3;phi3 = @(zeta) 3*zeta^2-2*zeta^3;phi4 = @(zeta) -l*zeta^2+l*zeta^3;Phi = @(zeta)[... phi1(zeta)*phi1(zeta) phi1(zeta)*phi2(zeta) phi1(zeta)*phi3(zeta) phi1(zeta)*phi4(zeta);... phi2(zeta)*phi1(zeta) phi2(zeta)*phi2(zeta) phi2(zeta)*phi3(zeta) phi2(zeta)*phi4(zeta);... phi3(zeta)*phi1(zeta) phi3(zeta)*phi2(zeta) phi3(zeta)*phi3(zeta) phi3(zeta)*phi4(zeta);... phi4(zeta)*phi1(zeta) phi4(zeta)*phi2(zeta) phi4(zeta)*phi3(zeta) phi4(zeta)*phi4(zeta) ];%PPhi = Phi(zeta)
%}% This function needs to be moved to where it is needed
% -------------------------------------------------------% -------------------------------------------------------% Displacement function handles
Qpo = @(po,l) [po*l/2 po*l^2/12 po*l/2 -po*l^2/12]';QP = @(zeta,P) [-P*phi1(zeta) -P*phi2(zeta) -P*phi3(zeta) -P*phi4(zeta)]';Qspring3 = @(k1,zeta) k1*Phi(zeta);Qspring4 = @(k1,zeta) k1*Phi(zeta);% This function needs to be moved to where it is needed% -------------------------------------------------------% -------------------------------------------------------% Calculation simplification
a=zeros(4,6);aa=zeros(4);aaa=zeros(4,2);I4=eye(4);A1=[I4 a];A2=[aaa I4 aa];A3=[aa I4 aaa];A4=[a I4];% -------------------------------------------------------% -------------------------------------------------------% load transformations
% p
% P
% -------------------------------------------------------% -------------------------------------------------------% Mass matrix, Stiffness matrix
M1 = (m*l/420)*[... 156 22*l 54 -13*l; 22*l 4*l^2 13*l -3*l^2 54 13*l 156 -22*l -13*l -3*l^2 -22*l 4*l^2];M2 = M1; % note, this can be different
M3 = M1; % note, this can be differentM4 = M1; % note, this can be differentK1 = (EI1/l^3)*[... 12 6*l -12 6*l 6*l 4*l^2 -6*l 2*l^2 -12 -6*l 12 -6*l 6*l 2*l^2 -6*l 4*l^2];K2 = (EI2/l^3)*[... 12 6*l -12 6*l 6*l 4*l^2 -6*l 2*l^2 -12 -6*l 12 -6*l 6*l 2*l^2 -6*l 4*l^2];K3 = (EI3/l^3)*[... 12 6*l -12 6*l 6*l 4*l^2 -6*l 2*l^2 -12 -6*l 12 -6*l 6*l 2*l^2 -6*l 4*l^2];K4 = (EI4/l^3)*[... 12 6*l -12 6*l 6*l 4*l^2 -6*l 2*l^2 -12 -6*l 12 -6*l 6*l 2*l^2 -6*l 4*l^2];% -------------------------------------------------------% calculations
K3 = K3 + Qspring3(k1,l/d2); % Khat
K4 = K4 + Qspring4(k1,l/d3); % KhatQpo4 = Qpo(po,l);QP4 = QP(l/d1,P);syms Q11 Q21 Q31 Q41syms Q12 Q22 Q32 Q42syms Q13 Q23 Q33 Q43syms Q14 Q24 Q34 Q44Q1 = [Q11 Q21 Q31 Q41]';Q2 = [-Q31 -Q41 Q32 Q42]';Q3 = [-Q32 -Q42 Q33 Q43]';Q4 = [-Q33 -Q43 0 0]'+QP4+Qpo4;Q11 = A1'*Q1;Q22 = A2'*Q2;Q33 = A3'*Q3;Q44 = A4'*Q4;Q = Q11+Q22+Q33+Q44;Q(1) = 0; Q(2) = 0;% -------------------------------------------------------% -------------------------------------------------------% global calculations
M = A1'*M1*A1+A2'*M2*A2+A3'*M3*A3+A4'*M4*A4;K = A1'*K1*A1+A2'*K2*A2+A3'*K3*A3+A4'*K4*A4;Minv = inv(M);% -------------------------------------------------------%qstatic = inv(K)*Q
time = (0:.001:22.5)';x0 = [0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0]; %[position,velocity]
[x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,tt]=ode23(@(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,tt) sys(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,tt,Minv,K), time, x0);
Then my sys.m looks like:
function [x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,tt] = sys(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,tt,Minv,K)f = zeros(20,1);f(1) = x(2);f(2) = x(3);f(3) = x(4);f(4) = x(5);f(5) = x(6);f(6) = x(7);f(7) = x(8);f(8) = x(9);f(9) = x(10);f(10) = x(11); delta = -Minv*K*[x(12) x(13) x(14) x(15) x(16) x(17) x(18) x(19) x(20) x(21)]'%+Minv*Q;
f(11) = delta(1)f(12) = delta(2)f(13) = delta(3)f(14) = delta(4)f(15) = delta(5)f(16) = delta(6)f(17) = delta(7)f(18) = delta(8)f(19) = delta(9)f(20) = delta(10)end
Here is the easier version of the code:
clc; clearfiguretime = (0:.001:22.5)';M = [1 2; 3 4];K = [5 6; 7 8];x0 = [0;0;0;0]; %[position,velocity][t1,x1,x2,x3,x4]=ode23(@(t1,x) trick(t1,x,M,K), time, x0);plot(t1,x(:,1));xlabel('Time(sec)');ylabel('Displacement');title('Stepped Response(Underdamped)');
Here is the code that it calls:
%2x2 sys
function f = trick(t,x,M,K)f = zeros(4,1);f(1) = x(2);f(2) = x(3);MinvnegK = -inv(M)*K;delta = MinvnegK*[x(4) x(5)]';f(3) = delta(1);f(4) = delta(2)
Best Answer