Smooth curved 2nd-order elements with narrow boundary layers

[REardleyBrunt]

I’m trying to generate a 2nd order (hex20) mesh with an interior “fluid” region surrounded by an exterior “solid” region separated by a curved surface, and require extremely narrow boundary layers on both sides of the interface to resolve the physics involved, and this generally requires high aspect ratio elements (~700 in the extreme) to avoid excessive element counts. The high-order solver this will be used with should be robust to curved elements and narrow boundary layers. The issue I’m having is that while I can get the 2nd order mid-nodes on the boundary to conform to the geometry (which requires setting set node constraint on), every other mid-node just forms a straight line with the other nodes of the edge rather than getting curvature, and this leads to elements on the outside of the curve intersecting themselves.

Without set node constraint on the elements at the boundary have no curvature (presumably the default smart setting is detecting the poor element shapes and aborting placement of mid-nodes on the curved surface):

Is there a way to smooth the placement of mid-nodes away from the curved interface? Ideally I would want the whole element to be deformed following the surface in the first layer, and as the boundary layers expand this would be smoothed out towards elements without curvature away from the surface.

The geometry I’m meshing currently is a cylindrical fluid region with annular solid region, so I’m generating a surface mesh on the end surfaces with boundary layers + pave and sweeping in the axial direction, so the but ideally a solution to this would also be applicable in an arbitrary 3D geometry with curvature.

An example script:

import sys
sys.path.append('/opt/Coreform-Cubit-2023.11/bin')
import cubit

### User Settings ###

L_r_fluid = 1           # pipe interior radius
L_r_solid = 1.1         # pipe exterior radius
L_z = 5                 # pipe length

delta_z_axial = 1e-1
delta_xy_core = 1e-1
fluid_first_row = 1e-3
fluid_growth_factor = 1.5
fluid_num_layers = 12
solid_first_row = 1e-3
solid_growth_factor = 1.8
solid_num_layers = 6

### Generate Geometry ###

cubit.cmd(f'cylinder height {L_z} radius {L_r_fluid}') # fluid region
cubit.cmd(f'cylinder height {L_z} radius {L_r_solid}') # solid+fluid region
cubit.cmd('subtract volume 1 from volume 2 keep_tool') # create solid-only region

# Imprint and merge
cubit.cmd('imprint volume all')
cubit.cmd('merge volume all')
cubit.cmd('compress')

### Set up boundary layers ###

# Create boundary layer in fluid region
cubit.cmd('create boundary_layer 1')
cubit.cmd(f'modify boundary_layer 1 uniform height {fluid_first_row} growth {fluid_growth_factor} layers {fluid_num_layers}')
cubit.cmd('modify boundary_layer 1 add surface 1 volume 1')
cubit.cmd('modify boundary_layer 1 continuity off')

# Create boundary layer in solid region
cubit.cmd('create boundary_layer 2')
cubit.cmd(f'modify boundary_layer 2 uniform height {solid_first_row} growth {solid_growth_factor} layers {solid_num_layers}')
cubit.cmd('modify boundary_layer 2 add surface 1 volume 2')
cubit.cmd('modify boundary_layer 2 continuity off')

### Generate Mesh ###

# create element blocks
cubit.cmd('block 1 add volume 1')
cubit.cmd('block 2 add volume 2')

# set element type
cubit.cmd('block 1 element type hex20')
cubit.cmd('block 2 element type hex20')

# force nodes to follow curved surfaces
cubit.cmd('set node constraint on')

# Set mesh size for axial resolution (on all volumes)
cubit.cmd(f'volume 1 2 size {delta_z_axial}')

# Set approximate mesh size for core (on fluid & solid inlet surfaces)

cubit.cmd(f'surface 2 5 size {delta_xy_core}')

# Mesh fluid inlet

cubit.cmd('surface 2 scheme pave')
cubit.cmd('mesh surface 2')

# Mesh solid inlet

cubit.cmd('surface 5 scheme pave')
cubit.cmd('mesh surface 5')

# Sweep mesh through volume

cubit.cmd('volume 1 2 redistribute nodes off')
cubit.cmd('volume 1 scheme Sweep source surface 2 target surface 3 sweep transform least squares')
cubit.cmd('volume 2 scheme Sweep source surface 5 target surface 6 sweep transform least squares')
cubit.cmd('volume 1 2 autosmooth target on fixed imprints off smart smooth off')
cubit.cmd('mesh volume 1')
cubit.cmd('mesh volume 2')

[strieben]

Hello,

I have posted the “same” or quite similar question already 10 years ago (cf. Meshing of a pipe for multi-physics simulations - Coreform Cubit / Meshing - Coreform LLC). I never got any reply!

I ended up doing this myself for perfectly cylindrical pipes by making tetrahedra out of the wedges or hexas in the boundary layer using Netgen and then adding curved mid-side nodes using a Python script and the VTK API.

My benchmark problem was to integrate the flow velocity field for semi-analytical Schlichting flow profiles at different Reynolds numbers (i.e. the higher the Reynolds number the steeper the gradient near the pipe wall that has to be resolved) over the 3D volume of the pipe and to compute the deviation from the theoretical mean velocity. With the fully rounded boundary layer I was able to integrate flow profiles “exactly” (i.e. with Deviation < 0.025% which is about one forth of the effects I am looking for) up to one order of magnitude more (Re=10^6) than with just the outermost boundary layer rounded to the cylindrical boundary (default behavior of Cubit, Re=10^5).

The ripple effects at lower Re numbers come from the implementation of the semi-analytic Schliching profiles, which are not perfectly smooth parabolic but piece-wise liner profiles. At higher Re numbers the overprediction of the mean velocities can be easily understood, because no velocity samples are taken in the high gradient region near the pipe wall. Only high velocities towards the middle of the pipe will be sampled!

It is obvious, that the deviations at high Re numbers will be large, if there are no suitable boundary layers for the velocity gradients at the wall at all (green graph). If one uses the default boundary layers of most meshing tools, where only mid-side nodes on the pipe wall will be placed on the curved geometry, but all other boundary layer elements will have straight edges, the situation improves but already above Re=10^5 the deviations grow again at a computing time of 15 minutes for 20 layers (pink graph). By smoothing all the 20 layers the deviation only become noticeable over Re=10^6 at the same 15 minutes computing time (blue graph). By reducing the the number of smoothed boundary layers and by playing a little bit with the bias or growth factor between layers, one can reduce the computing time by a factor of three while maintaining almost the same accuracy (orange graph).

[REardleyBrunt]

Thank you, this is a useful post to be aware of, the figure shows nicely what I’m after in terms of the hex20 elements. In my case it’s the hex20 elements themselves that I’m after, I think it should be feasible to add the curvature manually for a cylindrical geometry but it gets considerably more complicated for an arbitrary curved surface since you then can’t rely on a simple known function to describe the geometry. Any more tips on how you added the curved mid-side nodes would be greatly appreciated. Hopefully there will be a streamlined way to do this within cubit though.

[strieben]

I completely agree, that smoothed higher order boundary layers for arbitrary geometries would be a highly desirable feature. This would reduce the needed element count considerably!

But I am not aware of one single meshing tool, that currently provides these. The default in all the tools I know (e.g. COMSOL, ANSYS, STARCCM only supports linear FE types in its internal mesher, even Gmsh, etc.) is to put only midside nodes on the boundary onto the exact geometry. This produces exactly the self-intersection errors, that you encountered.

Perhaps, the marketing department of Coreform should conduct a market study on this topic? The IGA direction in which Coreform seems to be moving would justify such a feature. The implementation of it for arbitrary element types (mostly wedges and hexas but in my case also tetras) is, however, non-trivial, I guess…

[Norbert_Hofbauer]

Hi @REardleyBrunt and @strieben,
sorry for the delay, it took some time to write a script for this.

First of all a screenshot to show you that it’s possible to place the midnodes.

and the same just in paraview

So a short explanation how this can be done. I started with your example script from above.
Then i choose the volume and the boundary surface for which i want the adjust the hex nodes.
To get my wanted hexes near the boundary i can define a distance. The distance is used to check if the center point from a hex is near enough to the surface.

volume_id = 1
surface_id = 1
distance = 2.5e-1
hex_ids = get_elements(volume_id,surface_id,distance)

The get_elements checks the distance from the hexes and put them into a list if they are in range.
With the above values i can get the hexes that are within the boundary layer for volume 1.

When i got the hexes we can already process them. I turn off some output to speed things up. Otherwise the command line output will slow us down greatly.

cubit.cmd(f"info off")
cubit.cmd(f"warning off")
cubit.cmd(f"graphics off")
adjust_elements(hex_ids,surface_id,volume_id)
cubit.cmd(f"info on")
cubit.cmd(f"warning on")
cubit.cmd(f"graphics on")

The adjust_elements loops over all hexes. It computes the side of the hex that’s nearest to the surface and the average distance of the nodes in the vertices from that side. This average distance will be used to move the midnodes to the same distance to the surface. Only the nearest side of the hex will be adjusted, not the opposite side too.

def adjust_elements(hex_ids,surface_id,volume_id):
 cubit.cmd("set node constraint off")
 #check if hex20
 if (len(hex_ids) > 0):
  conn = cubit.get_expanded_connectivity("hex",hex_ids[0])
  if not (len(conn)==20):
   print("Not HEX20!")
   return False
 #get surface
 surface = cubit.surface(surface_id)
 #loop over hexes
 max_hexes = len(hex_ids)
 current_hex = 0
 counter = 0
 for hid in hex_ids:
   conn = cubit.get_expanded_connectivity("hex",hid)
   #adjust side
   side,distance = get_side_and_distance(conn,surface)
   adjust_side_distance(side,distance,conn,surface,volume_id)
   current_hex = current_hex + 1
   counter = counter + 1
   if counter==500:
    cubit.cmd(f"info on")
    print(f"adjusted {current_hex} hexes from {max_hexes} for volume {volume_id} and surface {surface_id}")
    counter = 0
    cubit.cmd(f"info off")

Here is the complete journal. Hopefully that’s what you were looking for. Also please note that’s it’s only working for hex20.

#!python
cubit.cmd("reset")

import sys
sys.path.append('/opt/Coreform-Cubit-2023.11/bin')
import cubit

### User Settings ###

L_r_fluid = 1           # pipe interior radius
L_r_solid = 1.1         # pipe exterior radius
L_z = 5                 # pipe length

delta_z_axial = 1e-1
delta_xy_core = 1e-1
fluid_first_row = 1e-3
fluid_growth_factor = 1.5
fluid_num_layers = 12
solid_first_row = 1e-3
solid_growth_factor = 1.8
solid_num_layers = 6

### Generate Geometry ###

cubit.cmd(f'cylinder height {L_z} radius {L_r_fluid}') # fluid region
cubit.cmd(f'cylinder height {L_z} radius {L_r_solid}') # solid+fluid region
cubit.cmd('subtract volume 1 from volume 2 keep_tool') # create solid-only region

# Imprint and merge
cubit.cmd('imprint volume all')
cubit.cmd('merge volume all')
cubit.cmd('compress')

### Set up boundary layers ###

# Create boundary layer in fluid region
cubit.cmd('create boundary_layer 1')
cubit.cmd(f'modify boundary_layer 1 uniform height {fluid_first_row} growth {fluid_growth_factor} layers {fluid_num_layers}')
cubit.cmd('modify boundary_layer 1 add surface 1 volume 1')
cubit.cmd('modify boundary_layer 1 continuity off')

# Create boundary layer in solid region
cubit.cmd('create boundary_layer 2')
cubit.cmd(f'modify boundary_layer 2 uniform height {solid_first_row} growth {solid_growth_factor} layers {solid_num_layers}')
cubit.cmd('modify boundary_layer 2 add surface 1 volume 2')
cubit.cmd('modify boundary_layer 2 continuity off')

### Generate Mesh ###

# create element blocks
cubit.cmd('block 1 add volume 1')
cubit.cmd('block 2 add volume 2')

# set element type
cubit.cmd('block 1 element type hex20')
cubit.cmd('block 2 element type hex20')

# force nodes to follow curved surfaces
cubit.cmd('set node constraint on')

# Set mesh size for axial resolution (on all volumes)
cubit.cmd(f'volume 1 2 size {delta_z_axial}')

# Set approximate mesh size for core (on fluid & solid inlet surfaces)

cubit.cmd(f'surface 2 5 size {delta_xy_core}')

# Mesh fluid inlet

cubit.cmd('surface 2 scheme pave')
cubit.cmd('mesh surface 2')

# Mesh solid inlet

cubit.cmd('surface 5 scheme pave')
cubit.cmd('mesh surface 5')

# Sweep mesh through volume

cubit.cmd('volume 1 2 redistribute nodes off')
cubit.cmd('volume 1 scheme Sweep source surface 2 target surface 3 sweep transform least squares')
cubit.cmd('volume 2 scheme Sweep source surface 5 target surface 6 sweep transform least squares')
cubit.cmd('volume 1 2 autosmooth target on fixed imprints off smart smooth off')
cubit.cmd('mesh volume 1')
cubit.cmd('mesh volume 2')

import numpy

# get all elements from a volume wich are in range from the surface
def get_elements(volume_id,surface_id,distance):
 #hex_ids = cubit.parse_cubit_list( "hex", f"all in volume {volume_id} with z_coord>2.4" )
 hex_ids = cubit.parse_cubit_list( "hex", f"all in volume {volume_id}" )
 surface = cubit.surface(surface_id)
 return_hex_ids=[]
 for hid in hex_ids:
   location = numpy.array(cubit.get_center_point("hex",hid))
   new_location = numpy.array(surface.closest_point_trimmed(location))
   delta = new_location - location
   current_distance = numpy.linalg.norm(delta)
   if current_distance<=distance:
    return_hex_ids.append(hid)
 #print(f"volume {volume_id} surface {surface_id}  distance {distance}")
 return return_hex_ids

#calc hex vertice distance to surface per side
def get_side_distance(side,connectivity,surface):
 distance = 0
 node_ids=[]
 if side == 1:
  node_ids.append(connectivity[0])
  node_ids.append(connectivity[1])
  node_ids.append(connectivity[4])
  node_ids.append(connectivity[5])
 elif side == 2:
  node_ids.append(connectivity[1])
  node_ids.append(connectivity[2])
  node_ids.append(connectivity[5])
  node_ids.append(connectivity[6])
 elif side == 3:
  node_ids.append(connectivity[2])
  node_ids.append(connectivity[3])
  node_ids.append(connectivity[6])
  node_ids.append(connectivity[7])
 elif side == 4:
  node_ids.append(connectivity[0])
  node_ids.append(connectivity[3])
  node_ids.append(connectivity[4])
  node_ids.append(connectivity[7])
 elif side == 5:
  node_ids.append(connectivity[0])
  node_ids.append(connectivity[1])
  node_ids.append(connectivity[2])
  node_ids.append(connectivity[3])
 elif side == 6:
  node_ids.append(connectivity[4])
  node_ids.append(connectivity[5])
  node_ids.append(connectivity[6])
  node_ids.append(connectivity[7])

 for nid in node_ids:
   location = numpy.array(cubit.get_nodal_coordinates(nid))
   new_location = numpy.array( surface.closest_point_trimmed(location) )
   delta = new_location - location
   current_distance = numpy.linalg.norm(delta)
   distance = distance + current_distance
 distance = distance/4
 return distance

#get side number that's facing the surface
def get_side_and_distance(connectivity,surface):
 distance = []
 distance.append(get_side_distance(1,connectivity,surface))
 distance.append(get_side_distance(2,connectivity,surface))
 distance.append(get_side_distance(3,connectivity,surface))
 distance.append(get_side_distance(4,connectivity,surface))
 distance.append(get_side_distance(5,connectivity,surface))
 distance.append(get_side_distance(6,connectivity,surface))
 side = distance.index(min(distance))+1
 return side,min(distance)
 
#calc adjust middle nodes distance to surface per side
def adjust_side_distance(side,distance,connectivity,surface,volume_id):
 node_ids=[]
 if side == 1:
  node_ids.append(connectivity[8])
  node_ids.append(connectivity[12])
  node_ids.append(connectivity[13])
  node_ids.append(connectivity[16])
 elif side == 2:
  node_ids.append(connectivity[9])
  node_ids.append(connectivity[13])
  node_ids.append(connectivity[14])
  node_ids.append(connectivity[17])
 elif side == 3:
  node_ids.append(connectivity[10])
  node_ids.append(connectivity[14])
  node_ids.append(connectivity[15])
  node_ids.append(connectivity[18])
 elif side == 4:
  node_ids.append(connectivity[11])
  node_ids.append(connectivity[12])
  node_ids.append(connectivity[15])
  node_ids.append(connectivity[19])
 elif side == 5:
  node_ids.append(connectivity[8])
  node_ids.append(connectivity[9])
  node_ids.append(connectivity[10])
  node_ids.append(connectivity[11])
 elif side == 6:
  node_ids.append(connectivity[16])
  node_ids.append(connectivity[17])
  node_ids.append(connectivity[18])
  node_ids.append(connectivity[19])

 for nid in node_ids:
   skip = False
   location = numpy.array(cubit.get_nodal_coordinates(nid))
   new_location = numpy.array(surface.closest_point_trimmed(location))
   current_delta = new_location - location   
   if distance > 0:
    if not ((current_delta[0]==0) and (current_delta[0]==0) and (current_delta[0]==0)):
     unit_delta = current_delta / numpy.linalg.norm(current_delta)
     contained = cubit.is_point_contained("volume",volume_id,location)
     if (contained==1) or (contained==2):
      delta = current_delta - unit_delta*distance
     elif contained==0:
      delta = current_delta + unit_delta*distance
     print(f"node_id {nid} distance {distance} unit_delta {unit_delta}")
     print(f"current_delta {current_delta} delta {delta}")
    else:
     skip = True
   else:
    delta = current_delta
   if not skip:
    cubit.cmd( f"node {nid} move x {delta[0]} y {delta[1]} z {delta[2]}")
 return True

# move middle nodes from hex20 to align to a surface
def adjust_elements(hex_ids,surface_id,volume_id):
 cubit.cmd("set node constraint off")
 #check if hex20
 if (len(hex_ids) > 0):
  conn = cubit.get_expanded_connectivity("hex",hex_ids[0])
  if not (len(conn)==20):
   print("Not HEX20!")
   return False
 #get surface
 surface = cubit.surface(surface_id)
 #loop over hexes
 max_hexes = len(hex_ids)
 current_hex = 0
 counter = 0
 for hid in hex_ids:
   conn = cubit.get_expanded_connectivity("hex",hid)
   #adjust side
   side,distance = get_side_and_distance(conn,surface)
   adjust_side_distance(side,distance,conn,surface,volume_id)
   current_hex = current_hex + 1
   counter = counter + 1
   if counter==500:
    cubit.cmd(f"info on")
    print(f"adjusted {current_hex} hexes from {max_hexes} for volume {volume_id} and surface {surface_id}")
    counter = 0
    cubit.cmd(f"info off")

 return True

volume_id = 1
surface_id = 1
distance = 2.5e-1
hex_ids = get_elements(volume_id,surface_id,distance)

cubit.cmd(f"info off")
cubit.cmd(f"warning off")
cubit.cmd(f"graphics off")
adjust_elements(hex_ids,surface_id,volume_id)
cubit.cmd(f"info on")
cubit.cmd(f"warning on")
cubit.cmd(f"graphics on")
#cubit.cmd(f"select hex {' '.join(str(hid) for hid in hex_ids)}")

volume_id = 2
surface_id = 1
distance = 5e-2
hex_ids = get_elements(volume_id,surface_id,distance)

cubit.cmd(f"info off")
cubit.cmd(f"warning off")
cubit.cmd(f"graphics off")
adjust_elements(hex_ids,surface_id,volume_id)
cubit.cmd(f"info on")
cubit.cmd(f"warning on")
cubit.cmd(f"graphics on")
#cubit.cmd(f"select hex {' '.join(str(hid) for hid in hex_ids)}")

#cubit.cmd(f"export mesh 'smoothed_boundary.e' overwrite")

[REardleyBrunt]

Hi @Norbert_Hofbauer, thank you very much - this works nicely, I was thinking it should be possible to manually adjust the node positions like this but since I’m not particularly familar with how mesh node manipulations work in cubit this would have taken me a very long time so this is a fantastic solution to have. It runs for me and seems to work nicely so I’ll put it in action. Thanks again!