Radiant energy exchange between neighboring surfaces

General
1

Dear All,

I was trying to figure out how to implement radiation between two surface regions.

In detail, I would like to calculate the view factors of each CVFID of the first region respect to each CVFID of the second region and then calculate the heat loss based on the radiosity calculated.

Skipping the particular implementation details, I’d like to figure out how to use GWFL to calculate relations between each single boundary element “k” and all the others “j”, obtaining a matrix F_kj to be used in the assembly procedure.

Thank you.
Lorenzo

2

Dear Lorenzo,

If I understand this correctly, you have to use the raytracing transformation in GWFL. There is a very detailed demo for contact here:
https://git.savannah.nongnu.org/cgit/getfem.git/tree/interface/tests/python/demo_large_sliding_contact.py
if you look at the code for direct_generic_assembly=True .

Basically you need to add a raytracing transformation to your model (you give it any name you like, let’s call it “radiation_trans”) and configure a source and target surface.

md.add_raytracing_transformation("radiation_trans", release_dist)
md.add_slave_contact_boundary_to_raytracing_transformation("radiation_trans", m, "u", SOURCE)
md.add_master_contact_boundary_to_raytracing_transformation("radiation_trans", m, "u", TARGET)

Then for any GWFL term you define on SOURCE region, at any integration point on SOURCE you can get the corresponding raytraced point on TARGET region with a syntax like:

md.add_nonlinear_term(mim, "(d-(X-Interpolate(X,radiation_trans)))*Test_d", SOURCE)

If d is a scalar field variable defined on SOURCE, this example term will evaluated the field d as the raytraced distance to TARGET.

Because there is no guarantee that the raytraced point always exist, GetFEM allows you to switch between different expressions based on whether the raytraced point exists or not, using the Interpolate_filter operator, e.g.:

md.add_nonlinear_term(mim,
  f"Interpolate_filter(radiation_trans,"
                     f"(d-(X-Interpolate(X,radiation_trans)))*Test_d,"
                     f"(d-{release_dist})*Test_d)", SOURCE)

This will evaluate d either as the raytraced distance, if raytracing succeed or as releast_dist if raytracing fails·

3

Hello Kostas,
thank you for your reply. Very interesting feature the raytracing transformation.

As far as I understood, the raytracing mathod maps each integration point on SOURCE (slave) to their projection on TARGET (master), and the projection (for raytracing) is perpendicular to the slave surface.

In this case, I’m not sure it is suitable for radiation purposes, because in my case I need to calculate the distance between each integration point on the SOURCE surface respect to each integration point on the TARGET surface. To calculate the view factors of each face on SOURCE respect to each face on TARGET i must calculate what shown in the following pictures

ViewFactorePicture.PNG
image.png

and then integrate the heat exchange as follows,
image.png

considering the calculated view factors F_k-j and also the reflections of radiation accounted for by 1-epsilon_j.

where epsilon_j is the emissivity of the surface material
Qj is the energy loss from each surface to the temperature of each surface
delta_kj is the unit tensor
sigma is the Stefan Boltzman constant
Aj the area of each element

Please, let me know if I understood correctly the raytracing features and if you have any other suggestion for my case.
Thank you,
Lorenzo

4

ok, now I understand your question better.

Then you should have a look at
https://git.savannah.nongnu.org/cgit/getfem.git/tree/interface/tests/python/check_secondary_domain.py

Basically you want to calculate nested integrals instead of simple integrals. This is done in GWFL with the Secondary_domain operator.

The sum equation you provide is not very helpful, you should use the corresponding integral that this sum originates from. GetFEM converts integrals to sums itself, you only work with integrals.

Basically you need to do something like:

mim_surface = .... # filtered mim on surface

md = gf.Model('real')
md.add_fem_variable('theta', mftheta)
md.add_fem_variable('Q', mfQ)
...
md.add_standard_secondary_domain('secdom', mim_surface)
md.add_macro("dist", "Norm(X-Secondary_Domain(X))") # this is "s" in your equations
md.add_macro("cos1", "-Normalized(X-Secondary_Domain(X)).Normal")
md.add_macro("cos2", "Normalized(X-Secondary_Domain(X)).Secondary_Domain(Normal)")
...
md.add_nonlinear_twodomain_term(mim, "...", SURFACE_RG, "secdom")

Basically if X is (an integration) point k, then Secondary_Domain(X) are all (integration) points j defined in mim_surface.

If you provide the Equations in pre-descritized form I could help you more.

5

Hello Kostas,
that solution looks very good! It seems to be the right method for my purpose.
Tomorrow I’ll try it out and let you know the outcomes.
Thank you thousands!
Lorenzo

6

Hello Kostas.
Maybe I was a little bit optimistic, because it seems that the problem requires a triple integral for grey bodies (the ones with emissivity < 1 and accounting also for the relevant radiation reflections).
So, I’d like to try before an approximation of the problem without reflections for which the double integral should work.
For this purpose, I’d kindly ask your help to figure out the correct implementation in GetFEM.

Referring to the following picture with index “i” for the mim point in mesh_1 and “j” for the mim point on mesh_2

image.png
image.png572×604 25.8 KB

the amount of heat radiated from the infinitesimal surface dAi to surface dAj is
image.png

In this equation, J_i is the radiosity, that can be approximated as just the emitted heat, J_i = E_i = epsilon * sigma * T_i^4.

Similarly, the heat from “j” toward “i” is dq_{j->i}, that has the same formula as above but replacing J_i with J_j.

Integrating with the secondary_domain technique the difference of dq_{i->j} and dq_{j->i} we achieve the net heat flowing out of mesh_1 surface A_i
image.png

swapping the indexes we can achieve q_{ji}, the net heat flowing out of surface A_j on mesh_2.

Considering that I have to calculate the temperature change on both mesh_1 and mesh_2 mutually radiating, can you suggest the correct way to implement it starting from your previous suggestion?

Lots of thanks,
Lorenzo

7

Dear @lorenzo.iron

When I run this script:

#!/usr/bin/env python3
# -*- coding: UTF8 -*-

import getfem as gf
import numpy as np

gf.util_trace_level(1)
gf.util_warning_level(1)

# Mesh size
NX = 5
NY = 10

LX = 2.       # Length
LY = 4.       # Height

WX = 0.5      # Gap length between block 1 an block 2

fem_order = 2

# Mesh generation
m1 = gf.Mesh("import", "structured",
             f"GT='GT_QK(2,1)';ORG=[0,0];SIZES=[{LX},{LY}];NSUBDIV=[{NX},{NY}]")
m2 = gf.Mesh("import", "structured",
             f"GT='GT_QK(2,1)';ORG=[{LX+WX},{0.5*LY}];SIZES=[{LX},{LY}];NSUBDIV=[{NX},{NY}]")

# Mesh region definition
XM1_RG = 11
XM2_RG = 21
XP1_RG = 12
XP2_RG = 22
YM1_RG = 13
YM2_RG = 23
YP1_RG = 14
YP2_RG = 24

m1.set_region(XM1_RG, m1.outer_faces_in_box(     [    -1e-5,-1e-5],          [1e-5,LY+1e-5]))
m2.set_region(XM2_RG, m2.outer_faces_in_box([LX+WX-1e-5,LY/2-1e-5],[LX+WX+1e-5,3*LY/2+1e-5]))

m1.set_region(XP1_RG, m1.outer_faces_in_box(         [LX-1e-5,-1e-5],         [LX+1e-5,LY+1e-5]))
m2.set_region(XP2_RG, m2.outer_faces_in_box([2*LX+WX-1e-5,LY/2-1e-5],[2*LX+WX+1e-5,3*LY/2+1e-5]))

m1.set_region(YM1_RG, m1.outer_faces_in_box(         [-1e-5,-1e-5],          [LX+1e-5,1e-5]))
m2.set_region(YM2_RG, m2.outer_faces_in_box([LX+WX-1e-5,LY/2-1e-5],[2*LX+WX+1e-5,LY/2+1e-5]))

m1.set_region(YP1_RG, m1.outer_faces_in_box(         [-1e-5,LY-1e-5],         [LX+1e-5,LY+1e-5]))
m2.set_region(YP2_RG, m2.outer_faces_in_box([LX+WX-1e-5,3*LY/2-1e-5],[2*LX+WX+1e-5,3*LY/2+1e-5]))

# Integration methods
mim1 = gf.MeshIm(m1, 5)
mim2 = gf.MeshIm(m2, 5)

# Finite element methods
mf1 = gf.MeshFem(m1)
mf1.set_classical_fem(fem_order)
mf2 = gf.MeshFem(m2)
mf2.set_classical_fem(fem_order)

# Model
md = gf.Model("real")

md.add_initialized_data("epsilon", 1.)            # Emissivity
md.add_initialized_data("sigma", 5.670374419e-8) # Stefan-Boltzmann constat

md.add_filtered_fem_variable("q1", mf1, XP1_RG)

md.add_macro("temp_dist1", f"200+100*({LY}-X(2))")
md.add_macro("temp_dist2", f"500-100*({LY/2}-X(2))")

md.add_standard_secondary_domain('block2', mim2, XM2_RG)


#md.add_macro("R", "Norm(X-Secondary_Domain(X))") # distance
md.add_macro("R2", "Norm_sqr(X-Secondary_Domain(X))") # distance squared
md.add_macro("cos1", "-Normalized(X-Secondary_Domain(X)).Normal")
md.add_macro("cos2", "Normalized(X-Secondary_Domain(X)).Secondary_Domain(Normal)")
md.add_linear_twodomain_term(mim1, "(q1-(epsilon*sigma/pi)*pow(temp_dist1,4)*cos1*cos2/R2)*Test_q1", XP1_RG, "block2")

md.variable_list()

md.solve("noisy")

print(md.variable("q1"))
mf1.export_to_vtu("results1.vtu", mf1, md.interpolation("temp_dist1", mf1), "Temperature",
                                  md.mesh_fem_of_variable("q1"), md.variable("q1"), "Local heat flux 1 to entire body 2")
mf2.export_to_vtu("results2.vtu", mf2, md.interpolation("temp_dist2", mf1), "Temperature")

I get

message from gf_model_get follow:
List of model variables and data:
Variable       q1                             1 copy   fem dependant       21 doubles.	           
Data           epsilon                        1 copy   constant size        1 double.	           
Data           sigma                          1 copy   constant size        1 double.	           

[  5.46448542   6.44897111   7.89445256   9.86730785  13.00680725
  17.73696775  26.05052058  40.94033532  67.79433132 118.70124878
 182.86410139 215.55721606 206.92865459 177.20697728 142.91517871
 112.47833449  85.78935932  64.1720862   46.64474521  32.99400788
  22.57582633]

It seems to work but you need to check the if the result about the heat flux q1 (W/m²) from body 1 to body 2 seems to be correct.

Regarding switching source and target, I do not think this is possible within the same getfem model, because AFAIK a getfem model only accepts one secondary domain. The solution would be two different getfem models:
In model 1: body 1 is the main domain and body 2 is the secondary domain
In model 2: body 2 is the main domain and body 1 is the secondary domain
Then you will have to transfer data between the models yourself manually.

image
image617×733 8.22 KB

image
image617×733 11.9 KB

8

The script uses this term

"(q1-(epsilon*sigma/pi)*pow(temp_dist1,4)*cos1*cos2/R2)*Test_q1"

which should be the energy density emitted from every point on body 1 to the entire body 2.

If you wan to find the energy density received at every point on body 1 from the entire body 2, I guess you can simply use

"(q1-(epsilon*sigma/pi)*pow(temp_dist2,4)*cos1*cos2/R2)*Test_q1"

instead (temp_dist1 replaced with temp_dist2).

9

Dear Kostas,
thank you very much for your support.
For my original problem, I’m following a simplified and semi-analytical way that doesn’t need double domain integrals.
On the other hand, for learning purposes, I tried playing a little bit with your example.

Attached you can find the file “Secondary_domain.py” where I tried your last suggestion, and it works!
I wanted to check your sentence about the possible need of multiple models. In this connection, GetFEM allowed me to add more than one “secondary_domain” (block1 and block2) and I used the add_linear_twodomain_term twice, and it worked! producing the same result as above.

from 1 to 2
[0.70212207 1.04776306 1.26374925 1.37371115 1.43252793 1.47311906
1.5005264 1.52567848 1.54619877 1.56702349 1.58512817 1.60396862
1.61993824 1.6361611 1.64703334 1.65519785 1.64728979 1.61718148
1.52254198 1.29220201 0.8873764 ]

from 2 to 1
[439.11839996 645.47448948 767.40724927 822.41554411 845.3822193
857.11864553 860.70618004 862.84107881 862.1512047 861.52532879
859.29270801 857.37472738 853.86713974 850.42971007 844.24884312
836.666985 821.26160296 795.07352684 738.34796458 618.03079265
418.3229001 ]

So far so good!
Maybe I didn’t fully get your point about this issue. Can you be so kind to check and confirm? Thank you.

(Sorry, I could not resist in changing the mesh boundary regions detection by mean of outer_faces_with_direction function, because it looks clearer to me and gets rid of LX and WX geometrical parameters)

Attached you can find the file “Secondary_domain_ASM_GENERIC.py”
Since in the “User Documentation → Double domain integrals or terms” it is mentioned
“For the utilisation with the Python/Scilab/Octave/Matlab interface, see the documentation on gf_asm command and the model object.”
I tried to replace the add_linear_twodomain_term with asm_generic, introducing the “Secondary_domain” parameter. The results were somehow different in modulus.

from 1 to 2
[-0.19414747 -1.10378408 -0.68364226 -1.45983641 -0.76744577 -1.56992072
-0.80123091 -1.62689632 -0.82503144 -1.67120159 -0.84569987 -1.71059364
-0.86440554 -1.74466777 -0.87956035 -1.76382997 -0.88285342 -1.71811025
-0.82476995 -1.35973701 -0.24593961 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0.
0. 0. 0. 0. 0. 0. ]

from 2 to 1
[ 0. 0. 0. 0. 0.
0. 0. 0. 0. 0.
0. 0. 0. 0. 0.
0. 0. 0. 0. 0.
0. -121.60081757 -679.50096694 -415.49533025 -873.82547428
-453.0140109 -913.39067345 -459.62855813 -920.06250829 -460.05032345
-918.78896459 -458.45991637 -914.36348446 -455.62690356 -906.83239076
-450.82962093 -891.61027478 -440.05210878 -844.82109677 -399.6484188
-650.76450196 -115.78694512]

My further question is, what is the difference between the two methods above? I got the first because you already explained it to me, but what happens in asm_generic when I set a Secondary Domain?

Thank you,
Lorenzo

(Attachment Secondary_domain-ASM_GENERIC.py is missing)

(Attachment Secondary_domain.py is missing)

10

In the following text the Attachments that were refused with the previous subittal

### file Secondary_domain.py

import getfem as gf
import numpy as np

gf.util_trace_level(1)
gf.util_warning_level(1)

# Mesh size
NX = 5
NY = 10

LX = 2.       # Length
LY = 4.       # Height

WX = 0.5      # Gap length between block 1 an block 2

fem_order = 2

# Mesh generation
m1 = gf.Mesh("import", "structured",
             f"GT='GT_QK(2,1)';ORG=[0,0];SIZES=[{LX},{LY}];NSUBDIV=[{NX},{NY}]")
m2 = gf.Mesh("import", "structured",
             f"GT='GT_QK(2,1)';ORG=[{LX+WX},0];SIZES=[{LX},{LY}];NSUBDIV=[{NX},{NY}]")

# Mesh region definition
XM1_RG = 11
XM2_RG = 21
XP1_RG = 12
XP2_RG = 22
YM1_RG = 13
YM2_RG = 23
YP1_RG = 14
YP2_RG = 24

m1.set_region(XM1_RG, m1.outer_faces_with_direction([-1,0], 1))
m2.set_region(XM2_RG, m2.outer_faces_with_direction([-1,0], 1))

m1.set_region(XP1_RG, m1.outer_faces_with_direction([1,0], 1))
m2.set_region(XP2_RG, m2.outer_faces_with_direction([1,0], 1))

m1.set_region(YM1_RG, m1.outer_faces_with_direction([0,-1], 1))
m2.set_region(YM2_RG, m2.outer_faces_with_direction([0,-1], 1))

m1.set_region(YP1_RG, m1.outer_faces_with_direction([0,1], 1))
m2.set_region(YP2_RG, m2.outer_faces_with_direction([0,1], 1))

# Integration methods
mim1 = gf.MeshIm(m1, 5)
mim2 = gf.MeshIm(m2, 5)

# Finite element methods
mf1 = gf.MeshFem(m1)
mf1.set_classical_fem(fem_order)
mf2 = gf.MeshFem(m2)
mf2.set_classical_fem(fem_order)

# Model
md = gf.Model("real")

md.add_initialized_data("epsilon", 1.)            # Emissivity
md.add_initialized_data("sigma", 5.670374419e-8)  # Stefan-Boltzmann constat

md.add_filtered_fem_variable("q1", mf1, XP1_RG)
md.add_filtered_fem_variable("q2", mf2, XM2_RG)

md.add_macro("temp_dist1", f"100+1.5*X(2)")
md.add_macro("temp_dist2", f"500-1.5*X(2)")

mf1.export_to_vtu("results.vtu", mf1, md.interpolation("temp_dist1", mf1), "Temperature 1")

md.add_standard_secondary_domain('block2', mim2, XM2_RG)
md.add_standard_secondary_domain('block1', mim1, XP1_RG)

#md.add_macro("R", "Norm(X-Secondary_Domain(X))") # distance
md.add_macro("R2", "Norm_sqr(X-Secondary_Domain(X))") # distance squared
md.add_macro("cos1", "-Normalized(X-Secondary_Domain(X)).Normal")
md.add_macro("cos2", "Normalized(X-Secondary_Domain(X)).Secondary_Domain(Normal)")

md.add_linear_twodomain_term(mim1, "(q1-(epsilon*sigma/pi)*pow(temp_dist1,4)*cos1*cos2/R2)*Test_q1", XP1_RG, "block2")
#md.add_linear_twodomain_term(mim1, "(q1-(epsilon*sigma/pi)*pow(temp_dist2,4)*cos1*cos2/R2)*Test_q1", XP1_RG, "block2")
md.add_linear_twodomain_term(mim2, "(q2-(epsilon*sigma/pi)*pow(temp_dist2,4)*cos1*cos2/R2)*Test_q2", XM2_RG, "block1")
md.variable_list()

md.solve("noisy")

print("from 1 to 2"); print(md.variable("q1"))
print()
print("from 2 to 1"); print(md.variable("q2"))



### file Secondary_domain_ASM_GENERIC.py

import getfem as gf
import numpy as np

gf.util_trace_level(1)
gf.util_warning_level(1)

# Mesh size
NX = 5
NY = 10

LX = 2.       # Length
LY = 4.       # Height

WX = 0.5      # Gap length between block 1 an block 2

fem_order = 2

# Mesh generation
m1 = gf.Mesh("import", "structured",
             f"GT='GT_QK(2,1)';ORG=[0,0];SIZES=[{LX},{LY}];NSUBDIV=[{NX},{NY}]")
m2 = gf.Mesh("import", "structured",
             f"GT='GT_QK(2,1)';ORG=[{LX+WX},0];SIZES=[{LX},{LY}];NSUBDIV=[{NX},{NY}]")

# Mesh region definition
XM1_RG = 11
XM2_RG = 21
XP1_RG = 12
XP2_RG = 22
YM1_RG = 13
YM2_RG = 23
YP1_RG = 14
YP2_RG = 24

m1.set_region(XM1_RG, m1.outer_faces_with_direction([-1,0], 1))
m2.set_region(XM2_RG, m2.outer_faces_with_direction([-1,0], 1))

m1.set_region(XP1_RG, m1.outer_faces_with_direction([1,0], 1))
m2.set_region(XP2_RG, m2.outer_faces_with_direction([1,0], 1))

m1.set_region(YM1_RG, m1.outer_faces_with_direction([0,-1], 1))
m2.set_region(YM2_RG, m2.outer_faces_with_direction([0,-1], 1))

m1.set_region(YP1_RG, m1.outer_faces_with_direction([0,1], 1))
m2.set_region(YP2_RG, m2.outer_faces_with_direction([0,1], 1))

# Integration methods
mim1 = gf.MeshIm(m1, 5)
mim2 = gf.MeshIm(m2, 5)

# Finite element methods
mf1 = gf.MeshFem(m1)
mf1.set_classical_fem(fem_order)
mf2 = gf.MeshFem(m2)
mf2.set_classical_fem(fem_order)

# Model
md = gf.Model("real")

md.add_initialized_data("epsilon", 1.)            # Emissivity
md.add_initialized_data("sigma", 5.670374419e-8)  # Stefan-Boltzmann constat

md.add_filtered_fem_variable("q1", mf1, XP1_RG)
md.add_filtered_fem_variable("q2", mf2, XM2_RG)

md.add_macro("temp_dist1", f"100+1.5*X(2)")
md.add_macro("temp_dist2", f"500-1.5*X(2)")

mf1.export_to_vtu("results.vtu", mf1, md.interpolation("temp_dist1", mf1), "Temperature 1")

md.add_standard_secondary_domain('block2', mim2, XM2_RG)
md.add_standard_secondary_domain('block1', mim1, XP1_RG)

#md.add_macro("R", "Norm(X-Secondary_Domain(X))") # distance
md.add_macro("R2", "Norm_sqr(X-Secondary_Domain(X))") # distance squared
md.add_macro("cos1", "-Normalized(X-Secondary_Domain(X)).Normal")
md.add_macro("cos2", "Normalized(X-Secondary_Domain(X)).Secondary_Domain(Normal)")

#md.add_linear_twodomain_term(mim1, "(q1-(epsilon*sigma/pi)*pow(temp_dist1,4)*cos1*cos2/R2)*Test_q1", XP1_RG, "block2")
#md.add_linear_twodomain_term(mim1, "(q1-(epsilon*sigma/pi)*pow(temp_dist2,4)*cos1*cos2/R2)*Test_q1", XP1_RG, "block2")
#md.add_linear_twodomain_term(mim2, "(q2-(epsilon*sigma/pi)*pow(temp_dist2,4)*cos1*cos2/R2)*Test_q2", XM2_RG, "block1")
#print("from 1 to 2"); print(md.variable("q1"))
#print()
#print("from 2 to 1"); print(md.variable("q2"))

md.variable_list()

#md.solve("noisy")




# TEST ASM_GENERIC

result_asm_1 = gf.asm_generic(mim1, 1, "(q1-(epsilon*sigma/pi)*pow(temp_dist1,4)*cos1*cos2/R2)*Test_q1", XP1_RG, md, 'Secondary_domain', "block2")
print("from 1 to 2"); print(result_asm_1)
print()


result_asm_2 = gf.asm_generic(mim2, 1, "(q2-(epsilon*sigma/pi)*pow(temp_dist2,4)*cos1*cos2/R2)*Test_q2", XM2_RG, md, 'Secondary_domain', "block1")
print("from 2 to 1"); print(result_asm_2)
print()
11

oh nice! My bad, I was not aware that multiple secondary domains are actually supported.

Another note, please remove the copyright and author header on the original snippet you downloaded. It was not correct and copied by mistake.

12

No problem, I actually prefer outer_faces_with_direction.

The asm_generic command calculates the right hand side vector of the equations to solve, the model solve method gives you directly the result of the equations, i.e. the nodal values of the field in W/m² or whatever the units are. The values of the RHS vector itself are not anything very useful in general. In classical mechanics, people use the RHS as a “Force vector” and add forces to it, but this is not really a very good approach.