FEM, CFD, CAE, and all that.

## Category: Numerical

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

### POD snippet in R

I’ve been working on proper orthogonal decomposition (POD) for a while. Though it’s nice to have packages such as modred out there, I cannot help myself writing a snippet using R, since lots of my previosu data processing was in that language. Without all the theory chitchat  like K-L transformation, from the the implementation-on-discrete-data perspective the (snapshot) POD method boils down to just singular value decomposition (SVD), and it’s out of box in R, which makes the task mostly straightforward. That’s why a POD function in R essentially calls svd in R (and it further calls LAPACK routines `DGESDD` and `ZGESVD`.):

```pod <- function(arr, nmod=dim(arr)[1], nv=nmod){
pod <- svd(arr, nmod, nv)
return(pod)
}```

Here “nmod” is the number of POD modes to be retrieved. Though it’s default to be the length of each snapshot, it’s rarely necessary this large. Mostly of the time a reasonable correlated set of snapshots requires less than 10 modes to obtain most of the system energy.

With pod of data “arr” obtained, we can do several thing immediately. We  can obtain the modes by retrieving the right singular vector:

```pod_u<- function(pod, nnod){
u <- matrix(unlist(pod["u"]), nrow = nnod )
return(u)
}```

or get the singular values:

```pod_sv<- function(pod){
sv <- as.vector(unlist(pod["d"]))
return(sv)
}```

or the coefficients of the pod modes for the approximation:

```pod_coef<- function(pod, nshot, nmod){
v <- matrix(unlist(pod["v"]), nrow = nshot )
sv <- pod_sv(pod)[1:nmod]
coef <- diag(sv) %*% t(v)
return(coef)
}```

and eventually the approximation itself

```pod_approx<- function(pod, nshot, nmod, nnod){
u <- pod_u(pod, nnod)
coef <- pod_coef(pod, nshot, nmod)
approx <- u %*% coef
return(approx)
}```

### 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)```

### Store a polynomial in non-symbolic language

Sometimes one wants to store a polynomial in a non-symbolic language, such as Fotran or C. A easy way to achieve this is to take advantage of the isomorphism between polynomial space $P(n)$ and space $R^{n+1}$. Namely, we can build an array with entries as the coefficients of the polynomial to be stored. Furthermore, considering these coefficients are just the coordinates of the polynomial in the space spanned by the canonical basis $1$,$x$,$x^2$,$\dots$, we can express the polynomial using any basis expansion.

### Legendre polynomial generation in Fortran

In spectral-based numerical method, orthogonal polynomial basis is commonly used to generate “nice” linear systems. A popular candidate of such polynomial family is Legendre polynomial.  As other Jacobi polynomials,  recursive definition makes it easy to generate a high order member.  The following is  a Fortran snippet to generate Legendre polynomial, following Bonnet’s recursion formula.

```! Copyright (C) 2012 Yi Zhang

! Author: Yi Zhang <zhangy2@onid.orst.edu>
! Keywords: numerical, tools, polynomial

! This file is free software; you can redistribute it and/or modify
! it under the terms of the FreeBSD License. So essentially do
! whatever you want.

module legendre_poly

type polycoef
real, pointer :: p(:)
end type polycoef

contains

subroutine gen_legendre_poly()
! q_n = p_(n-1) = 1/(n-1)*[(2n-3)xp_(n-2)-(n-2)p_(n-3)]
!     = 1/(n-1)*[(2n-3)xq_(n-1)-(n-2)q_(n-2)]
implicit none
integer, parameter :: n = 11
type(polycoef) :: a(n)
integer :: i, j

DO i = 1, n
allocate(a(i) % p(i))
a(i) % p = 0
ENDDO
a(1) % p(1) = 1
a(2) % p    = (/0,1/)
do i = 3, n
do j = 1, i-2
a(i) % p(j) = -(i-2) * a(i-2) % p(j)
enddo
do j = 1, i-1
a(i) % p(j+1) = a(i) % p(j+1) + (2*i-3) * a(i-1) % p(j)
enddo
a(i) % p = a(i) % p / (i-1)
write(*,*) a(i) % p
enddo
end subroutine gen_legendre_poly

end module legendre_poly```

DG methods usually can use up to order 10 for high order scheme. So this snippet can be used to generate polynomial for such applications’ basis function, usually projected to a canonical element.