FEM, CFD, CAE, and all that.

## Tag: FEM

### Soft shell structure under hydrodynamical load

In the middle of constructing a reduced-order model (ROM) of a highly flexible structure interacting with surround flow, I built a structure model and loaded it with oscillating pressure at the lower part. The animation from the simulation is as follows.

There are several challenges to this type of FEA problems, mostly are related to the highly flexibility nature of the material. The structure is first subject to a quasi-static hydrostatic load, which makes it buckle to a new equilibrium position. After that, ambient oscillating pressure with a much smaller magnitude than the initial buckling load is applied to make the fabric vibrate.

This is a perfect example for using explicit-implicit switching for time integration. Since it’s usually difficult for an implicit solver to find the bulking configuration, as a result of the fact that the Newton-Ralphson method does not gaarantee  convergence, one needs to issue an explicit analysis for the initial loading stage. Some solvers provide specific “bulking simulation” controls for this type of the problem. After the initial loading the structure to a new equilibrium (passing a bifurcation point, as one clearly see in the animation), we proceed to impose the oscillating load. This type of load is supposed to be of a much smaller magnitude than the initial load, so that the structural vibration would be a perturbation near the new equilibrium.

Since explicit method suffers from a general stability problem (notorious forward Euler blow-up in every toy example in every textbook), it is only applicable to short-duration problems such as impact/contact, and it is very easy for the blow-up to happen during a cyclic loading problem, due to change of sign of the derivatives. Thus we need to turn off the explicit simulation and replace it with an implicit one. The benefit? (Unconditional) stability and hence greater allowable time step. The cost, on the other hand, is much higher since now we are actually inverting the stiffness matrix and perform in-time iterations.

### 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    ##
## 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”.

### Domain Decomposition for Poisson equation in 1D

It’s not really Poisson equation, just a second order ODE in that sense. But it serves a good example to demonstrate the idea of domain decomposition.

Using Dirichlet-Neumann method, one of the two subdomains gives out Neumann boundary condition and receives Dirichlet boundary data at the interface, and the other one does the opposite. The Mathematica modul would be somthing like this:

DNDecomp:   Main module;
DBCLaplace: linear FEM module for Dirichlet data ua at xa, and ub at xb;
MBCLaplace: linear FEM module for Mixed data ua at xa, and flux db at xb;
(xi,ui):    interface position and data;
theta:      relaxation coefficient to update ui in each iteration.
DNDecomp[xa_, ua_, xb_, ub_, xi_, ui0_, n1_, n2_, iter_, \[Theta]_] :=
Module[{ui, di, x1, u1, x2, u2, h1, h2, sol},
f[t_] := 2;
ui = ui0;
h1 = (xi - xa)/n1; h2 = (xb - xi)/n2;
sol = {};
Do[{
{x2, u2} = DBCLaplace[xi, ui, xb, ub, n2];
di = (u2[[2]] - u2[[1]])/h2;
{x1, u1} = MBCLaplace[xa, ua, xi, di, n1];
ui = \[Theta] u1[[-1]] + (1 - \[Theta]) u2[[1]];
AppendTo[sol, {{x1, u1}\[Transpose], {x2, u2}\[Transpose]}]},
{iter}
];
sol
]

In above example force function $f=2$ gives a parabolic exact solution. The effect of the relaxation coefficient $\theta$ is obvious in plots below.

Above: theta=0.7, below: theta=0.5

Below is another example with exact solution $u(x)=-x \sin(x)$.

Iteration = 6, theta = 0.5

### Three faces of FEM

Finite Element Method (FEM) acquires several interpretations through its developments over decades. Assuming we are solving a PDE of form

$L(u)=0$

where $L$ is partial-differential operator, $u=u(x)$ is unknown function. FEM could be interpreted in the following ways (let me know if you have extra answers):

• Rayleigh-Ritz
• Galerkin method
• Method of Weighted Residuals (MWR)

The first one, Rayleigh-Ritz method, is actually how structure engineers were inspired to devise FEM. The connection of this method and FEM is the two way street of variational formulation of PDE, which gives the Euler equation of certain functional over $u(x)$. The derivation is done through searching for the minimum of the functional (usually some form of energy ). Rayleigh-Ritz method addresses the PDE problem by assuming the form of $u$ as

$u(x)=\sum_i^N u_i \Phi_i(x)$

then bringing this formulation into PDE’s functional formulation, in which $u_i$‘s are acquired by solving a minimization problem. Original Rayleigh-Ritz method adopts functions with global support as $\Phi_i(x)$. A textbook example is the boundary value problem of a beam’s modes, in which we can assume $\Phi(x)$ as sinusoidal functions. $\Phi_i(x)$ are called trial functions.

In order to have a PDE’s variational formulation in a general way, and in the same time remove some regularity requirement of the solution, Galerkin method is used. The so called test function are the ones with expected regularity (usually infinitive). By the procedure of Galerkin method, test functions $\Psi_i(x)$ are applied to $L(u)=0$ and the integral formulation is acquired. If $L$ is linear, this integral formulation is the same as functional formulation in minimization problem by Rayleigh-Ritz method. Here we have introduced two sets of functions, trial functions which define the finite dimensional approximation of $u(x)$, and test functions which defines how close (accurate) the PDE is solved. The name of “test” can be understood in the following way. How close a variable $v$ is close to zero could be measured by test it with other variables:

$(v,w)=0$

where $w$ is known variable used for testing, and $( , )$ is proper inner product. Conceptually, the greater the variable domain consisting of $w$, the closer to zero $v$ is. And if $v$ can make every $w$ satisfy above equation, it should be zero itself by definition.

So we can see, the space of test function $\{\Psi_i\}$ defines the extension, and in turn, accuracy, the PDE $L(u)=0$ stands by having

$(\Psi_i, L(u))=0$

A subcategory of Galerkin method is Bubnov-Galerkin method, where test functions and trial functions are the same. Otherwise, it’s called Petrov-Galerkin method. Historically, trial function is also commonly referred as shape function.

Another abstraction of FEM is by MWR. Since all the numerical methods for solving PDE relies on moving from infinite dimensional solution spaces to finite ones, errors are introduced in all the methods. So instead of having $L(u)$ equal to exact zero, $u(x)$‘s numerical approximation $\tilde{u}(x)$ introduces residual $R(x)$:

$L(\tilde{u}(x))=R(x)$

If numerical solution is in the form of trial function’s expansion, the coefficients (coordinates) $u_i$ could be acquired by restricting $R$. In the light of above philosophy, $R$ could be restricted as to be zero in some sense, specifically,

$(R(x),\Psi_i(x))=0$

for some $\Psi_i$ and $(\cdot,\cdot)$. This can be looked as putting weight on some locations defined by $\Psi_i$: the greater $\Psi_i$ is on certain location, the more restricted $R$ is there. By this $\Psi_i$‘s are also called weighted functions.

Though mathematically Galerkin method and MWR are in the same look, they are actually on different emphasis. In the former, we are specifically looking for integral form, with shifting regularity to test functions in mind, while in MWR we are focusing on minimizing residuals, and following integral form is just the math “trick” played on certain weight functions. In fact, by using different weight function in MWR, we can recover other numerical discretizations. Besides Bubnov/Petrov Galerkin methods, $\Psi_i=\delta(x-x_i)$ gives collocation method, $\Psi_i=1$ within certain cell and $0$ otherwise gives FVM, and $\Psi_i=\partial R/\partial u_i$ gives least squares method.