[laserConvBC] Boundary condition with convection, heat source and radiation

Tobi · February 15, 2017, 07:16

Dear everybody,

last year I build a new BC named laserConvectionBC (https://bitbucket.org/shor-ty/laserconvectionbc). Here we can model a heat source (LASER) with additional convection. However, today I wanted to implement the radiation into the boundary condition and focused a lot of problems. Okay starting from the scratch. At the beginning we make an equilibrium of all the heat fluxes:

Thermal conductivity from Cell center to Face center
Convection
Heat source
Radiation

In addition $\Delta = 1/\delta$ while the capital delta is the one used in FOAM.

To identify all terms we get:

$q_1 = \lambda \nabla T = \lambda (T_c - T_f) / \delta$
$q_2 = \alpha (T_f - T_\infty)$

$q_4 = C (T_f^4 - T_\infty^4)$

Here the pre-factor in the last term includes the emission factor + Stefan Bolzmann Constant. However, now we equilibrate all fluxes to get:

which lead to:

$\lambda \nabla T = \lambda (T_c - T_f) / \delta +\alpha (T_f - T_\infty) -Pf +C (T_f^4 - T_\infty^4)$

For a own mixed BC, we have to get the following equation:
$T_f = \omega T_f^0 + (1-\omega)\left[T_c + \frac{\nabla^0 T}{\Delta}\right]$

The quantities are identified as:

given value by the user (Dirichlet)
$\frac{\nabla^0 T}{\Delta} :=$ Gradient expression (Neumann)
$\omega :=$ fraction which blends between Dirichlet and Neumann

The problem that I am focusing now is that we cannot get the above equation out of the combination of the fluxes. For example, the pre-factor of

and the gradient expression, which correspond to $(1-\omega)$ can be evaluated - after linearization - as:

$(1-\omega) = \frac{\lambda}{\lambda + \alpha\delta + \delta C T_f^{3*}}$

While the $\omega$ infront of $T_\infty$ is equal to:

$\omega = \frac{\alpha\delta + \delta C T_\infty^3}{\lambda + \alpha\delta + \delta C T_f^{3*}}$

The quantity $T_f^{3*}$ comes from the linearization and denotes the T field of the last iteration. It is obvious that $\omega \neq (1-\omega)$ . So the question is, how can I implement that in FOAM. I am out of ideas but my colleague will teach me (I hope so). However, I think it might be also interesting for other people. In principle I think I am not anymore able to use the mixed BC and have to make my own BC with the correction of everything ... ?

Each feedback is welcomed.

gtarang · May 31, 2017, 15:10

I have written this condition for my own CFD code (not OpenFOAM) and I did it by numerically solving the first equation and providing the value of Tf. Although I have started using OpenFOAM very recently, I will try to write it in OpenFOAM and share it here.

Thanks for the BC
-
Tarang

Tobi · May 31, 2017, 15:21

The easiest way is to put the radiation into the HTC and done.

kera · October 2, 2017, 04:08

Hallo Tobi and FOAMers,

Somehow even I am in need to use the convection radiation boundary condition, but unlike you I would like to use it with groovyBC. So I tried to simplify the radiation term and incorporated in the equations as follows:

Radiation:

$q_{rad} = \epsilon \sigma (T_f^4 - T_\infty^4) = \epsilon \sigma (T_f^2 + T_\infty^2)(T_f^2 - T_\infty^2)$

$q_{rad} = h_{rad} (T_f - T_\infty)$

where $h_{rad} = \epsilon \sigma (T_f^2 + T_\infty^2)(T_f + T_\infty)$

Convection:

$q_{conv} = h_{conv} (T_f - T_\infty)$

Now if I make the energy balance for all 3 modes of heat transfer I have:

$q_{cond} + q_{conv} + q_{rad} = 0$

$k \frac{T_f - T_c}{\delta} + h_{conv} (T_f - T_\infty) + h_{rad} (T_f - T_\infty) = 0$

After isolating

and dividing by $(h_{conv}+h_{rad})\delta$ I get:

$\frac{k}{(h_{conv}+h_{rad}) \delta} T_c + T_\infty = T_f \Big[1+ \frac{k}{(h_{conv}+h_{rad}) \delta}\Big]$

on further simplification I get

$T_f = f T_\infty + (1 - f) T_c$

where f is $\Big[1+ \frac{k}{(h_{conv}+h_{rad}) \delta}\Big]^{-1}$

Now implementing it in OpenFOAM with groovyBC

Code:

    myPatch
    {
        type            groovyBC;
        refValue        uniform 933.15;
        refGradient     uniform 0;
        valueFraction   uniform 1;
        value           uniform 933.15;
        valueExpression "T_inf";
        gradientExpression "0";
        fractionExpression "1.0/(1.0+k/(mag(delta())*(h_conv +h_rad)))";
        variables       
6
(
"k= x ;"
"sig_rad=5.6e-8;"
"emmi=0.7;"
"T_inf=933.15;"
"h_conv =  y ;
"h_rad=sig_rad*emmi*(pow(T,2)+pow(T_inf,2))*(T+T_inf);"
)
;
}

Questions:

1) I am not sure if this is the right way of doing it ---> 2nd Opinion is always useful

2) Are there any reference case with an analytical solution for this kind of problems?
3) Were you (Tobi) able to implement your boundary condition for any of your simulations with sane results?

Cheers!
Ricky

PS: The Idea was taken from here.

Tobi · October 2, 2017, 04:50

Hi,

actually the post is not related to the boundary condition itself (: ... I could not find any tutorial about modeling the radiation term in the way I did or you do. However, I know it is a common method to put the radiation term into the HTC because of the high non-linearity. Your equations are correct and the modeling and implementation are correct too. Keep in mind that based on the explicit term T_f in the HTC, you should make a few loops to converge everything. That was actually the reason why I put the radiation term directly into a temperature depended HTC (by the way, nice that you hide your k and conv part using x and y :P ). Based on the fact that you derived everything yourself, I think you know what I am talking about.

Again you have

as explicit term in the boundary condition which makes your problem stiff.

kera · October 2, 2017, 05:03

Hallo Tobi,

The reason for posting it here, was because of the high non-linearity term and therefore I was not sure if at all my approach is in the right direction. --> Besides I was assuming that one need very small time steps (w.r.t. the grid) to see the effects, which could have been lot more easier if there were any analytical solution for such problems or perhaps a reference case.

Quote:

Originally Posted by Tobi

you should make a few loops to converge everything.

If I understand it correctly did you mean using loops at the boundary condition? Is it possible to do that using groovyBC? if yes how? I would really appreciate your help.

Regards,
Ricky

Tobi · October 2, 2017, 06:08

Hi Ricky,

actually, I made a FOAM case where I was calculating the HTC for convection. Thus, I got an temperature depended heat flux and corresponding HTC. The same I did analytically for the radiation. Actually, I added the convective and radiation heat flux and build the temperature depended HTC which I was fitting with an polynomial which is simple to use as BC.

If you want to have the real problem, then you have to make outer-corrections. I have no idea how much you need because of the high nonlinearity of the explicit

term. A boundary loop cannot be used and is not possible in FOAM (as far as I know). Just as an hint: The face values depend on the internal field too (gradients). Thus, you have to update everything; the internalField and boundaryField. The time step itself is questionable. I agree that larger time-steps will make the problem harder to solve. I can imagine that you need under-relaxation in combination of the outer-loops (PIMPLE) and residual control.

kera · October 2, 2017, 09:09

Hallo Tobi,

Thank you for your suggestions, and I will look in to it, indeed the simulation settings need to be more stiffer for such problems. I will let you know if I come up with something else or get stuck in between

Quote:

Originally Posted by Tobi

actually, I made a FOAM case where I was calculating the HTC for convection. Thus, I got an temperature depended heat flux and corresponding HTC. The same I did analytically for the radiation. Actually, I added the convective and radiation heat flux and build the temperature depended HTC which I was fitting with an polynomial which is simple to use as BC.

I also think that your approach is also interesting, I would give it a try as well.

Regards,
Ricky

Tobi · October 2, 2017, 09:25

I used that approach in my Ph.D thesis. So what I did actually was:

cht solver case: cooling an aluminum alloy from solidus to room temperature via free convection
Calculating the average wall temperature and the corresponding heat flux
Using the average wall temperature for calculating the radiation heat flux
Add both heat flux
Calculate the HTC for the given average wall temperature
Fitting the temperature depended HTC using a polynomial
Using the polynomial to calculate the HTC which includes convection and radiation and thus, the radiation term can be removed from the energy balance

That was my approach.
Good luck.

Tobi

wolfindark · August 3, 2020, 07:32

Tobias, do you have an updated version of your laserBC for the new versions of openfoam?

gtarang · August 3, 2020, 13:43

Hi,
I used inspiration from Tobi's boundary condition and it works in OpenFOAM-6 and it should work in v7 (I suspect no major changes in of7) but it will not work in v8. Well, it is not documented as exhaustive as Tobi's, but basic documentation is there in the header file. If you need help, pls reach out to me.

https://github.com/gargtarang/OpenFO...bulentHeatFlux

Tobi · August 3, 2020, 14:07

Hi, up to now I did not port it to OpenFOAM 8 as I was not aware that people are using the boundary condition, as it is very limited (only in one axis).

Furthermore, I did not programmed it in a nice way

I might update it somewhen but first, there are other things to update first. I do have too much material that wants to be updated

wolfindark · August 3, 2020, 21:10

Quote:

Originally Posted by Tobi

Hi, up to now I did not port it to OpenFOAM 8 as I was not aware that people are using the boundary condition, as it is very limited (only in one axis).

Furthermore, I did not programmed it in a nice way

I might update it somewhen but first, there are other things to update first. I do have too much material that wants to be updated

Thanks Tobias,
I am trying to use it icoReactingMultiphaseInterFoam (OF v2006). Following compilation error comes out:

Code:

In file included from laserConvectionFvPatchField.H:450,
                 from laserConvectionFvPatchFields.H:44,
                 from laserConvectionFvPatchFields.C:43:
laserConvectionFvPatchField.C: In constructor ‘Foam::laserConvectionFvPatchField<Type>::laserConvectionFvPatchField(const Foam::fvPatch&, const Foam::DimensionedField<Type, Foam::volMesh>&, const Foam::dictionary&)’:
laserConvectionFvPatchField.C:346:66: error: no matching function for call to ‘Foam::List<Foam::Vector<double> >::append(Foam::ITstream&)’
  346 |             sourceCenters_.append(sourceDict.lookup("spotCenter"));
      |                                                                  ^

hoai nguyen · August 19, 2020, 04:12

Hi, I am using OpenFOAM to simulate laser bonding of two materials (glasses). I have already compiled laserconvectionBC of Dr. Tobi. Nevertheless, this code is applyed for one material. I do not know which solvers I enable use my case. Could you please give me advise?
Thank you so much!

ramonkey · September 2, 2020, 05:04

Hallo Everyone,

i have a problem with decomposing with the laser BC.
When i decompose and let the Simulation solve with the following command:
"mpirun -np 8 icoReactingMultiphaseInterFoam -parallel"
it gives me the following output:
"For patch Top only convection is applied, no laser source found."
And it don't uses the laser source.
also the source disappears by the decomposing process in the T-Files in the Processor Folders.

When i let the Simulation run without decomposing the source is used and there are no Problems.

The Version of OpenFoam is v1912.

If somebody could help me, that would be very nice!

Best regards!

Tobi · September 2, 2020, 15:26

Interesting point. I had no problems decomposing the guy during my PhD. however, things might change

can you provide a "simple" test case so that I can check it?

Tobi · September 2, 2020, 15:27

Quote:

Originally Posted by gtarang

Hi,
I used inspiration from Tobi's boundary condition and it works in OpenFOAM-6 and it should work in v7 (I suspect no major changes in of7) but it will not work in v8. Well, it is not documented as exhaustive as Tobi's, but basic documentation is there in the header file. If you need help, pls reach out to me.

https://github.com/gargtarang/OpenFO...bulentHeatFlux

Just for me , ... what changes do you have in the condition?

ramonkey · September 3, 2020, 04:31

Quote:

Originally Posted by Tobi

Interesting point. I had no problems decomposing the guy during my PhD. however, things might change

can you provide a "simple" test case so that I can check it?

https://github.com/ramonwagner1997/SLM-Simulation
thats the whole case yet. I think it`s not that complex.
i build it up on the solidMelting2D Case.

And Thanks for the quick response!

ramonkey · September 22, 2020, 05:59

Quote:

Originally Posted by Tobi

Interesting point. I had no problems decomposing the guy during my PhD. however, things might change

can you provide a "simple" test case so that I can check it?

Hey Tobi,
just wanted to ask if you had the Time to look at the Case i sent.

gtarang · September 22, 2020, 10:30

Quote:

Originally Posted by Tobi

Just for me , ... what changes do you have in the condition?

Sorry for replying late
1. 3 types of sources -> circular, ellipse, linear
2. Evaporation heat loss
3. Radiative heat loss
4. Sinusoidal movement of heat spot

Radiative and evaporative heat loss is not linearised which I will solve using Ridder's method later.

-
gtarang

February 15, 2017, 07:16	Boundary condition with convection, heat source and radiation	#1
Tobi Super Moderator Tobias Holzmann Join Date: Oct 2010 Location: Tussenhausen Posts: 2,708 Blog Entries: 6 Rep Power: 51	Dear everybody, last year I build a new BC named laserConvectionBC (https://bitbucket.org/shor-ty/laserconvectionbc). Here we can model a heat source (LASER) with additional convection. However, today I wanted to implement the radiation into the boundary condition and focused a lot of problems. Okay starting from the scratch. At the beginning we make an equilibrium of all the heat fluxes: Thermal conductivity from Cell center to Face center Convection Heat source Radiation In addition $\Delta = 1/\delta$ while the capital delta is the one used in FOAM. To identify all terms we get: $q_1 = \lambda \nabla T = \lambda (T_c - T_f) / \delta$ $q_2 = \alpha (T_f - T_\infty)$ $q_4 = C (T_f^4 - T_\infty^4)$ Here the pre-factor in the last term includes the emission factor + Stefan Bolzmann Constant. However, now we equilibrate all fluxes to get: which lead to: $\lambda \nabla T = \lambda (T_c - T_f) / \delta +\alpha (T_f - T_\infty) -Pf +C (T_f^4 - T_\infty^4)$ For a own mixed BC, we have to get the following equation: $T_f = \omega T_f^0 + (1-\omega)\left[T_c + \frac{\nabla^0 T}{\Delta}\right]$ The quantities are identified as: given value by the user (Dirichlet) $\frac{\nabla^0 T}{\Delta} :=$ Gradient expression (Neumann) $\omega :=$ fraction which blends between Dirichlet and Neumann The problem that I am focusing now is that we cannot get the above equation out of the combination of the fluxes. For example, the pre-factor of and the gradient expression, which correspond to $(1-\omega)$ can be evaluated - after linearization - as: $(1-\omega) = \frac{\lambda}{\lambda + \alpha\delta + \delta C T_f^{3}}$ While the $\omega$ infront of $T_\infty$ is equal to: $\omega = \frac{\alpha\delta + \delta C T_\infty^3}{\lambda + \alpha\delta + \delta C T_f^{3}}$ The quantity $T_f^{3}$ comes from the linearization and denotes the T field of the last iteration. It is obvious that $\omega \neq (1-\omega)$ . So the question is, how can I implement that in FOAM. I am out of ideas but my colleague will teach me (I hope so). However, I think it might be also interesting for other people. In principle I think I am not anymore able to use the mixed BC and have to make my own BC with the correction of everything ... ? Each feedback is welcomed. dats and hoai nguyen like this. __________________ Keep foaming, Tobias Holzmann Last edited by Tobi; September 19, 2017 at 04:21.*

May 31, 2017, 15:21		#3
Tobi Super Moderator Tobias Holzmann Join Date: Oct 2010 Location: Tussenhausen Posts: 2,708 Blog Entries: 6 Rep Power: 51	The easiest way is to put the radiation into the HTC and done. dats likes this. __________________ Keep foaming, Tobias Holzmann

October 2, 2017, 04:08		#4
kera Member Ricky Join Date: Jul 2014 Location: Germany Posts: 78 Rep Power: 11	Hallo Tobi and FOAMers, Somehow even I am in need to use the convection radiation boundary condition, but unlike you I would like to use it with groovyBC. So I tried to simplify the radiation term and incorporated in the equations as follows: Radiation: $q_{rad} = \epsilon \sigma (T_f^4 - T_\infty^4) = \epsilon \sigma (T_f^2 + T_\infty^2)(T_f^2 - T_\infty^2)$ $q_{rad} = h_{rad} (T_f - T_\infty)$ where $h_{rad} = \epsilon \sigma (T_f^2 + T_\infty^2)(T_f + T_\infty)$ Convection: $q_{conv} = h_{conv} (T_f - T_\infty)$ Now if I make the energy balance for all 3 modes of heat transfer I have: $q_{cond} + q_{conv} + q_{rad} = 0$ $k \frac{T_f - T_c}{\delta} + h_{conv} (T_f - T_\infty) + h_{rad} (T_f - T_\infty) = 0$ After isolating and dividing by $(h_{conv}+h_{rad})\delta$ I get: $\frac{k}{(h_{conv}+h_{rad}) \delta} T_c + T_\infty = T_f \Big[1+ \frac{k}{(h_{conv}+h_{rad}) \delta}\Big]$ on further simplification I get $T_f = f T_\infty + (1 - f) T_c$ where f is $\Big[1+ \frac{k}{(h_{conv}+h_{rad}) \delta}\Big]^{-1}$ Now implementing it in OpenFOAM with groovyBC Code: myPatch { type groovyBC; refValue uniform 933.15; refGradient uniform 0; valueFraction uniform 1; value uniform 933.15; valueExpression "T_inf"; gradientExpression "0"; fractionExpression "1.0/(1.0+k/(mag(delta())(h_conv +h_rad)))"; variables 6 ( "k= x ;" "sig_rad=5.6e-8;" "emmi=0.7;" "T_inf=933.15;" "h_conv = y ; "h_rad=sig_rademmi(pow(T,2)+pow(T_inf,2))(T+T_inf);" ) ; } *Questions:* 1) I am not sure if this is the right way of doing it ---> 2nd Opinion is always useful 2) Are there any reference case with an analytical solution for this kind of problems? 3) Were you (Tobi) able to implement your boundary condition for any of your simulations with sane results? Cheers! Ricky PS: The Idea was taken from here. mbookin likes this. __________________ If it is easy, then something is fishy! Last edited by kera; October 2, 2017 at 10:25.

October 2, 2017, 04:50		#5
Tobi Super Moderator Tobias Holzmann Join Date: Oct 2010 Location: Tussenhausen Posts: 2,708 Blog Entries: 6 Rep Power: 51	Hi, actually the post is not related to the boundary condition itself (: ... I could not find any tutorial about modeling the radiation term in the way I did or you do. However, I know it is a common method to put the radiation term into the HTC because of the high non-linearity. Your equations are correct and the modeling and implementation are correct too. Keep in mind that based on the explicit term T_f in the HTC, you should make a few loops to converge everything. That was actually the reason why I put the radiation term directly into a temperature depended HTC (by the way, nice that you hide your k and conv part using x and y :P ). Based on the fact that you derived everything yourself, I think you know what I am talking about. Again you have as explicit term in the boundary condition which makes your problem stiff. __________________ Keep foaming, Tobias Holzmann

October 2, 2017, 06:08		#7
Tobi Super Moderator Tobias Holzmann Join Date: Oct 2010 Location: Tussenhausen Posts: 2,708 Blog Entries: 6 Rep Power: 51	Hi Ricky, actually, I made a FOAM case where I was calculating the HTC for convection. Thus, I got an temperature depended heat flux and corresponding HTC. The same I did analytically for the radiation. Actually, I added the convective and radiation heat flux and build the temperature depended HTC which I was fitting with an polynomial which is simple to use as BC. If you want to have the real problem, then you have to make outer-corrections. I have no idea how much you need because of the high nonlinearity of the explicit term. A boundary loop cannot be used and is not possible in FOAM (as far as I know). Just as an hint: The face values depend on the internal field too (gradients). Thus, you have to update everything; the internalField and boundaryField. The time step itself is questionable. I agree that larger time-steps will make the problem harder to solve. I can imagine that you need under-relaxation in combination of the outer-loops (PIMPLE) and residual control. __________________ Keep foaming, Tobias Holzmann

May 31, 2017, 15:10		#2
gtarang Member Tarang Join Date: Feb 2011 Location: Delhi, India Posts: 47 Rep Power: 15	I have written this condition for my own CFD code (not OpenFOAM) and I did it by numerically solving the first equation and providing the value of Tf. Although I have started using OpenFOAM very recently, I will try to write it in OpenFOAM and share it here. Thanks for the BC - Tarang

October 2, 2017, 09:25		#9
Tobi Super Moderator Tobias Holzmann Join Date: Oct 2010 Location: Tussenhausen Posts: 2,708 Blog Entries: 6 Rep Power: 51	I used that approach in my Ph.D thesis. So what I did actually was: cht solver case: cooling an aluminum alloy from solidus to room temperature via free convection Calculating the average wall temperature and the corresponding heat flux Using the average wall temperature for calculating the radiation heat flux Add both heat flux Calculate the HTC for the given average wall temperature Fitting the temperature depended HTC using a polynomial Using the polynomial to calculate the HTC which includes convection and radiation and thus, the radiation term can be removed from the energy balance That was my approach. Good luck. Tobi arvindpj and kera like this. __________________ Keep foaming, Tobias Holzmann

August 3, 2020, 07:32		#10
wolfindark Member Neilson Whit Join Date: Aug 2011 Posts: 74 Rep Power: 14	Tobias, do you have an updated version of your laserBC for the new versions of openfoam?

August 3, 2020, 13:43		#11
gtarang Member Tarang Join Date: Feb 2011 Location: Delhi, India Posts: 47 Rep Power: 15	Hi, I used inspiration from Tobi's boundary condition and it works in OpenFOAM-6 and it should work in v7 (I suspect no major changes in of7) but it will not work in v8. Well, it is not documented as exhaustive as Tobi's, but basic documentation is there in the header file. If you need help, pls reach out to me. https://github.com/gargtarang/OpenFO...bulentHeatFlux Tobi and wolfindark like this.

August 3, 2020, 14:07		#12
Tobi Super Moderator Tobias Holzmann Join Date: Oct 2010 Location: Tussenhausen Posts: 2,708 Blog Entries: 6 Rep Power: 51	Hi, up to now I did not port it to OpenFOAM 8 as I was not aware that people are using the boundary condition, as it is very limited (only in one axis). Furthermore, I did not programmed it in a nice way I might update it somewhen but first, there are other things to update first. I do have too much material that wants to be updated __________________ Keep foaming, Tobias Holzmann

August 19, 2020, 04:12	Laser bonding two materials	#14
hoai nguyen New Member Đà Nẵng Join Date: Jun 2019 Posts: 4 Rep Power: 6	Hi, I am using OpenFOAM to simulate laser bonding of two materials (glasses). I have already compiled laserconvectionBC of Dr. Tobi. Nevertheless, this code is applyed for one material. I do not know which solvers I enable use my case. Could you please give me advise? Thank you so much!

September 2, 2020, 05:04	Problems with Decomposing	#15
ramonkey New Member Ramon Join Date: May 2020 Posts: 3 Rep Power: 5	Hallo Everyone, i have a problem with decomposing with the laser BC. When i decompose and let the Simulation solve with the following command: "mpirun -np 8 icoReactingMultiphaseInterFoam -parallel" it gives me the following output: "For patch Top only convection is applied, no laser source found." And it don't uses the laser source. also the source disappears by the decomposing process in the T-Files in the Processor Folders. When i let the Simulation run without decomposing the source is used and there are no Problems. The Version of OpenFoam is v1912. If somebody could help me, that would be very nice! Best regards!

September 2, 2020, 15:26		#16
Tobi Super Moderator Tobias Holzmann Join Date: Oct 2010 Location: Tussenhausen Posts: 2,708 Blog Entries: 6 Rep Power: 51	Interesting point. I had no problems decomposing the guy during my PhD. however, things might change can you provide a "simple" test case so that I can check it? __________________ Keep foaming, Tobias Holzmann

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Radiation in semi-transparent media with surface-to-surface model?	mpeppels	CFX	11	August 22, 2019 07:30
Question about adaptive timestepping	Guille1811	CFX	25	November 12, 2017 17:38
Monte Carlo Simulation: H-Energy is not convergating & high Incident Radiation	volleyHC	CFX	5	April 3, 2016 05:41
Radiation heat transfer boundary condition	natrask	OpenFOAM Programming & Development	0	February 8, 2015 09:05
GETVAR Error in Multiband Monte Carlo Radiation Simulation with Directional Source	silvan	CFX	3	June 16, 2014 09:49