rayleigh_coef.m

结构力学中的有限元例子,包含了7个分类文件夹

function [alph, bet] = rayleigh_coef (K,M,damp,N)

% identify coeff (alpha, beta) for Rayleigh damping matrix 
% INPUT:
%         K     -   stiffness matrix
%         M     -   mass matrix
%         damp  -   damping for 1st & N modes (vector)
%         N     -   the 1-st and N-th modes will be taken for coeff calculation
% OUTPUT:
%         alph  -   coeff for Rayleigh damping matrix
%         bet   -   coeff for Rayleigh damping matrix
% uses: eigFEM_s.m
% 2004

% "N" smallest eigenvalues and eigenvectors for stiffness & mass: K, M
[ E_Vec, eig_Val] = eigFEM_s (K, M, N); % eig_Val gives vector of omega from high to low

bet  = (2*damp(2)*eig_Val(1)-2*damp(1)*eig_Val(N))/(eig_Val(1)^2-eig_Val(N)^2);
alph = 2*damp(2)*eig_Val(1) - bet*eig_Val(1)^2;


% MDOF system under externally applied time-dependent force P:
%        M*\ddot X + C*\dot X + K*X = P                                                (1)
% by orthogonal transformation:
%        eig_v'*M*eig_v*\ddot x + eig_v'*C*eig_v*\dot x + eig_v'*K*eig_v*x = eig_v'*P  (2)
% (1) is reduced to "n" uncoupled equations of the form:
%        \ddot x_i + 2*damp_i*omega_i*\dot x_i + omega_i^2*x_i = P_i(t)                (3)
% where:
% x     - displacements in the transformed coordinates,
% damp  - damping ratio in uncoupled mode,
% omega - natural frequency of the system,
% P(t)  - modified force vector in transformed coordinates,
% eig_v - normalized eigen vectors of the system.
%
% orthogonal transformation (2) is valid only when C is some function of M & K. The
% orthogonal transformation of damping in (2) reduces C to:
%       eig_v'*C*eig_v = [a+b*omega_1^2        0                0          ...;
%                               0          a+b*omega_2^2        0          ...;
%                               0              0            a+b*omega_3^2  ...];
% which is a system of equations:
%     2*damp_i*omega_i = a + b*omega_i^2                                               (4)
% It is obvious that "a" & "b" could not satisfy all modes. Solving (4) for damping:
%     damp_i = a/(2*omega_i) + b*omega_i/2                                             (5)
% (5) gives non-linearity of damping for low frequencies (first term in (5) dominates) 
% and for higher frequencies damping is linear (2nd term in (5) dominates).
%
% First few modes usually give significant mass participation.
% Engineering structures usually give 95% of mass participation for the first 3 modes.
% It is realistic to assume that damping ratio is linearly proportional to frequency.
%     damp_i = (damp_m-damp_1)/(omega_m-omega_1)*(omega_i-omega_1)+damp_1              (6)
% where: 
% damp_m - damping ratio for "m"th significant mode,
% damp_1 - damping ratio for the first mode,
% damp_i - damping ratio for the "i"th mode (i<=m).
% (6) could be used to extrapolate for interval m<i<2.5*m:  ... (omega_(m+i)-omega_m)...  .

% Reference: Indrajit Chowdhury & Shambhu P. Dasgupta. Computation of Rayleigh 
%            Damping Coefficients for Large Systems.