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

1 view (last 30 days)
Nicholas Morris on 2 Aug 2021
Commented: Ravi Kumar on 5 Aug 2021
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)
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;
Does anyone have any idea why I am not getting any Thermal STRESS?
##### 2 CommentsShowHide 1 older comment
Nicholas Morris on 2 Aug 2021
Hi Ravi,
Thanks for the quick reply. I have
tmodel.Geometry = gm
in my code. I must have missed it in my copy and pasting. I have also used
smodel.ReferenceTemperature = Tref
as well with no luck. However, I admit I don't totally understand what the Tref value is supposed to be. I tried 273.15 % [K], because my temperature profile is in Kelvin. I can try altering this value.
Also, you can try shorterning the length to 350 to speed up the code. Atleast until results are non-zero.
Thanks agian,
Nick

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.
Ravi Kumar on 5 Aug 2021
Axisymmetric might be a slightly better option than solid model, even then it is very long compared to the thickness. Would you be able to do some unit length analysis and use it to judge the full length?

R2021a

### Community Treasure Hunt

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

Start Hunting!

Translated by