Eigen Function for modes of natural frequency

32 views (last 30 days)
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)))
ans = 6×1
2.8545 21.4957 63.7919 135.2966 242.8566 957.7791
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
Torsten
Torsten on 16 Feb 2024
Edited: Torsten on 16 Feb 2024
Seriously: How should anybody not involved in the physical background of your problem be able to answer this ?
And your RelTol and AbsTol values are much too low as to any solver could be able to satisfy these error tolerances.
Syed Abdul Rafay
Syed Abdul Rafay on 17 Feb 2024
I have already tried different tolerance values but result is same. I have tried more than 100 different combinations. The hand calculations give correct 4th 5 th values but when transferred to code they change only first 3 values are correct rest diverge. The problem here is not the equations it has to do something with matlab implementations that's why I am asking question here.

Sign in to comment.

Answers (1)

Rupesh
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!

Products


Release

R2022b

Community Treasure Hunt

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

Start Hunting!