BHE simulation in a more complicated geometry

saturngreen67 · March 25, 2024, 6:39pm

Hi everyone,

I am trying to do a BHE simulation in a complicated environment. It includes 2 layers with a cavern inside the second layer (made by a Boolean operator) but every time I make a mesh I encounter the following error when I try to run the BHE simulation:

critical: /project/ProcessLib/HeatTransportBHE/LocalAssemblers/LocalDataInitializer.h:111 operator()()
error: You are trying to build a local assembler for an unknown mesh element type (N7MeshLib15TemplateElementINS_8TriRule3EEE). Maybe you have disabled this mesh element type in your build configuration, or a mesh element order does not match shape function order given in the project file.

Can anyone tell me how can I tackle this issue?

Best regards,
Farid

dmitri.naumov · March 25, 2024, 7:06pm

Hi Farid,

it looks like you have a 2D triangle elements in your 3D domain, could it be? If so, try to remove these boundary/interface elements and maybe it will already work.

– d

saturngreen67 · March 25, 2024, 7:22pm

Thank you very much for your prompt response. I am fairly new to this stuff. But I don’t have enough time already. How can I find and remove them?

Here are the codes:

import numpy as np

import pandas as pd

import gmsh

import sys

import os

from pathlib import Path

gmsh.initialize()

gmsh.model.add(“t3”)

out_dir = Path(os.environ.get(“OGS_TESTRUNNER_OUT_DIR”, “_out”))

if not out_dir.exists():

out_dir.mkdir(parents=True)

exe_dir = os.environ.get(“OGS_BINARY_DIR”)

mesh file name

bhe_mesh_file_name = “test3”

#geometry parameters

width = 1300

length = 1300

depth = 700

bhe_depth = 1100

#element sizes

bhe_radius = 0.07

alpha = 6.134 # see Diersch et al. 2011 Part 2

delta = alpha * bhe_radius # meshsize at BHE and distance of the surrounding optimal mesh points

elem_size_corner = 40

elem_size_BHE_relax = 2

#Creating Earth layers

gmsh.model.occ.addPoint(0.0, 0.0, 0.0, elem_size_corner, 1)

gmsh.model.occ.addPoint(1300, 0.0, 0.0, elem_size_corner, 2)

gmsh.model.occ.addPoint(1300, 1300, 0.0, elem_size_corner, 3)

gmsh.model.occ.addPoint(0.0, 1300, 0.0, elem_size_corner, 4)

gmsh.model.occ.addLine(1, 4, 1)

gmsh.model.occ.addLine(2, 1, 2)

gmsh.model.occ.addLine(4, 3, 3)

gmsh.model.occ.addLine(3, 2, 4)

gmsh.model.occ.synchronize()

gmsh.model.occ.addCurveLoop([1, 3, 4, 2], 1)

gmsh.model.occ.addPlaneSurface([-1], 1)

gmsh.model.occ.synchronize()

gmsh.model.occ.addPoint(650, 650, 0, delta, 5)

gmsh.model.occ.addPoint(650, 649.6, 0, delta, 6)

gmsh.model.occ.addPoint(650, 650.4, 0, delta, 7)

gmsh.model.occ.addPoint(650.3464, 650.2, 0, delta, 8)

gmsh.model.occ.addPoint(649.6536, 650.2, 0, delta, 9)

gmsh.model.occ.addPoint(650.3464, 649.8, 0, delta, 10)

gmsh.model.occ.addPoint(649.6536, 649.8, 0, delta, 11)

gmsh.model.occ.synchronize()

gmsh.model.mesh.embed(0,[5],2,1)

gmsh.model.mesh.embed(0,[6],2,1)

gmsh.model.mesh.embed(0,[7],2,1)

gmsh.model.mesh.embed(0,[8],2,1)

gmsh.model.mesh.embed(0,[9],2,1)

gmsh.model.mesh.embed(0,[10],2,1)

gmsh.model.mesh.embed(0,[11],2,1)

R1 = gmsh.model.occ.extrude([(2, 1)], 0, 0, -depth, [6], [1], True)

gmsh.model.occ.synchronize()

gmsh.model.addPhysicalGroup(3, [R1[1][1]], 1)

R2 = gmsh.model.occ.extrude([(2, 6)], 0, 0, -1300, [12], [1], True)

gmsh.model.occ.synchronize()

#Creating Salt Cavern

gmsh.model.occ.addPoint(650, 650, -1000, delta, 105)

gmsh.model.occ.addPoint(650, 649.6, -1000, delta, 106)

gmsh.model.occ.addPoint(650, 650.4, -1000, delta, 107)

gmsh.model.occ.addPoint(650.3464, 650.2, -1000, delta, 108)

gmsh.model.occ.addPoint(649.6536, 650.2, -1000, delta, 109)

gmsh.model.occ.addPoint(650.3464, 649.8, -1000, delta, 110)

gmsh.model.occ.addPoint(649.6536, 649.8, -1000, delta, 111)

gmsh.model.occ.addPoint(650, 650, -1000, delta, 112)

gmsh.model.occ.addPoint(650, 625, -1000, delta, 113)

gmsh.model.occ.addPoint(650, 675, -1000, delta, 114)

gmsh.model.occ.addPoint(625, 650, -1000, delta, 115)

gmsh.model.occ.addPoint(675, 650, -1000, delta, 116)

gmsh.model.occ.synchronize()

gmsh.model.occ.addCircleArc(113,105,116,106)

gmsh.model.occ.addCircleArc(116,105,114,107)

gmsh.model.occ.addCircleArc(114,105,115,108)

gmsh.model.occ.addCircleArc(115,105,113,109)

gmsh.model.occ.synchronize()

gmsh.model.occ.addCurveLoop([106, 107, 108, 109], 101)

gmsh.model.occ.addPlaneSurface([101], 101)

gmsh.model.occ.synchronize()

gmsh.model.mesh.embed(0,[105],2,101)

gmsh.model.mesh.embed(0,[106],2,101)

gmsh.model.mesh.embed(0,[107],2,101)

gmsh.model.mesh.embed(0,[108],2,101)

gmsh.model.mesh.embed(0,[109],2,101)

gmsh.model.mesh.embed(0,[110],2,101)

gmsh.model.mesh.embed(0,[111],2,101)

gmsh.model.mesh.embed(0,[112],2,101)

gmsh.model.mesh.embed(0,[113],2,101)

gmsh.model.mesh.embed(0,[114],2,101)

gmsh.model.mesh.embed(0,[115],2,101)

gmsh.model.mesh.embed(0,[116],2,101)

gmsh.model.occ.synchronize()

cyc1=gmsh.model.occ.extrude([(2,101)], 0, 0, -200, [12],[1],True)

gmsh.model.geo.synchronize()

gmsh.model.occ.synchronize()

#gmsh.model.removeEntities([(3,3)], recursive=False)

cavern,other = gmsh.model.occ.fragment([R2[1]], [cyc1[1]])

#gmsh.model.occ.cut([R2[1]], [cyc1[1]])

gmsh.model.geo.synchronize()

gmsh.model.occ.synchronize()

#cyc2=gmsh.model.occ.extrude([(2,101)], 0, 0, -200, [12],[1],True)

gmsh.model.geo.synchronize()

gmsh.model.occ.synchronize()

gmsh.model.addPhysicalGroup(3, [3], 3)

gmsh.model.geo.synchronize()

gmsh.model.occ.synchronize()

gmsh.model.addPhysicalGroup(3, [4], 4)

gmsh.model.geo.synchronize()

gmsh.model.occ.synchronize()

#Adding BHE to the model

bhe = gmsh.model.occ.extrude([(0, 5)], 0, 0, -bhe_depth, [110], [1], True)

gmsh.model.occ.synchronize()

gmsh.model.addPhysicalGroup(1, [bhe[1][1]], 501)

gmsh.model.addPhysicalGroup(2,[1], name=‘top’)

gmsh.model.addPhysicalGroup(2,[111], name=‘bottom’)

gmsh.model.occ.synchronize()

gmsh.option.setNumber(“Mesh.MshFileVersion”, 2.2)

gmsh.model.mesh.generate(3)

gmsh.model.mesh.removeDuplicateNodes()

gmsh.write(f"{out_dir}/{bhe_mesh_file_name}.msh")

if “CI” not in os.environ:

gmsh.fltk.run()

gmsh.finalize()

saturngreen67 · March 26, 2024, 8:06am

test10-codes.doc (45.5 KB)

dmitri.naumov · March 26, 2024, 2:22pm

Hi Farid, if you would upload your vtu file, I could take a look. There is a utility we use to remove certain elements: Remove Mesh Elements — take a look. Also the checkMesh utility helps with analysis.

– d

saturngreen67 · March 26, 2024, 2:44pm

Dear Dmitri,

Thank you very much for your response. Here is the link to the VTU file:

M_Jaeschke · March 26, 2024, 7:21pm

Dear Farid,

thanks for your interest to the Heat_Transport_BHE-Process in OGS and the interesting model. I looked at your creation und generation of the geometry/mesh and I see some critical points. The BHE nodes should be on already existing soil nodes, your extrusion of the BHE don’t fulfill this specification. This could be the reason, why OGS detected an unknown mesh element type.

How the mesh is translated to OGS? Do you use GMSH2OGS or msh2vtu? For mesh2vtu see this small BHE–example in ogstools - from my point of view mesh2vtu is the better way to do it.

I hope my thoughts will help!

Best regards,
Max

saturngreen67 · March 27, 2024, 9:11am

Dear Max,

Thank you for your response. I used GMSH2OGS but after I got your message, I tried it again with mesh2vtu as well but still couldn’t run the simulation. It appears that (only) in the BHE module of OGS, only volumes created by extrusion of a 2D surface can work. So I tried it with a little simpler geometry but more volumes and it worked.

Many thanks to both you and Dmitri for your support.

Best regards,
Farid