MATLAB Answers

matrix multiplication and transposing zeroing out problem

3 views (last 30 days)
Melissa
Melissa on 19 Jun 2013
Good Afternoon All, I am trying to calculate the stiffness matrix of a powertrain mount. I am having issues when applying my transformation matrices where they keep “zeroing” out and I cannot figure out where in my formulation the problem exists. Given a CG location [xc,yc,gc] and a mount location [x,y,z] , orientation of the mount, [phi_x,phi_y,phi_z] and principal stiffness values [kx,ky,kz] the stiffness matrix can be calculated as follows : (reference paper: http://www.mecheng.osu.edu/adl/files/adl/journal_papers/singh_article_J109.pdf)
The initial stiffness matrix is K_i=diag([kx ky kz])
and when applying to the system we must use a transformation matrix. The transformation matrix R is just the directional cosines (DCM) with respect to the CG. Thusly producing a new stiffness matrix of: K_s=R*K_r*R’
Then proceeds to use the rotation matrix (B) for the cross vector product which is: [0, –(zc-z), (yc-y); (zc-z), 0, -(xc-x); -(yc-y), (xc-x) 0];
Combining these to form the overall stiffness matrix is: K_overall=[K_s, K_s.*B’; (K_s*B’)’, B.*K_s.*B’]
The matrix zeros out at all elements of the K_overall other than the first one of K_s. Any advice on where I am going wrong would greatly be appreciated.
Here is the code I have written:
CG=[0,0,0]; Mount=[-251,173,62.7];K=[223.667,44.73,44.73]; Phi=[0,-45,0];
% User Input
% CG: Center of Gravity Position Matrix expressed as [xc,yc,zc]
% Mount: Any Mount Position Matrix expressed as [x,y,z]
% K: Stiffness Matrix of each Mount expressed as [Kx,Ky,Kz]
% Phi: Orientation Angels of each Mount expressed as [Phix,Phiy,Phiz]
%
% Output
% K_local: the local stiffness matrix
%K_s Static Mount Stiffness Matrix Caclulation
%Note: K_s is the static input stiffness by user
%Obtaining Initial Individual Stiffness Values (Initial Input)
Kx=K(1);Ky=K(2); Kz=K(3);
%Setting up initial matrix
K_i=zeros(3,3);
K_i(1,1)=Kx;
K_i(2,2)=Ky;
K_i(3,3)=Kz;
%R_i: Rotational Transformation Matrix of K_i (DCM)
%Obtaining the Initial Orientation Value of Mounts
Phi_x=Phi(1); Phi_y=Phi(2); Phi_z=Phi(3);
%Rotational Angles
angles = [Phi_z(:) Phi_y(:) Phi_x(:)];
%Direction Consine Matrix (DCM) Formulation
DCM = zeros(3,3,size(angles,1));
cang = cos( angles );
sang = sin( angles );
% Mathetmatical Representation
% [ cy*cz, cy*sz, -sy]
% [ sy*sx*cz-sz*cx, sy*sx*sz+cz*cx, cy*sx]
% [ sy*cx*cz+sz*sx, sy*cx*sz-cz*sx, cy*cx]
DCM(1,1,:) = cang(:,2).*cang(:,1);
DCM(1,2,:) = cang(:,2).*sang(:,1);
DCM(1,3,:) = -sang(:,2);
DCM(2,1,:) = sang(:,3).*sang(:,2).*cang(:,1) - cang(:,3).*sang(:,1);
DCM(2,2,:) = sang(:,3).*sang(:,2).*sang(:,1) + cang(:,3).*cang(:,1);
DCM(2,3,:) = sang(:,3).*cang(:,2);
DCM(3,1,:) = cang(:,3).*sang(:,2).*cang(:,1) + sang(:,3).*sang(:,1);
DCM(3,2,:) = cang(:,3).*sang(:,2).*sang(:,1) - sang(:,3).*cang(:,1);
DCM(3,3,:) = cang(:,3).*cang(:,2);
R_i=DCM;
%K_s: Static Stiffness Matrix of Mount
K_s=R_i.*K_i.*R_i.';
%Rotation Matrix for the cross-vector product (B) Formulation
B=zeros(3,3);
B(1,2)=-(CG(3)-Mount(3));
B(1,3)=CG(2)-Mount(2);
B(2,1)=CG(3)-Mount(3);
B(2,3)=-(CG(1)-Mount(1));
B(3,1)=-(CG(2)-Mount(2));
B(3,2)=CG(1)-Mount(1);
B;
%Local Stiffness Matrix Formulation
K_local_11=K_s;
K_local_12=K_s.*B.';
K_local_21=(K_s.*B.').';
K_local_22=B.*K_s.*B.';
K_local=[K_local_11, K_local_12; K_local_21, K_local_22]

Accepted Answer

Iain
Iain on 19 Jun 2013
Your problem is that you're using ELEMENTWISE multiplication. Use normal, matrix multiplication.
Also, if you have a DCM "T", and an initial set of co-ordinates, "X" (column vector), you simply
X_new_axis = T * X;
I'm not sure why you're using it as "ks = T*K*T' ", but I think you might be thinking of quaternions.

More Answers (1)

Sean de Wolski
Sean de Wolski on 19 Jun 2013
If you put a breakpoint on the first line and then step through this with the debugger (i.e. dbstep or the dbstep button), what is the first line where the defined variable is not what you would expect?
  3 Comments
Melissa
Melissa on 19 Jun 2013
Ah okay. Yeah that would do it. Wonder where I went wrong with my derivation then. Thanks a lot Sean de Wolski :)

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!