Calculating Stiffness of T-Base Lathe

In this notebook we walk through the calculations transformation of stiffness results for my TBase lathe using MATLAB functions from my error budgeting program

%Execute this code from working directory FEEB/notebook ONCE so that the working directory changes to src
%cd ..\..
%cd FEEB\src

Model Structure

Currently the model is stored as a number of structures in MATLAB. As MATLAB is primarily a scripting language, the program is very flat in its structure with several dozen scripts containing functions that act on the structures. This methodology is prone to error and I will soon refactor the code in a way that takes advantage of the object-oriented programming features available in MATLAB.

clc; clear all;
  Name            Size            Bytes  Class     Attributes

  Beams           1x11            90536  struct              
  CS              1x10            24726  struct              
  DOFs            1x72            78102  struct              
  Elements        1x11            93544  struct              
  Joints          1x12             9026  struct              
  K              72x72            41472  double              
  Nodes           1x12             6688  struct              
  Reordered       1x1             44960  struct              
  Results         1x1             90496  struct              

Example 1: Calculating Stiffness of T-Base Lathe

To demonstrate the utility of the MATLAB program I created the following example. Here we want to determine the stiffness of the tool tip of a TBase lathe when a force is applied to the tip of the TBase lathe. Calling ‘SolveStiffness’ automatically sets up 36 loading scenarios to determine the stiffness matrix of the displacement at some node relative to a force applied at another node. In each load scenario, a prescribed displacement boundary condition is applied in the direction of interest, and the reaction forces at the node of interest are recorded. Dividing the reaction force by the applied force gives a measure of the stiffness at the point of interest.

%Perform stiffness calculation study
%Measure stiffness at joint1 relative to force applied at joint2


[StiffnessMatrix,~] = SolveStiffness(Joints,Beams,jointapply,jointmeasure);

%Measured the diagonal components so those are the ones we need for now

%Format output so it is easy to read
for i=1:6
    if i<=3
        div=1000; %to account for Rad->mRad
    disp([names(i),sprintf(' = %.3f ',StiffnessMatrix(i)/div),unit_out])


    "Kxx"    " = 93.805 "    "n/mm"

    "Kyy"    " = 562.939 "    "n/mm"

    "Kzz"    " = 99.396 "    "n/mm"

    "Krxx"    " = 1071.483 "    ""

    "Kryy"    " = 199.716 "    ""

    "Krzz"    " = 1197.025 "    ""

Matrix transformation function operates on elements so we transform the node we are using to represent the TBase model into a simple spring with only diagonal elements.

StiffnessMatrix = diag(repmat(StiffnessMatrix,1,2));

Example 2: Stiffness Coordinate System Transformation

We have successfully calculated the stiffness matrix of the TBase lathe, but the actual stiffness was measured with respect to a reference frame that was more convenient for the metrology setup. In order to compare results we must transform the stiffness estimate from the model coordinate system to the metrology coordinate system. Both Coordinates are shown in the plot below.

%Generate transformation matrix
tx=0;ty=0;tz=0;rx=0; ry=-90; rz=0;
RotationHTM = RotationTransform([tx,ty,tz,rx,ry,rz])

%Plot related
ModelHTM=eye(4); ModelHTM(:,4)=ones(4,1);
RotationHTM_Tensor(:,:,1) = ModelHTM;
RotationHTM_Tensor(:,:,2) = RotationHTM;
text(0.9,0.9,0.9,"Model Orientation")
text(-0.1,-0.1,-0.1,"Measured Orientation")
axis equal
RotationHTM =

     0     0    -1     0
     0     1     0     0
     1     0     0     0
     0     0     0     1

Illustration of the orientation of the MATLAB model CS and the metrology CS
%Perform coordinate transformation
AlignedMatrix = LocalToGlobal(StiffnessMatrix,RotationHTM);
%Extract diagonal components
AlignedMatrix = AlignedMatrix(1:6,1:6);
StiffnessComponents = diag(AlignedMatrix);
stringvals=sprintf('%.3f \n',StiffnessComponents);

for i=1:6
    if i<=3
        div=1000; %to account for Rad->mRad
    disp([names(i),sprintf(' = %.3f ',StiffnessComponents(i)/div),unit_out])


    "Kxx"    " = 99.396 "    "n/mm"

    "Kyy"    " = 562.939 "    "n/mm"

    "Kzz"    " = 93.805 "    "n/mm"

    "Krxx"    " = 1197.025 "    ""

    "Kryy"    " = 199.716 "    ""

    "Krzz"    " = 1071.483 "    ""