FEM, CFD, CAE, and all that.

## Category: Engineering

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

### Moving to ggplot2 from gnuplot

Finally decided to move to ggplot2, after years’ use of gnuplot. I love gnuplot’s simple and straightforward, though its syntax could be a pain sometimes. All the plots in my previous latex docs are generated from gnuplot, using its nice latex terminal for publication (this is still unmatched by ggplot2, about this in a minute).

One thing that motivates me to use ggplot is Hadley’s book. The idea of grammar of graphics looks very sweet. Another reason is that by moving in R ecosystem I can use sweave to produce dynamic docs for some previosu docs (new docs will be written in org-mode so that’s not necessary). As to this,  another option is to use org-mode (babel) + matplotlib + latex. I love the elegant idea of writing whole thing in org,  and there are plenty examples to check out, see here and here. On the other hand, sweave is just that simple but has been having some long-time problems (some are addressed by knitr). As to matplotlib, it does the work alright, but in my opinion the plots are not as nice looking as those from ggplot.

However, gnuplot and matplotlib have one great advantage against ggplot, which is latex text rendering in plots. Latex rendering is supported by activating the latex option in matplotlib, similar to using latex terminal in gnuplot. While in ggplot, there is no such thing out of box. Workaround includes using tikzdevice (seems the project is in stall) or extrafont, but neither is satisfactory to me. That’s why org-mode gives a straightforward (though not satisfactory) solution by allowing me use gnuplot when needed, at least for a while, because better solution appears.

### Prohibit the metric system!

Through Irreal I found this eloquent argument against metric system in favor of US system ( I don’t call it imperial since UK has officially adopted metric). 8 signatures have been gathered toward 25000 goal, astonishingly less than the ratio of nations endorsing the almighty US system (which is just US, essentially).

### Check Fortran 77 line length

Another issue with F77 is its length limit: each line should be less than 72 char, otherwise the rest of the content will be ignored. Of course I can demand compiler work in free-form mode, but that would cause other compatibility problems. This happens occasionally after I’ve done some “sed s/” stuff over the variable names. So I need to remember checking line length, by adding an awk line in the Makefile:

awk '/^[^cC]/ && length > 72 {print FILENAME,FNR}' *.f


### Harmonic vibration of a damped system

Unlike in previous post that phase change could only be $0$ or $\pi$, in the damped system, as the ratio between dynamics response and static response is expressed in the same way as the undamped case
$\frac{u_d}{u_{s0}}=R_d\sin(\omega t-\phi),$
the phase lag $\phi$ varies in $[0,\pi]$. The following Mathematica outputs the displacement response factor, phase lag and normalized time history as the function of damping ratio $\zeta$ and frequency ratio $w=\omega/\omega_0$.

Rd[w_, \[Zeta]_] := 1/Sqrt[(1 - w^2)^2 + (2 \[Zeta] w)^2]
\[Phi][w_, \[Zeta]_] := -ArcTan[1 - w^2, -2 \[Zeta] w]
ud[w_, \[Zeta]_, t_] :=
Rd[w, \[Zeta]] Sin[2 Pi (t - \[Phi][w, \[Zeta]]/(2 Pi))]
us[t_] := Sin[2 Pi t]
Manipulate[
GraphicsRow[{GraphicsColumn[{Plot[Rd[w, \[Zeta]], {w, 0, 3},
AxesOrigin -> {1, 1}, PlotRange -> {{0, 3}, {0, 5}},
AxesLabel -> {"\!$$\*FractionBox[\"\[Omega]\", SubscriptBox[\"\ \[Omega]\", \"0\"]]$$", "\!$$\*SubscriptBox[\"R\", \"d\"]$$"},
Epilog -> {PointSize[Large], Red, Point[{w, Rd[w, \[Zeta]]}]},
Ticks -> {{0, 1, 2, 3}, {0, 1, 2, 3, 4, 5}}],
Plot[\[Phi][w, \[Zeta]], {w, 0, 3}, AxesOrigin -> {1, Pi/2},
PlotRange -> {{0, 3}, {0, Pi}},
AxesLabel -> {"\!$$\*FractionBox[\"\[Omega]\", SubscriptBox[\"\ \[Omega]\", \"0\"]]$$", "\[Phi]"},
Epilog -> {PointSize[Large], Red,
Point[{w, \[Phi][w, \[Zeta]]}]},
Ticks -> {{0, 1, 2, 3}, {0, Pi/2, Pi}}]}, Frame -> All],
Plot[{ud[w, \[Zeta], t], us[t]}, {t, 0, 2},
PlotStyle -> {Thick, Dashed},
AxesLabel -> {"t/T",
"\!$$\*SubscriptBox[\"u\", \"d\"]$$/\!$$\*SubscriptBox[\"u\", \ \"st\"]$$"}]}, Frame -> All], {{w, 0.5, "Frequency ratio"}, 0,
3}, {{\[Zeta], 0.2, "Damping ratio"}, 0, 1}]

Special cases include

1. Long period excitation, i.e., $\omega\ll\omega_0$, gives a pseudo-static response. In this case, the system “waits” until it “feels” the excitation completely. $R_d$ is greater but very close to $1$, and the displacement is essentially in phase with excitation force, in other words, dynamic effect is near to none.

2. Short period excitation gives very small $R_d$, though the phase lag is $\pi$. Here the system barely reacts when the load is reversed, thus leads to small displacement.

3. Resonant period, i.e., $\omega\approx\omega_0$, leads to $\phi\approx\pi/2$. Now $R_d$ is very sensitive to damping change, namely, the response is controlled by the damping: a small change of damping ratio $\zeta$ leads to great reaction of the structure. When $\omega=\omega_0$, we have $R_d=1/(2 \zeta)$.