Mesh Mess

FEM, CFD, CAE, and all that.

Tag: CFD

A water tank mesh generater in R

I have been using LS-DYNA ICFD solver for several years. Basically, it’s a incompressible CFD solver addon over the platform of LS-DYNA, which has its original flow solver based on compressible flows. It’s nice to work with people in LSTC, the develop company of LS-DYNA, but their pre-post processor really doesn’t work well with me. For my own research, during the development of a numerical wave tank, most of the time what I need is just a tank mesh. In ICFD, mesh generation is done automatically within the flow domain, so the mesh information that is requried at the file input is just the meshed surface: the walls, the free surface, the bottom, and the other  domain boundaries. Since open lsprepost to build a mesh file from scratch is the last thing I’d like to do when every time certain mesh parameters are changed, such as in convergence test, I just write a script to do that work. The output plots for mesh example from this R script are like this:

The first plot shows all the faces of generated. Plot 2 shows the mesh for bottom, free surface and the “ceiling”, together with “side walls”. Plot 3 shows free surface mesh and two ends of the tank, where the wave are generated and dissipated.

R turns out to be handy for this job, for its way of handling data, especially the recycling and the mapping due to its embracing of functional paradigm. As the script below shows. The code may be fragile, still saves a lot of time from GUI.

##########################################################################
## Copyright (C) 2012-2013 Yi Zhang                                     ##
##                                                                      ##
## Author: Yi Zhang                                                     ##
## Created: 30 Nov 2012                                                 ##
##                                                                      ##
## Purpose: This file generates a surface mesh as used by ICFD          ##
## in LS-DYNA keyfile format.                                           ##
##                                                                      ##
## This file is free software; you can redistribute it and/or modify    ##
## it under the terms of the GNU General Public License as published by ##
## the Free Software Foundation; either version 2, or (at your option)  ##
## any later version.                                                   ##
##########################################################################

rm(list=ls())
meshfile <- "~/tmp/tank3d_mesh.k"

d <- c(10., 1., 0.5, 0.3) # dim: (x, y, z, z-air)
m <- c(75,14,10,10)           # nb. of elements in (x,y,z,z1)
n <- m + 1                # nb. of nodes
h <- d/m

## total nb. of elems. & nodes.
mtol <- 3*m[1]*m[2] + 2*(m[2]*(m[3]+m[4]) + m[1]*(m[3]+m[4]))
ntol <- 3*n[1]*n[2] + 2*(n[1]*(m[3]+m[4]-2)++(m[2]-1)*(m[3]+m[4]-2))

## beginning and the end id's of the each face
facegen <- function (members) {
  faces <- matrix(c(1, members[1]), nrow=1, ncol=2)
  for (i in 2:length(members)) {
    faces <- rbind(faces, c(tail(faces,1)[2]+1,
                            tail(faces,1)[2]+members[i]))}
  return(faces)}
edgegen <- function (i, k) {
  ## j = (1, 2, 3, 4)  (S, E, N, W)
  edge <- list(nface[i,1] : (nface[i,1]+m[k]),
               seq(nface[i,1]+m[k], nface[i,2], n[k]),
               (nface[i,2]-m[k]) : nface[i,2],
               seq(nface[i,1], nface[i,2], n[k]))
  return(edge) }
nbn   <- c(n[1]*n[2], n[1]*n[2], n[1]*n[2], n[1]*n[3],
           n[1]*n[3], n[1]*n[4], n[1]*n[4], n[2]*n[3],
           n[2]*n[3], n[2]*n[4], n[2]*n[4]) 
nbm   <- c(m[1]*m[2], m[1]*m[2], m[1]*m[2], m[1]*m[3],
           m[1]*m[3], m[1]*m[4], m[1]*m[4], m[2]*m[3],
           m[2]*m[3], m[2]*m[4], m[2]*m[4]) 
nface <- facegen(nbn)
mface <- facegen(nbm)
nedge <- list(edgegen(1, 1), edgegen(2, 1), edgegen(3, 1),
              edgegen(4, 1), edgegen(5, 1), edgegen(6, 1),
              edgegen(7, 1), edgegen(8, 2), edgegen(9, 2),
              edgegen(10, 2), edgegen(11, 2))

## nodes
x   <- rep(seq(0, d[1], h[1]), n[2])
y   <- as.vector(sapply(seq(-d[2]/2, d[2]/2, h[2]), function (x) rep(x, n[1])))
nod <- rbind(cbind(x, y, -d[3]),
             cbind(x, y, 0.   ),
             cbind(x, y,  d[4]))

x   <- rep(seq(0, d[1], h[1]), n[3])
z   <- as.vector(sapply(seq(-d[3], 0, h[3]), function (x) rep(x, n[1])))
nod <- rbind(nod,
             cbind(x, -d[2]/2, z),
             cbind(x,  d[2]/2, z))

x   <- rep(seq(0, d[1], h[1]), n[4])
z   <- as.vector(sapply(seq(0, d[4], h[4]), function (x) rep(x, n[1])))
nod <- rbind(nod,
             cbind(x, -d[2]/2, z),
             cbind(x,  d[2]/2, z))

y   <- rep(seq(-d[2]/2, d[2]/2, h[2]), n[3])
z   <- as.vector(sapply(seq(-d[3], 0, h[3]), function (x) rep(x, n[2])))
nod <- rbind(nod,
             cbind(0.,   y, z),
             cbind(d[1], y, z))

y   <- rep(seq(-d[2]/2, d[2]/2, h[2]), n[4])
z   <- as.vector(sapply(seq(0, d[4], h[4]), function (x) rep(x, n[2])))
nod <- rbind(nod,
             cbind(0.,   y, z),
             cbind(d[1], y, z))

nod <- cbind(1:(dim(nod)[1]), nod)

## elements
f <- function (i, j) {
  p <- nface[i,1] : nface[i,2]
  a <- p[(p-nface[i,1]+1) %% n[j] != 0];
  b <- p[(p-nface[i,1]+1) %% n[j] != 1]
  return(cbind(a[a <  nface[i,2] - n[j]],
               b[b   nface[i,1] + n[j]],
               a[a >  nface[i,1] + n[j] - 1]))
}
ele <- cbind(1:mtol,
             unlist(sapply(1:11, function (i) rep(i,
                                                  mface[i,2]-mface[i,1]+1))), 
             rbind(f(1, 1), f(2, 1), f(3, 1), f(4, 1),
                   f(5, 1), f(6, 1), f(7, 1),
                   f(8, 2), f(9, 2), f(10, 2), f(11, 2)))

## Eliminate redundant nodes (at corners/edges)
e <- ele[,3:6]
nr <- function (ma, a, b) {
  for(i in 1:length(a)) ma[ma == a[i]] <- b[i]
  return(ma) }
a <- rbind(c(11, 2, 7, 2), c(11, 4, 6, 2),
           c(10, 2, 7, 4), c(10, 4, 6, 4),
           c(9 , 2, 5, 2), c(9 , 4, 4, 2),
           c(8 , 2, 5, 4), c(8 , 4, 4, 4),
           c(11, 1, 2, 2), c(11, 3, 3, 2),
           c(10, 1, 2, 4), c(10, 3, 3, 4),
           c(9 , 1, 1, 2), c(9 , 3, 2, 2),
           c(8 , 1, 1, 4), c(8 , 3, 2, 4),
           c(7 , 1, 2, 3), c(7 , 3, 3, 3),
           c(6 , 1, 2, 1), c(6 , 3, 3, 1),
           c(5 , 1, 1, 3), c(5 , 3, 2, 3),
           c(4 , 1, 1, 1), c(4 , 3, 2, 1))
for(i in 1:dim(a)[1]){
  f1 <- a[i,1]; s1 <- a[i,2]; f2 <- a[i,3]; s2 <- a[i,4]
  e <- nr(e, unlist(nedge[[f1]][s1]), unlist(nedge[[f2]][s2]))}

nid <- sort(unique(as.vector(e)))
nod <- nod[nid,]
nod[,1] <- 1:ntol
e <- nr(e, nid, 1:ntol)
ele[,3:6] <- e

## Write file in LS-DYNA keyfile format
write("*KEYWORD", file=meshfile) 
write("*MESH_SURFACE_ELEMENT", file=meshfile, append=TRUE) 
write.table(cbind(format(ele[,1], width=8,justify="right"),
                  format(ele[,2:6],
                         width=7,justify="right")),
            file=meshfile, append=TRUE, quote=FALSE,
            row.names=FALSE, col.names=FALSE)  
write("*MESH_SURFACE_NODE", file=meshfile, append=TRUE) 
write.table(cbind(format(nod[,1],
                         width=8, digits=1, justify="right"),
                  format(nod[,2:4], width=15, digits=7,
                         justify="right")),
            file=meshfile, append=TRUE, quote=FALSE,
            row.names=FALSE, col.names=FALSE)  
write("*END", file=meshfile, append=TRUE)

You get what you mesh for

In finite element modeling, in fact, in all the numerical modeling, the resolution of the results depends on the space where the problem is projected. As to PDE’s numerical solution, the projection is from space with infinite dimension to space with finite dimension. So no matter it’s finite element or finite volume or finite difference, the size of the mesh is what you can see at the finest scale. Beyond that, it’s either physical modeling (such as turbulence subscale modeling) or some other educated guess.

Here is an example.

Water surface as wave coming in

This plot shows a finite element solution of Navier-Stokes equation, in a water waves problem. The free surface is tracked using Level Set method. Taking closer look reveals the following not-so-smooth “water surface”:

Free surface, enlarged

The resolution of the free surface is restricted by the size of the element. For any mesh, examing close enough reveals similar wiggles. To put it another way, in the numerical simulation, free surface is not really a surface, it’s more like a layer with finite thickness. If I plot the centroid of the elements that reside on the free surface, it looks like this:

Elements on the free surface

This is what a non-polished numerical result should look like. Smoothing such as with cubic splines sweepts everything under the table, including the message: “you get what you mesh for”.

Guidelines for CFD papers

An interesting editorial from Journal of Fluids Engineering by Malcolm J. Andrews, essentially about is worth of showing to people in numerical modeling of flows, contains some guidelines for submitting numerical results to JFE, such as:

Providing complex figures may be nice for a presentation to an audience or sponsor, but are often not scientific or quantitative unless great care is taken in their presentation/discussion or they illustrate a novel aspect of the work, and rarely provide much detailed insight into the problem under consideration. The usual x-y plots are often more illuminating, but also require the author to spend more time thinking about what is important and why which helps make the article more archival.

Simply reporting one parameter when, in fact, there are multiple parameters that self-interact suggests that the author does not understand the diagnostic, or its proper use, and also the basic elements of the flow itself e.g., simply reporting pressure, and not associated velocity fields, might indicate a lack of basic understanding. The article should include a detailed description of the results, their consequences, and their importance i.e., simply stating values or shapes does not warrant archival.

Nondimensional parameters serve not only to collapse data but they demonstrate an understanding of the basic parameters that control the processes of interest and form the basis of generality that can underlie resulting archival value formulas. Not expressing results in nondimensional form substantially weakens the archival value, suggests that the author does not understand the fundamental flow physics, and also suggests that the results have no generality or archival value.

It is crucial to provide the applicable parameter ranges for the commercial software or diagnostic and ensure that they are met in the current application e.g., this might mean answering the question about an appropriate use of a turbulence model, the Reynolds number range of the experiment, or the Stokes relaxation time of particle in the flow relative to the time scale of interest in the flow.

Log law at a smooth plate

Log-law is THE most important and elegant analytical result in turbulence, considering many others are dauntingly complicated,  though justified. What log law says is  the fact that for a turbulent flow over flat plate, there is this section of flow away from the wall, within the turbulent boundary layer,  in which the relation between the distance and mean flow velocity is simply of log function. This conclusion could be drawn by an asymptotic analysis matching outer flow and inner flow.   Any CFD trying tempted to address turbulence models should recover the “law of the wall“, among which log law is a part of, before moving forward to more advanced problems.

The delicacy of the modeling in this case is the boundary condition (shocking!) at the wall, because log law is only valid when y+ is approximately above 30 ( the other cap depends on Reynolds number ), and the boundary conditions there for unknown variables in turbulence models are generally not known a priori, where is all kinds of wall function come into play.

Channel Flow

Plot of CFD result against log law

Across an interface

I was reading this interesting paper by Idelsohn et al when I noticed that for the Navier-Stokes equation the velocity boundary condition in normal and tangential direction have different physical meanings. We know that the NSE requires complete (against normal) velocity profile to determine the problem. Through a point in incompressible flow domain, we draw an arbitrary straight (to eliminate tension’s effect) line as interface separating two sides of the domain. In order to determine the problem on side A, velocity at interface should be applied as boundary condition. Now, normal component of  velocity at that point should be continuous across the interface (since it’s arbitrary instead of physical), and it corresponds to the incompressibility condition: what goes in should equal to what goes out. It is easy to see when we imagine a slice of spacial volume (against material volume) along the interface, and its thickness is much less than its span. The match of normal velocity indicates the conservation of mass within that volume.

On the other  hand, the match of tangential component of the velocity means something else. Imagine that they does not match across the interface, what would happen? Well, the gradient of tangential velocity would be infinite (remember that interface has no thickness, so the dx in the derivative goes to zero), and the shear stress in that direction would be infinite (Newtonian fluid), and apparently it’s physically impossible. This is like a boundary layer with zero thickness. For a boundary layer, viscous effect dominates and causes relatively large viscous force. In our imagined case this means the infinite viscous force and unbalance the momentum conservation.

So here are two conservation laws hidden in one boundary condition. Mathematically it suffices to see this when we degenerate from NSE to Euler equation. By losing one order of the PDE, the equation can not meet both normal and tangential boundary condition, and this corresponds to the slip condition when only normal velocity is prescribed (usually zero) out side the boundary layer.