Eigen Function for modes of natural frequency
32 views (last 30 days)
Show older comments
clear all
clc
syms lambda
syms x
syms i
rhoEpoxy = 1200;
Em = 3e+09;
R=5;
E_L=70e+9;
E_r=200e+9;
rho_L=2702;
rho_r=5700;
nu=3;
A_K=zeros(R+1,R+1);
A_M=zeros(R+1,R+1);
G=zeros(R+1,R+1);
H1_K=zeros(R+1,R+1);
H1_M=zeros(R+1,R+1);
H2_M=zeros(R+1,R+1);
for ri=1:R+1
r = ri-1;
for mi=1:R+1
m = mi-1;
fun1 = @(ksei1) (ksei1.^(r+m)).*((1-(exp(nu*ksei1)-1)./(exp(nu)-1))+E_r/E_L*(exp(nu*ksei1)-1)./(exp(nu)-1));
%G(mi,ri)=integral(fun1,0,1);
fun_h1=@(ksei2,s2) (-2*nu/(exp(nu)-1)*(E_r/E_L-1)*exp(nu*s2)).*(s2.^r).*(ksei2.^m)+...
(nu^2/(exp(nu)-1)*(E_r/E_L-1)*exp(nu*s2)).*(s2.^r).*(ksei2.^(m+1))-...
(nu^2/(exp(nu)-1)*(E_r/E_L-1)*exp(nu*s2)).*(s2.^(r+1)).*(ksei2.^m);
fun_h1_lambda=@(ksei3,s3) (-1/6*((ksei3-s3).^3).*((1-(exp(nu*s3)-1)./(exp(nu)-1))+...
rho_r/rho_L*(exp(nu*s3)-1)./(exp(nu)-1))).*(s3.^r).*(ksei3.^m);
fun_h2_lambda=@(ksei4,s4) (1/6*(ksei4.^3).*((1-(exp(nu*s4)-1)./(exp(nu)-1))+...
rho_r/rho_L*(exp(nu*s4)-1)./(exp(nu)-1))-1/2*(ksei4.^2).*s4.*((1-(exp(nu*s4)-1)./(exp(nu)-1))+...
rho_r/rho_L*(exp(nu*s4)-1)./(exp(nu)-1))).*(s4.^r).*(ksei4.^m);
options = {'RelTol', 1e-22, 'AbsTol', 1e-24};
G(mi,ri) = integral(fun1, 0, 1, options{:});
H1_K(mi,ri) = integral2(fun_h1, 0, 1, 0, @(ksei2) ksei2, options{:});
H1_M(mi,ri) = integral2(fun_h1_lambda, 0, 1, 0, @(ksei3) ksei3, options{:});
H2_M(mi,ri) = integral2(fun_h2_lambda, 0, 1, 0, 1, options{:});
A_K(mi,ri)=G(mi,ri)+H1_K(mi,ri);
A_M(mi,ri)=H1_M(mi,ri)+H2_M(mi,ri);
end
end
sort(sqrt(eig(A_K,-A_M)))
For some reason my 4th, 5th 6th mode values are not correct they don't match the results. My first 3 mode values for vibration are correct adn they match the research paper. after 3rd mode no matter what boyndary condition I am using they don't match after 3rd mode. Can someone tell me why.
3 Comments
Answers (1)
Rupesh
on 25 Mar 2024
Edited: Rupesh
on 26 Mar 2024
Hi Syed Abdul Rafay,
It looks like the problem you're having with finding the higher mode values might be because of some tricky parts in the math calculations MATLAB does. To address the discrepancy in the 4th mode value in your vibration analysis problem, you can focus on two main modifications: refining the numerical integration technique and enhancing the accuracy of the eigenvalue problem solver. Below are some thoughts and potential adjustments that might guide you towards resolving the issue:
1] Numerical Integration Precision
Here the goal is to improve the precision for numerical inetgration, especially for integrands that exhibit complex behavior over the integration range. A useful strategy is to segment the integration range into parts where the integrand behaves more uniformly. This approach allows for more accurate integration over each segment, which is particularly beneficial for capturing the nuances of higher modes.
% Assuming 'fun1' is your integrand and it behaves differently across the range
splitPoint = 0.5; % This is an illustrative value; adjust based on your integrand's behavior
options = {'RelTol', 1e-6, 'AbsTol', 1e-8}; % Fine-tune these tolerances as needed
integralValue = integral(fun1, 0, splitPoint, options{:}) + integral(fun1, splitPoint, 1, options{:});
2] Eigenvalue Solver Accuracy
The choice and configuration of the eigenvalue solver are crucial for accurately determining the vibration modes. When dealing with large matrices or seeking a subset of eigenvalues, the eigs function can offer a more efficient and potentially more accurate approach than the standard eig function, especially for higher modes.
% Assuming A_K and A_M are your stiffness and mass matrices, respectively
numModes = 6; % Number of modes you're interested in
opts.tol = 1e-6; % Adjust the tolerance to improve accuracy
[eigenvectors, eigenvalues] = eigs(A_K, A_M, numModes, 'smallestabs', opts);
sorted_eigenvalues = sort(sqrt(diag(eigenvalues)));
3] Matrix Conditioning
The condition of the matrices used in the eigenvalue problem significantly affects the accuracy of the results. Ensuring that the stiffness and mass matrices (A_K and A_M) are well-conditioned .This might involve revisiting how these matrices are constructed and ensuring they are symmetric when theoretically expected to be so.
% Check for symmetry and condition number of A_K and A_M
if ~isequal(A_K, A_K.') || ~isequal(A_M, A_M.')
warning('Matrices A_K or A_M are not symmetric. Check their construction.');
end
if cond(A_K) > 1e10 || cond(A_M) > 1e10
warning('One of the matrices is poorly conditioned. Consider revising its construction.');
end
You can also refer to below document for the usage of above-mentioned functions .
Hope this helps!
0 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!