This tutorial illustrates how to couple a heat transfer model with linear elasticity by manually defining the multiphysics equations with the FEATool PDE syntax. The custom model examines how temperature rise in a brake disk under braking results in heat induced displacements with resulting increased stresses and strains. Using the assumption that all variations in the rotational direction can be neglected it is possible to use an axisymmetric formulation in cylindrical coordinates, saving a significant amount of computational time and resources compared to full 3D simulations.

Previous tutorial articles have introduced the axisymmetric equations for stresses and strains as well as heat transfer and general transport equations. The temperature equation for heat conduction in axisymmetric form can simply be defined as

```
eqn_T = 'r*rho*cpd*T'' - r*K*(Tr_r + Tz_z) = 0';
```

where *rho* is the density, *cp* the heat capacity, *K* the thermal
conductivity, and *r* and *z* the space dimension
coordinates. (Although this tutorial uses the MATLAB m-script and
custom equation approach to set up and run the model, the FEATool GUI
with predefined axisymmetric physics modes can be used as well.)

The PDE equations used for calculating the displacements *u* and *w*
are here reformulated for stability by multiplying the *r*-equation
with *r ^{2}* and substituting and solving for

*u/r*instead of just

*u*, and additionally multiplying the

*z*-equation with

*r*. This eliminates divisions by

*r*at the center

*r = 0*line. Moreover, the equations also here incorporate the temperature strain relation $\epsilon_{th} = \alpha(T-T_{ref})$. The resulting axisymmetric stress-strain equations read

```
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 - T_ref )_r - 2*r*c3*( T - Tref )_t = 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 - T_ref )_z = 0'];
```

with the equation coefficients $\lambda = E\frac{\nu}{(1+\nu)(1-2\nu)}$, $c_1 = \frac{1-\nu}{\nu}$, $c_2 = \frac{1-2\nu}{2\nu}$, and $c_3 = \alpha\frac{E}{1-2\nu}$.

When a constant reference temperature *T_ref* is used implicit
bilinear forms involving the temperature *T*, as well as the weak
equation linear right hand side (RHS) source terms $-\int
r^2c_3T_{ref}\frac{\partial \phi_u}{\partial r} d\Omega$ and $-\int
rc_3T_{ref}\frac{\partial \phi_u}{\partial z} d\Omega$ will
automatically be added.

The temperature *T* is here coupled to the displacements *u* (here
redefined as *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 command statements

```
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 );
```

take 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).

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

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].

The complete custom PDE FEATool brake disk simulation MATLAB m-script model using the pre-defined axisymmetric stress-strain physics mode is available as the ex_axistressstrain4 model in the FEATool examples directory.

[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.