internalHeatSource for cylindrical geometry

1 view (last 30 days)
Alexander on 5 Mar 2020
Commented: Alexander on 10 Mar 2020
Dear all,
I have the following problem using the Matlab PDE toolbox for thermal analysis:
If I call the Matlab function internalHeatSource for different segments it looks like as if the latest call overwrites the previous calls.
Problem: I have a solid rod inside a solid cylindrical shell separated with a layer of gas. All three segments have body heating proportional to their density (gamma heating). The shell is cooled with water from outside. The geometry has a cylindrical symmetry. I’m looking for a steady state solution.
The boundary conditions at the axis of symmetry (y=0), top and bottom (x=x1,x2) have thermal flux zero. I provide the Matalab code. At the end of the code the exact analytical solution of the problem is calculated and compared to the fem results.
Considering that the segments have different heatings I have to call internalHeatSourceat least two times. It looks like that the latest call of the Matlab function internalHeatSource overwrites the previous calls.
Case 1:
rho_solid = 1.e3; % kg/m3
rho_gas = 1.; % kg/m3
The fem calculates higher central temperature as if in the calculations it uses rho_gas = 1.e3.
Case 2:
Use the same heating value for solid and gas segments
rho_solid = 1.e3; % kg/m3
rho_gas = 1.e3; % kg/m3
Both solutions, fem and analytical are identical.
Case 3:
Use the correct values for the densities
rho_solid = 1.e3; % kg/m3
rho_gas = 1.; % kg/m3
but swap the sequence of the internalHeatSource calls: Let the call for segment 2 (with rho_gas = 1.;) be the last in the code.
The central temperature produced by fem is incredibly low as if for all segments the gas density is used.
Could you please help me with finding a solution for this strange behavior.
Thank you in advance,

Answers (2)

Ravi Kumar
Ravi Kumar on 5 Mar 2020
Having an axisymmetric analsysis, which is not supported yet, would simplify your setup.Here is an alternative approach which might solve your problem:
with a single call:
hsFcn = @(region,state) assignHS(region,state,f_gas, f_solid);
where you define assignHS to assigns different heat source values based on the subdmains as:
function q = assignHS(region,state,f_gas, f_solid)
if region.subdomain == 2
q = f_gas(region,state);
q = f_solid(region,state);
You can past assignHS function at the end of your script or as a seperate file. The way you have the setup should also work, it might be missing to detect the variation in value of heat source based on the query. I will investigate this further. Sorry for the inconvenince this has caused you.
Ravi Kumar
Ravi Kumar on 7 Mar 2020
Hi Alexander,
I did notice the difference. The mesh is highly resolved, so I would agree with you the differece is large to be attributed for numerical error. When you say exact solution, do you have a reference that I can read and understand the code that you have used to compte the exact solution at the end of the code you shared.

Sign in to comment.

Alexander on 9 Mar 2020
Dear Ravi,
I did a small study which you can see in the attached matlab file.
I did the following calculations:
  1. Axi-symmetrical model with a function for the heat source as you proposed.
  2. A 2D model
  3. A 2D model but calculated with 2016b version. This version I use routinely.
  4. Analytic solution (explanation of the equations you can find in pfd file)
What I see is that only the 2D model from 2016b produces the same result as the analytic solution. Both fem calculations from 2019b produce different result. I have intention to switch to the 2019 release but as you see for the moment I have reasons to stay with the old version.
I will appreciate a lot if you could find a solution.
Thank you,
Alexander on 10 Mar 2020
Dear Ravi,
I think I know why the fem gives a wrong answer. The code you sent produces the same mistake because if the assignHS function is called for multiple nodes from different segments, only the value for the last segment is assigned for the whole array. I adjusted the code now for both assignHS and assignK, and now all solutions are in excellent agreement.( I checked that only for 2016b for the moment.)
Thank you for your assistance and I hope that issue will be resolved in the next Matlab revision.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!