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

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?

[Ravi Kumar]

Hi Nicholas,

I tried to run the script. In the thermal portion, I added the missing line tmodel.Geometry = gm. But the thermal simulation is taking very long time. Were you able to get the expected thermal results? Do the temperature contours look correct?

You also might want to assign reference temperature in structural model, using smodel.ReferenceTemperature = Tref. Use a Tref that is consistent with the thermal simulation.

Regards,

Ravi

[Nicholas Morris]

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

Answer 1

[Ravi Kumar]

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.

[Nicholas Morris]

Geometry: [mm]

Specific heat: [J/gK]

Thermal Conductivity: [W/mmK]

Mass Density: [g/mm^3]

Youngs Modulus: [Pa]

Poisons Ratio: [unitless]

Coef. of Thermal Exp. : [1/K]

Is this correct?

[Nicholas Morris]

OK, great! I got non-zero thermal stress :)

Ravi, I would love to follow up on "the geometry is pretty much like a line, I am unsure solving this as a solid model is the right approach". I was thinkining about us the axis-symmetric approach. What would you suggest?

Also, what do you mean by "I would suggest you verify with a linear model if the results are acceptable given the aspect ratio of the geometry. "

Thanks!

[Ravi Kumar]

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?