The following post will break down and explain the steps used to create the images and video for the axisymmetric stress-strain brake disk analysis model. By using the Octave or Matlab command line and exported finite element fea data we can create custom postprocessing and visualizations that is not possible in the FEATool GUI.

The first step is to run the FEATool model script file to generate and return the finite element struct fea and solution stored in fea.sol.u (in this case the automatic pre-defined plots can be turned off, iplot=0, since we will do postprocessing manually)

fea = axi_stress_strain_heat_brake_disk( 'iplot', 0 );

Next we prescribe the model constants and parameters used in the postprocessing

deltad = 5.5e-3;     % Disk thickness.
rd     = 66e-3;      % Disk inner radius.
rp     = 75.5e-3;    % Pad inner radius.
Rd     = 113.5e-3;   % Disk outer radius.

To continue, the solution vector for all time steps is saved in a temporary variable u (the solution at each time step corresponds to a column in u). The parameters cmax and cmin are also calculated which later are used to fix the colormap range (otherwise we would not be able to see much difference between hot and cold maxima)

u    = fea.sol.u;
cmax = max(u(:)) - 273.15;
cmin = 0;

Next, coordinates to plot and visualize the countours and boundaries of the brake disk are created (here the assumption is that the disk is centered at the origin with the radial coordinate spanning the x-z plane)

% Points around the circumference.
npth = 72;
th   = linspace( 0, 2*pi, npth );

% Inner coordinates.
x    = rd*cos(th);
z    = rd*sin(th);

% Outer coordinates.
xx   = Rd*cos(th);
zz   = Rd*sin(th);

y    = deltad*ones(size(x));

We also use the ringgrid function to generate coordinates around the disk in the grid struct g and evaluation points in rz. These will be used to interpolate the temperature field from axisymmetric coordinates to three dimensions (the ringgrid can be visualized by itself with using the plotgrid(g) command)

npr  = 25;
g    = ringgrid( npr-1, npth-1, rd, Rd );
r    = sqrt( g.p(1,:).^2 + g.p(2,:).^2 );
rz   = [rd+sqrt(eps) rd+(Rd-rd)/2 Rd-sqrt(eps);deltad/2*[1 1 1]];

Now we can begin by creating a figure to plot in, an empty array for the stresses at the three points (inner radius, outer radius, and the mid point), and starting the postprocessing loop

st = zeros(3,1);
for i=1:size(u,2)

We use subplot to hold two horizontal images. First is the 3D temperature field, starting with drawing the black outline of the brake disk

  hold on

  plot3(x,  y, z,'k-')
  plot3(x, -y, z,'k-')
  plot3(xx, y,zz,'k-')

The i:th solution is selected as the only solution in fea.sol.u (since evalexpr per default uses the solution at the last time step), after which the temperature is evaluated along the radius at the disk center line coordinate (z=0).

  fea.sol.u = u(:,i);
  T = evalexpr( 'T-273.15', [r;zeros(size(r))], fea );

We use the fact that the solution is axisymmetric to extrapolate this data all along the brake disk using the annular grid struct g created earlier. The Matlab patch command is used to plot the x-z temperature field

  h = patch( 'faces',     g.c', ...
             'vertices',  [g.p(1,:)' zeros(size(g.p,2),1) g.p(2,:)'], ...
             'facevertexcdata', T, ...
             'facecolor', 'interp', ...
             'linestyle', 'none' );

Furthermore, we also plot the temperature for the original patch in the positive and mirrored (negative) z-directions

  postplot( fea, 'surfexpr', 'T-273.15', 'colorbar', 'off' )
  fea.grid.p(2,:) = -fea.grid.p(2,:);
  postplot( fea, 'surfexpr', 'T-273.15', 'colorbar', 'off' )
  fea.grid.p(2,:) = -fea.grid.p(2,:);

Lastly, we set a title, fix the color map range, and choose a suitable viewing angle

  title( 'Temperature' )
  caxis( [cmin cmax] )
  axis off
  axis tight

For the second subplot we want to follow the von Mieses stress as a function of time for the three points in rz. Just as for the temperature, the stress can be evaluated with the evalexpr command where the stress expression svm has been defined earlier in the model fea.const field

  hold on
  grid on
  if( i>1 )   % Skip first zero solution
    st = [st evalexpr( '(svm)*1e-6', rz, fea )];

The stress is plotted as a function of time, with additional text and legends added for identification

    plot( tlist(1:i), st(1,:) )
    plot( tlist(1:i), st(2,:) )
    plot( tlist(1:i), st(3,:) )

    text( tlist(i), st(1,end), 'r = rd' )
    text( tlist(i), st(2,end), 'r = (rd+Rd)/2' )
    text( tlist(i), st(3,end), 'r = Rd' )

    axis( [0 tmax 0 max(st(:))] )
  xlabel( 'Time [s]' )
  ylabel( 'von Mieses stress [MPa]' )

Finally, we force drawing with the drawnow command, optionally write out and save a jpeg image file, and close the loop

  print( '-r200', '-djpeg', sprintf('img%03i.jpg',i) )

After processing has finished the images can easily be converted to a video with a suitable tool, such as for example ffmpeg

ffmpeg -i img%03d.jpg -c:v libx264 -vf "fps=25,format=yuv420p,crop=1600:1198:0:0" out.mp4

From the video we can see that the temperature is the highest on the outer two-thirds of the disk consistent with where the break pad is acting. In contrast the maximum stress is achieved at the inner radius r=rd around t=3, and is significantly higher than at the two points measured further out on the disk.