I have a 2D rectangular tube filled with particles vtu.mfx (13.7 KB)
I set the boundary condition as no slip wall on the top and bottom walls. In that case, the velocity of gas at the boundary should be 0. However, when I extracted velocity data from the vtu files, I saw that the velocity was not zero. Here is the code how I extracted the velocity field and the corresponding coordinate.
import vtk
from vtk.util.numpy_support import vtk_to_numpy
import matplotlib.pyplot as plt
def get_coord(file_name):
#read vtu
reader = vtk.vtkXMLUnstructuredGridReader()
reader.SetFileName(file_name)
reader.Update()
output = reader.GetOutput()
data = output.GetCellData()
#convert cell to point
converter = vtk.vtkCellDataToPointData()
converter.ProcessAllArraysOn()
converter.SetInputConnection(reader.GetOutputPort())
converter.Update()
#extract coordinate info
Point_cordinates = reader.GetOutput().GetPoints().GetData()
coord = vtk_to_numpy(Point_cordinates)
return coord
def get_sim_res(file_name):
#read vtu
reader = vtk.vtkXMLUnstructuredGridReader()
reader.SetFileName(file_name)
reader.Update()
output = reader.GetOutput()
data = output.GetCellData()
#convert cell to point
converter = vtk.vtkCellDataToPointData()
converter.ProcessAllArraysOn()
converter.SetInputConnection(reader.GetOutputPort())
converter.Update()
#extract fields info
sim_res = {}
for i in range(data.GetNumberOfArrays()):
field = converter.GetOutput().GetPointData().GetArray(i)
field = vtk_to_numpy(field)
sim_res[data.GetArrayName(i)] = field
return sim_res
file_name = 'BACKGROUND_IC_0010.vtu'
coord = get_coord(file_name)
sim_res = get_sim_res(file_name)
plt.scatter(coord[:,0],coord[:,1],c=sim_res['U_G'])
plt.show()
print(sim_res['U_G'])
As Jeff mentioned, the velocity was actually stored in the center of each face (or edge in this case). In that case, it makes sense that the velocities are not 0 at the boundaries, because the coordinate that I extracted with the get_coord() function above did not really get the coordinate of where the velocity filed was measured. However, if that is the case, it doesn’t seem like I can access the real coordinate of the velocity measurement just using a simple translation, because the coordinate tuple list actually goes from 0 (y_min) to 0.05 m (y_max) on the y direction. So did I misunderstand how the velocity field was stored?
The coordinates of each velocity component is not stored in the vtk file. VTK treats this as u,v,w being located at the cell center where in reality, u is located at the east face center, v at the north face center and w at the bottom face center (here east/west refer to x-direction, North/South refer to y-direction and top/bottom refer to z-direction). Once you apply the vtkCellDataToPointData, you get (u,v,w) values at each cell corner but these are based on cell centered (assumed) values. I am not aware of a vtk filter that can handle staggered velocity components so you may have to do the conversion manually.
In the example above, if all u velocities are stored in the center of the east face, then there should be only 16 data points (the red dots). However, when I used vtk.vtkCellDataToPointData(), it gave me 25 data points of u, meaning each lattice point has a data attached to it. That said, if I do any simple translation, the data will exceed the physical domain of the simulation. I am wondering if I used the wrong data converter?
Thanks for the clarification. What I was trying to say is that the vtk.vtkCellDataToPointData() will not help you because it does not know the velocity components are staggered. You will get values at each cell corner but they are incorrect and doing a translation will not help. To get correct values at the wall, you need to actually know what boundary condition is applied. See the picture below. u is stored at the wall and set to zero so this is good to go (red). To apply a zero v component at the wall (blue), we use a ghost cell and reverse the value of v (from the fluid cell adjacent to the wall) so the interpolated value at the wall is zero. This is what your script needs to do and this is not easy since the BC type is not stored in the vtk file.
Thank you Jeff! I have much better understanding now, but still have some follow up questions. Could you help me explain more? Thanks!
In your sketch, the v vector (blue one going up) is not located at the wall, so it doesn’t need to go to 0, is that correct?
If I want to just get a field value (e.g. v) and exactly where the value was measured (e.g. x and y), what function should I use? It seems vtk.vtkCellDataToPointData() will extract value at each corner, which is not correct for u and v filed. Do you know how to extract information from the center of each face in python? I guess for the (x,y) information I can just do the point data filter, translate them to the center and discard points that are in ghost cells. But I’m not sure how to access the field values at each face center…
Correct, the blue vector is not zero for a no-slip wall. I don’t know of a filter that can do what you want directly. My best guess is you can get the coordinates of each corner and do the calculation manually. Please see How to have access to all cell data? - #18 by jeff.dietiker
I see. Thanks Jeff! I see your point that VTK is just for visualizing, and the actual simulation doesn’t care about coordinates. So I guess there is no automatic way to extract information of field/vector values at each point with corresponding coordinates in any programming language? At last just making sure, the ghost cells are always attached when there is a no-slip wall boundary condition, right? Or is it more complicated than that?
If by automatic, you mean without writing code, then I am not aware of a way to do that. I believe there will be some code writing involved. Now if you start from the vtk file, you don’t have all the information needed. It would be easier to do it directly in MFiX because MFiX knows about fluid and ghost cells. There is always ghost cells next to the BC planes, whether it is a no slip, free slip or inlet/outlet. It gets a lot more complicated with cut cells but I won’t go there at this point.
