internalHeatSource for cylindrical geometry

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,

Alexander

Answer 1

[Ravi Kumar]

Having an axisymmetric analsysis, which is not supported yet, would simplify your setup.Here is an alternative approach which might solve your problem:

Replace:

internalHeatSource(thermalmodel,f_gas,'Face',2); 
internalHeatSource(thermalmodel,f_solid,'Face',1);
internalHeatSource(thermalmodel,f_solid,'Face',3);

with a single call:

hsFcn = @(region,state) assignHS(region,state,f_gas, f_solid);
internalHeatSource(thermalmodel,hsFcn)

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);
else
    q = f_solid(region,state);
end
end

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.

Regards,

Ravi

[Alexander]

Dear Ravi,

thank you for the advise. I did implement your code but I still see a diffeence between the fem and exact solutions (shich is now much smaller). It cannot be because of unaccuracy of numeric solution, so there is still something wrong with the heat sources.

Please, let me know if you find a different solution of this problem.

You mentioned some unsupported axisymmetric analsysis. Is that a separate tool?

thank you

Alexander

[Ravi Kumar]

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.

Regards,

Ravi

Answer 2

[Alexander]

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

[Alexander]

Hallo Ravi,

I do not know if that can help you, but I did another experiment:

I wrote a code using procedures from the Matlab2016b where you again have three models:

• Axi-symmetrical
• 2D
• analytical

If I run the code with Matlab2016b I get agreement between 2D and analytical, but the axi-symmetrical solution is different.

Surprisingly, I can run the same code with Matlab2019b and I get all three in excellent agreement.

That can be seen as a solution for time being but I do not have much confidence in this story.

Best regards,

Alexander

[Ravi Kumar]

Hi Alexander,

Thanks for the additional validation information. I will go through it and get back to you.

Regards,

Ravi

[Alexander]

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.

Alexander