Creating SPH particles in CUBIT

Coreform Cubit
#1

Hi. I’m interested in created SPH particles for use in LS-DYNA in a fluid analysis. DYNA’s own tool does the job, but creates particles on a cartesian grid (i.e. on the points of hexagonal elements), which is fine except that spherical particles don’t pack that way under gravity. This means that once you start the simulation, you need to add more particles to actually fill the domain as particles settle into their close-packed hexagonal configuration.

So, I’m interested to see if CUBIT can create nodes in a hexagonal arrangement to minimize the settling. The nodes of a regular, uniform tet mesh ought to do the trick. But a typical tet meshing operation will not produce a regular uniform tet mesh. I’m thinking about an approach where a bounding box is filled with regular tets and then nodes outside the domain of interest (an .stl shape) are deleted, like a sculpting process.

Any thoughts on how to do this?
Thanks

#2

This paper speaks to the issue. I’m aiming for the HCP configuration

http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.756.9006&rep=rep1&type=pdf

#3

Welcome @JeffBodner!

I had some old code lying around that was close, so I adapted it slightly. Notice that here I use the creation of CAD vertices to represent the particle locations. This is probably not the most efficient approach, as directly creating mesh nodes is likely more efficient. But hopefully this is enough to get you started!

piston.stp (56.4 KB)

import math
import numpy

def compute_fcc_packing_in_volume( targetVolumeID, sphereDiameter, bbox ):
  sphereVertexIDs = initialize_fcc_packing( sphereDiameter, bbox )
  targetBody = cubit.body( targetVolumeID )
  for v in sphereVertexIDs:
    vertexCoordinates = cubit.vertex( v ).coordinates()
    vertexInVolume = targetBody.point_containment( vertexCoordinates )
    if vertexInVolume <= 0:
      cubit.cmd( f"delete vertex {v}" )

def initialize_fcc_packing( sphereDiameter, bbox ):
  create_fcc_background_grid( sphereDiameter, bbox )
  latticeVolumeID = cubit.get_last_id( "volume" )
  faceNodeIDs = get_element_faces( latticeVolumeID )
  faceNodeIDs = extract_unique_element_faces( faceNodeIDs )
  beginVertexID = cubit.get_last_id( "vertex" ) + 1
  add_vertices_at_mesh_nodes( latticeVolumeID )
  add_vertices_at_face_centers( latticeVolumeID )
  endVertexID = cubit.get_last_id( "vertex" )
  sphereVertexIDs = list( range( beginVertexID, endVertexID + 1 ) )
  cubit.cmd( f"delete volume {latticeVolumeID}" )
  return sphereVertexIDs

def create_fcc_background_grid( sphereDiameter, bbox ):
  latticeConstant = 2 * math.sqrt( 2.0 ) * ( sphereDiameter / 2.0 )
  xLength = bbox[3] - bbox[0]
  yLength = bbox[4] - bbox[1]
  zLength = bbox[5] - bbox[2]
  nX = math.ceil( xLength / latticeConstant )
  nY = math.ceil( yLength / latticeConstant )
  nZ = math.ceil( zLength / latticeConstant )
  xLength = nX * latticeConstant
  yLength = nY * latticeConstant
  zLength = nZ * latticeConstant
  xCenter = ( bbox[0] + bbox[3] ) / 2.0
  yCenter = ( bbox[1] + bbox[4] ) / 2.0
  zCenter = ( bbox[2] + bbox[5] ) / 2.0
  cubit.cmd( f"brick x {xLength} y {yLength} z {zLength}" )
  latticeVolumeID = cubit.get_last_id( "volume" )
  cubit.cmd( f"move volume {latticeVolumeID} x {xCenter} y {yCenter} z {zCenter}" )
  cubit.cmd( f"volume {latticeVolumeID} scheme map" )
  cubit.cmd( f"volume {latticeVolumeID} size {latticeConstant}" )
  cubit.cmd( f"mesh volume {latticeVolumeID}" )

def get_element_face_nodal_connectivity( elementID ):
  elemNodeIDs = cubit.get_connectivity( "hex", elementID )
  nConn = [ [ elemNodeIDs[ 2 ], elemNodeIDs[ 1 ], elemNodeIDs[ 0 ], elemNodeIDs[ 3 ] ],
            [ elemNodeIDs[ 2 ], elemNodeIDs[ 3 ], elemNodeIDs[ 7 ], elemNodeIDs[ 6 ] ],
            [ elemNodeIDs[ 2 ], elemNodeIDs[ 6 ], elemNodeIDs[ 5 ], elemNodeIDs[ 1 ] ],
            [ elemNodeIDs[ 4 ], elemNodeIDs[ 5 ], elemNodeIDs[ 6 ], elemNodeIDs[ 7 ] ],
            [ elemNodeIDs[ 4 ], elemNodeIDs[ 7 ], elemNodeIDs[ 3 ], elemNodeIDs[ 0 ] ],
            [ elemNodeIDs[ 4 ], elemNodeIDs[ 0 ], elemNodeIDs[ 1 ], elemNodeIDs[ 5 ] ] ]
  return nConn
  
def get_element_faces( volumeID ):
  elemIDs = cubit.parse_cubit_list( "hex", f"in volume {volumeID}" )
  faceNodeIDs = [ ]
  for e in elemIDs:
    elemFaceNodeIDs = get_element_face_nodal_connectivity( e )
    for f in range( 0, len( elemFaceNodeIDs ) ):
      faceNodeIDs.append( elemFaceNodeIDs[ f ] )
  return numpy.array( faceNodeIDs )

def extract_unique_element_faces( faceNodeIDs ):
  tmp = faceNodeIDs
  for i in range( 0, len( tmp ) ):
    tmp[i].sort()
  _, indices = numpy.unique(faceNodeIDs, axis = 0, return_index = True)
  return faceNodeIDs[ indices ]

def compute_node_list_centroid( nodeIDList ):
  centroid = numpy.array( [ 0.0, 0.0, 0.0 ] )
  for n in nodeIDList:
    centroid += cubit.get_nodal_coordinates( int( n ) )
  centroid /= len( nodeIDList )
  return centroid    

def add_vertices_at_mesh_nodes( volumeID ):
  nodeIDs = cubit.parse_cubit_list( "node", f"in volume {volumeID}" )
  for n in nodeIDs:
    x, y, z = cubit.get_nodal_coordinates( int( n ) )
    cubit.cmd( f"create vertex location {x} {y} {z}" )

def add_vertices_at_face_centers( volumeID ):
  faceNodeIDs = get_element_faces( volumeID )
  faceNodeIDs = extract_unique_element_faces( faceNodeIDs )
  for nodeIDList in faceNodeIDs:
    x, y, z = compute_node_list_centroid( nodeIDList )
    cubit.cmd( f"create vertex location {x} {y} {z}" )

def list_to_str( input_str ):
  return " ".join( [ str( val ) for val in input_str ] )

###### SCRIPT ######
cubit.cmd( "reset" )
cubit.cmd( "echo off" )
cubit.cmd( "undo off" )
cubit.cmd( "journal off" )
cubit.cmd( "info off" )

cubit.cmd( "graphics perspective off ")
cubit.cmd( "graphics parallel scale  7.0398914 ")
cubit.cmd( "from  6.4723108 -17.271496 -16.445939 ")
cubit.cmd( "at            0      -1.75          0 ")
cubit.cmd( "up   -0.084038083 -0.74097628  0.66625201 ")
cubit.cmd( "geometry visibility on" )

cubit.cmd('import step "piston.stp" heal')
targetVolumeID = 1
sphereDiameter = 0.25
bbox = cubit.volume( targetVolumeID ).bounding_box()

compute_fcc_packing_in_volume( targetVolumeID, sphereDiameter, bbox )
cubit.cmd("delete volume 1")
cubit.cmd("block 1 vertex all")
cubit.cmd(f"vertex all size {sphereDiameter}")
cubit.cmd("mesh vertex all")

cubit.cmd( "info on" )
cubit.cmd( "echo on" )
cubit.cmd( "undo on" )
cubit.cmd( "journal on" )
cubit.cmd( "geometry visibility on" )
#4

Hi Jeff,

The following method is not how I would create this for production code. It is not scalable for a large number of nodes, but it is a nice demonstration of the math. In production, I would just create nodes directly and not worry about the tetrahedrons, but they are hard to visualize.

I talked with Greg and he said having more examples is always good. This is just another way of looking at the problem and it let me easily visualize the problem along the way.

I create a single tetrahedron with unit length sides. I then copy the volume to create the pattern
for the base. Then I reflect the base about a top tetrahedron vertex.
image image

I then just replicate this pattern through the rest of the layers and imprint and merge all the vertices.
image

Here I created unit spheres at every vertex.

image

Below is the sample code to create the equilateral tetrahedra (edge length = 1). Given the code here it wouldn’t be that hard to convert this to just creating nodes.

Karl

from math import *
cubit.cmd('reset')
cubit.cmd('undo off')

cubit.cmd( f'create pyramid height {sqrt(6)/3} sides 3 radius {.5/cos(radians(30))} top 0 ')
cubit.cmd('rotate Volume 1  angle 30  about Z include_merged ')
cubit.cmd('move Vertex 3  location 0 0 0 include_merged ')
cubit.cmd('Volume 1  copy move x 1 repeat 9 ')
cubit.cmd(f'Volume 1 to 10 copy move x {cos(radians(60))} y {sin(radians(60))} ')
cubit.cmd(f'Volume 1 to 20 copy move y {2*sin(radians(60))} repeat 4')

vertex_4 = cubit.get_center_point("vertex", 4)
cubit.cmd(f'volume 1 to 100 copy reflect origin {vertex_4[0]}   {vertex_4[1]}  {vertex_4[2]} direction 0 0 1')
cubit.cmd(f'volume 1 to 200 copy move z {2*sqrt(6)/3} repeat 4')

cubit.cmd('imprint all'')
cubit.cmd('merge vertex all')
cubit.cmd('mesh vertex all')

cubit.cmd('block 1 node all')
cubit.cmd('block 1 element type SPHERE')

cubit.cmd('undo on')