Problem with transient heat transfer

[mabr]

Hello,

I am trying to solve a transient heat equation with space dependent heat flux. The case setup is shown in figure below. The top surface has changing heat flux as boundary condition. Hole surface has heat flux set to 0.1, and side/bottom surface are kept fixed at 110 °C.







I am interested in temperature distribution after 1 second. After specifying the time steps (10) and size (0.1) and running the case, I am getting temperature distribution that does not seem right. Even if I increase number of steps (100) and size (0.01), the result remains the same. Below are temperature distributions for both timesteps.







Next thing I remeshed the geometry and observed both result for different mesh densities. There is difference in temperature but the result still looks incorrect. The temperature distribution for both meshes is seen below







If I run the case up until 800 seconds (steps=200, size=4) the temperatures start to approach towards the steady state distribution. These two look correct. So it looks like the error is only in the early seconds and I am not sure how to tackle this.







Do you maybe have any advices on how to improve this study. Is there some explanation for the strange temperature distribution and how to solve it? I feel I have some errors in my transient parameters or am using the wrong approach here. The simulation and solver blocks are shown below.

Code: Select all

Simulation
  Max Output Level = 5
  Coordinate System = Cartesian
  Coordinate Mapping(3) = 1 2 3
  Simulation Type = Transient !steady state
  timestep intervals(1) = 10
  timestep sizes(1) = 0.1
  Output Intervals(1) = 1
  Steady State Max Iterations = 1
  Timestepping Method = BDF
  BDF Order = 3
  Solver Input File = new_sif.sif
  Post File = coarse_transient10_bicgstabl.vtu
  !Use Mesh Names = Logical True
  !element="p:3"
End

Code: Select all

Solver 1
  equation = Heat Equation
  procedure = "HeatSolve" "HeatSolver"
  variable = Temperature
  exec solver = Always
  stabilize = True
  bubbles = False
  lumped mass matrix = True
  optimize bandwidth = True
  steady state convergence tolerance = 1.0e-5
  nonlinear system convergence tolerance = 1.0e-9
  nonlinear system max iterations = 200
  nonlinear system relaxation factor = 0.5
  nonlinear system newton after iterations = 3
  nonlinear system newton after tolerance = 1.0e-5
  linear system solver = Iterative
  linear system iterative method = BiCGStabl
  linear system max iterations = 500
  linear system convergence tolerance = 1.0e-10
  bicgstabl polynomial degree = 2
  linear system preconditioning = ILU3
  linear system ilut tolerance = 1.0e-3
  linear system abort not converged = False
  linear system residual output = 10
  linear system precondition recompute = 1
End

Thank you very much for your time,
mabr

[kevinarden]

Heat takes time to transfer. If you touch a cold block with a hot source for 1 second, very little will happen, unless it is a very large extremely hot source. So the solution at 1 second is trivial.

Did you specify the initial conditions of the block? If not the default is zero degrees. A cold frozen block.

The steady state and 800 second solutions look good because the heat has had time to transfer.

[raback]

Hi,

Yes, you cannot just hit the object in constant temperature with a flux at t=0. In real life it would take some time for a temperature gradient to settle. This effect depends on the formulation. In finite volumes you would conserve energy by construction and you would probably see the first layer of cells accomodating your energy. In finite elements you try to capture a continuous temperature profile on a coarse mesh so that it weakly fulfills the equations making the solution seem unphysical in this case.

Would it be possible to rather use "external temperature" and "heat transfer coefficient". This often is a more physical way to introduce a temperature flux.

-Peter

[mabr]

Dear Kevin and Peter,

Thank you very much for your reply, I appreciate the time .

I agree with you that 1 second of heating is not a lot for large temperature evolution. The problems that I am solving have a ramp up phase (around 0.01s) from 0 W/m^2 to such heat flux profile, which then lasts for a few seconds and then there is a ramp down phase where this heat flux is reduced back to 0, so I am interested in temperature distribution only after a few seconds. For these cases, the peak power deposition is quite large (in the case presented here around ~2 MW/m^2) and there should be some temperature change on surface after a few seconds and this is what I want to detect.

I have set the initial temperature to 110 °C.

Peter, for specifying "external temperature" and "htc", I assume you meant changing boundary condition of temperature=110 °C. I changed it, but unfortunately the problem remains.

For example below is evolution after 5 seconds (50 steps and size is 0.1s). I also tried to use higher order of elements (element="p:3"), but there seem to be no much difference. The difference between neighboring nodes after 5 seconds is still too high.



For reference, I am attaching the sif file.
coarse_case.sif


As you suggest the approach I am using is probably not the best. Do you maybe have some further ideas to better physically introduce the flux in FEM methods that I could try out? I read about flux corrected transport methods, but I don't know if they are suitable for FEM methods and how easy it would be to implement them. If not, do you then suggest to check out finite volumes for these types of problems?

thank you for your help,
mabr

[kevinarden]

What are the dimensions of the block and the hole size?

[mabr]

The dimensions of the block are

length=1m
width=0.4m
height=0.2m
hole_radius=0.08m

for reference, I attach the mesh.
case_mesh.tar.gz

Do you think this plays a significant role?

regards,

[kevinarden]

Wanted to check units. They are correct.

Problem appears that boundary condition 2 is defined twice, in this case Elmer takes the second occurrence, so you don't have the BC on the hole that is expected. Changing this apparent typo improves the results on the first two steps

step1.png

(538.81 KiB) Not downloaded yet

step2.png

(77.53 KiB) Not downloaded yet

[mika]

Do you maybe have some further ideas to better physically introduce the flux in FEM methods that I could try out?

If not, do you then suggest to check out finite volumes for these types of problems?

I believe the mixed FEM formulation presented in the section "Mixed Approximation of the Poisson equation" of Elmer Models Manual might come quite close to equations obtained by a finite volume approximation. In this mixed FE approximation the temperature is approximated as discontinuous, while the normal component of the flux is continuous. The solver is written to handle a model problem and material parameters are not directly given physical interpretations. Despite its probable merits, I hesitate to give a strong recommendation to try it as it is a relatively new addition to Elmer and has not subject to extensive testing yet.

I also tried to use higher order of elements (element="p:3"), but there seem to be no much difference

Although the following will not solve the actual problem, I also mention that if one seeks a better spatial resolution by applying p-elements, files for better visualization of p-FEM results can be generated if the keyword "Lagrange Element Degree" is used. For more information on this fresh option see the commit message

https://github.com/ElmerCSC/elmerfem/co ... 8197959604

-- Mika

[mabr]

Wanted to check units. They are correct.

Problem appears that boundary condition 2 is defined twice, in this case Elmer takes the second occurrence, so you don't have the BC on the hole that is expected. Changing this apparent typo improves the results on the first two steps
step1.png
step2.png

Yes, you are correct, good catch. A rookie error from my side, thank you. I am using script that reads data from another code and maps it to Elmer case, so there is an error inside.

I fixed the error:

Code: Select all

Boundary Condition 2
Target Boundaries(1) = 1
  !name = 'plasma'
  heat flux bc = true
  heat flux = 0.1
End

Boundary Condition 3
Target Boundaries(1) = 2
  name = 'bc2'
  heat flux = 0.1
End

I am a bit surprised by the result then. Isn't default BC in Elmer specified as heat flux = 0? So if BC is not determined, Elmer will automatically set it to heat flux=0. Maybe I am wrong on this one.

After I change this, I still get the same result. I compared the same two nodes from both cases 1:1 and I get the same values. So maybe I am missing something here. Also your color map is a bit different than mine, so hard to say.




Do you mind posting your case, because I somehow cannot reproduce your results. Thank you.

regards,

[kevinarden]

This is the sif I used with your mesh. I edited the color map in paraview to improve the look of the temperature field since the values of temperature at each node are actually close in value to each other.

If you have BC 2 defined twice Elmer uses the second one.

coarse_case.sif

(92.83 KiB) Downloaded 115 times