This post shows how to couple a heat transfer model with linear elasticity. The custom PDE equation model examines how temperature rise in break disk under breaking results in stresses and strains. Using the assumption that all variations in the rotational direction can be neglected we can use an axisymmetric formulation saving computational time compared to full 3D.

Previous posts have introduced the axisymmetric equations for stresses and strains as well as heat transfer and general transport equations. The equations used for stresses and strains have in this model been reformulated for stability (eliminating divisions by r), and also to incorporate the temperature strain relation .

eqn_T = '  r*rhod*cpd*T'' - r*Kd*(Tr_r + Tz_z) = 0';
eqn_r = '- lambda   *( c1*r^3*ur_r + r^2*wz_r ) - lambda*c2*( r^3*uz_z + r^2*wr_z  ) - lambda/nu*r^2*u_r + lambda/nu*r^2*ur_t + lambda*2*r*wz_t + lambda/nu*2*r*u_t + r^2*c3*T_r - 2*r*c3*( T_t - Tref ) = 0';
eqn_z = '- lambda*c2*( r^2*uz_r    + r*wr_r   ) - lambda   *( r^2*ur_z + c1*r*wz_z ) - lambda*2*r*u_z + r*c3*T_z = 0';

Here the temperature T is coupled to the displacements u (equivalent to u/r) and w, that is a one-way multiphysics coupling (as the temperature field is assumed to not depend on the deformation of the disk). The problem is also quasi-static in that the temperature field is time dependent, while the stresses are static (in each time step). The following commands

fea.sdim = { 'r' 'z' };
fea.dvar = { 'u' 'w' 'T' };
fea.sfun = { 'sflag2' 'sflag2' 'sflag2' };
fea.eqn  = parseeqn( { eqn_r eqn_z eqn_T }, fea.dvar, fea.sdim );

takes the FEATool strong PDE equation form above, parses it, and fills out the fea.eqn weak equation struct (that is the eqn.(). form, coef, and sfun fields). In addition two source term contributions must be manually added in order to account for the temperature coupling when a reference temperature Tref is used.

fea.eqn.f.form{1} = vertcat( fea.eqn.f.form{1}, 2 );
fea.eqn.f.form{2} = vertcat( fea.eqn.f.form{2}, 3 );
fea.eqn.f.coef{1} = vertcat( fea.eqn.f.coef{1}, '-r^2*c3*Tref' );
fea.eqn.f.coef{2} = vertcat( fea.eqn.f.coef{2}, '-r  *c3*Tref' );

This will add the weak right hand side source terms and , respectively.

Simulations have been performed with the geometry and material parameters given in reference [1] resulting in the final temperature and stress fields shown in the figure below.


FEATool Brake Disk Matlab Multiphysics Simulation Plots

The resulting temperature and stress curves on the disk surface at various times also agree well with the results computed by Nastran FEM solver in reference [1].


FEATool Brake Disk Matlab Multiphysics Simulation Curves

The complete FEATool brake disk simulation Octave and Matlab m-script model can be downloaded from the link below.

FEATool Brake Disk Stress Simulation Matlab m-script


[1] A. Adamowicz, Axisymmetric FE Model to Analysis of Thermal Stresses in a Brake Disk, Journal of Theoretical and Applied Mechanics, 53, 2, pp. 357-370, Warsaw 2015.