# Why do I get zero Thermal Stress? Am I using the PDE Toolbox incorrectly?

3 views (last 30 days)

Show older comments

I am attempting to determine the Von Mises Stresses that occur in a long tunsten tube when applied to a thermal load. I don't know if it is because of the matieral properties of Tungsten coupled with the thin walls of the tube, or if it is something I am during incorrectly in my code, but everytime I solve the structural model given the results of the thermal model, I get ZERO STRESSES. Not even small numbers...ZERO.

Here is my code:

tmodel = createpde('thermal','transient'); % Create a Transient Thermal Model to out thermal results to Structural Model

%% Geometry

D = 2.5; % Inner Diameter of Flow Tube [mm]

thick = 0.125; % Thickness of Flow Tube [mm]

rin = D/2; % Inside Radius of Flow Tube [mm]

rout = rin + thick; % Outside Radius of Flow Tube [mm]

L = 1318; %Length of Flow Tube [mm]

gm = multicylinder([rin,rout],L,'void',[1,0]); % Create flow tube geometry

msh = generateMesh(tmodel,'Hmax',0.5); % Generates a mesh using tetrahedrons no larger than Hmax

%% Thermal Properties

% Specific Heat f(T)

SHfun =@(location,state) (21.868372 + 8.068661e-3*state.u - 3.756196e-6*state.u.^2 + 1.075862e-9*state.u.^3 + (1.406637e4./state.u.^2))./183.84; % [J/gK]

% Thermal Conductivity f(T)

TCfun =@(loaction,state) (149.441 - 45.466e-3*state.u + 13.193e-6*state.u.^2 - 1.484e-9*state.u.^3 + (3.866e6./state.u.^2))./1000; % [W/mmK]

% Mass Density f(T)

MDfun =@(location,state) (19.25 - 2.66207e-4*(state.u - 293.15) - 3.0595e-9*(state.u - 293.15).^2 - 9.5185e-12*(state.u - 293.15).^3)./1000; % [g/mm3]

% Apply

thermalProperties(tmodel,'ThermalConductivity',TCfun, ...

'MassDensity',MDfun, ...

'SpecificHeat',SHfun); % Apply Thermal Properties

%% Thermal Initial Conditions

% For the sake of this questions, lets say the thermal profile is as follows:

T0_fun =@(location) location.z.^2 % Non linear Temperature profile applied axially

thermalIC(tmodel,T0_fun)

%% Solve

tlist = [0 0.1]

Tresults = solve(tmodel,tlist); % Solves Thermal model for arbitary time steps

The point of the code up to this point is to get the thermal profile that I want to apply to my structural model into the proper PDE format of Tresults. Now that that is down, I apply it to the Structural model.

smodel = createpde('structural','static-solid'); % Create Structural Model

smodel.Geometry= gm;

smodel.Mesh = msh;

Unfortunatly it doesn't allow the Young's Modulus to change with Temperature.

% YMfun =@(location,state) 391.448 - 1.3160e-2*(state.u - 273.15) - 1.4838e5*(state.u - 273.15).^2

structuralProperties(smodel,'YoungsModulus',370e9,'PoissonsRatio',0.28)

structuralBodyLoad(smodel,'Temperature',Tresults,'Timestep',1);

I have also tried apply the results of all Timesteps.

structuralBC(smodel,'Face',[1 3 4],'Constraint','fixed')

Sresults = solve(smodel);

pdeplot3D(smodel,'ColorMapData',Sresults.VonMisesStress)

Running this yeilds ZERO STRESS. I have also tried a similar Axis-Symmetric model with the same results of ZERO STRESS.

The only time I get any STRESS results is when I apply a Pressure load:

% pLoad =@(location,state) pft_press(1).*location.z.^5 + pft_press(2).*location.z.^4 + pft_press(3).*location.z.^3 + pft_press(4).*location.z.^2 + pft_press(5).*location.z + pft_press(6) + location.x-location.x;

% structuralBoundaryLoad(smodel,'Face',3,'Pressure',pLoad)

Does anyone have any idea why I am not getting any Thermal STRESS?

### Accepted Answer

Ravi Kumar
on 2 Aug 2021

Got it. You don't have CTE specified along with the material properties. Refer to structuralProperties() call in the bi-metal example.

I also think you have unit mismatch. Your comments say geometry is in mm, but your properties are in SI (meters), check it. In addition, the geometry is pretty much like a line, I am unsure solving this as a solid model is the right approach. I got non-zero thermal stress using constant properties. I would suggest you verify with a linear model if the results are acceptable given the aspect ratio of the geometry.

##### 3 Comments

Ravi Kumar
on 5 Aug 2021

### More Answers (0)

### See Also

### Community Treasure Hunt

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

Start Hunting!