FEM, CFD, CAE, and all that.

## Tag: FSI

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

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

### Vortex pair: a mesmerizing photo

Physics Today has this beautiful and mesmerizing photo of vortex pair visualized by laser-induced fluorescence. The description of the photo:

Reproducing the phenomenon at the laboratory scale, Charles Williamson and coworkers at Cornell University have developed a novel laser-induced fluorescence technique for visualizing how vortex pairs interact with the ground. They pool dye on a horizontal surface prior to launching the primary vortices from rotating flaps above it, and then illuminate the system with a laser light sheet. In this image, green dye visualizes the secondary vortices, which have been generated by the ground interaction and swept up above the red primary vortices. The vortex pairs are about 8 cm apart. The lower half of the photo shows the vortex system optically reflected in the surface below.

Charles Williamson’s review of VIV was actually my entry portal to the subject.

### Being irrotational and viscous

I have been working on something involving both potential flows and turbulent flows,  and this is how I found fascinating work by Daniel D. Joseph. In one of his papers, Potential flow of viscous fluids: Historical notes, one can get a peek of the missing link in our textbooks. I have also learned a lot in his web page about the same topic. His other two books I can put my hands on, Stability of Fluid Motions, Elementary Stability and Bifurcation Theory, are the standard textbooks for those topics as well. Shame on me, even though I read those two books, I am not aware of this potential viscous flow, making me wonder how fast I am drifting away from theoretical fluid dynamics with only the computing stuff in my mind.

### ALE formulation

There is a old post for Lagrangian configuration, which is popular in solid mechanics but not in CFD. ALE is dedicated for combining both Eularian and Lagrangian configuration. In Lagrangian, one has the map between Eularian configuration region $\Omega$ and Lagrangian reference configuration domain $\Omega-0$.

$x=\phi(\xi,t), \forall x\in \Omega,\xi\in \Omega_0$

Now assume a new reference region $\Omega_{\ast}$, one have another map $f$:

$x=f(c,t)=\phi(\xi,t),\forall c\in \Omega_{\ast}$

Here $f$ is identity mapping if $\Omega_{\ast}=\Omega$ while is $\phi$ if $\Omega_{\ast}=\Omega_0$. Similar to material velocity in spatial configuration:

$u(x,t)=u(\phi(\xi,t),t)=\phi_{,t}(\xi,t)|_{\xi}$

new reference configuration also has a velocity, called “mesh velocity” in spatial configuration:

$u_{\ast}(x,t)=u_{\ast}(f(c,t),t)=f_{,t}(c,t)|_c$

With proper boundary condition $\phi(\xi,0)=\xi_0$ or $f(c,0)=c_0$ above two equations can be seen as PDE of $\phi$ or $f$ when corresponding velocities are given. Mesh motion $f$ is devised to reduce the mesh distortion during time evolution.

For a function $g(x),\forall x\in \Omega$, there is counterpart of $g$ in other reference configurations:

$g(x,t)=g_0(\xi,t)=g_{\ast}(c,t)$

Take material derivatives:

$\frac{Dg(x,t)}{Dt}|_{\xi}=g_{,t}(x,t)+g_{,i}(x,t)x_{i,t}|_{\xi}=g_{\ast,t}(c,t)+g_{\ast,i}(c,t)c_{i,t}|_{\xi}$

in which $x_{i,t}|_{\xi}\equiv u_i$ is the velocity of material particle $\xi$ in spatial configuration, and $c_{i,t}|_{\xi}\equiv w_i$ is the velocity of material particle $\xi$ in arbitrary reference configuration. On the other hand, fixing material particle marker while taking time derivative of both sides of

$\phi_j(\xi,t)=f_j(c,t)$

gives

$u_j=f_{j,t}|_c+f_{j,i}c_{i,t}|_{\xi}=u_{\ast j}+x_{j,i}c_{i,t}|_{\xi}\Longrightarrow \frac{\partial x_j}{\partial c_i}w_i=u_j-u_{\ast j}$

where $u_{\ast}$ is the velocity of fixed mesh point $c$ in spatial configuration.

Combining last two conclusion we have

$Dg(x,t)/Dt=\dot{g_0}=g_{\ast,t}|_c+g_{,j}(u_j-u_{\ast j})$

Define $b\equiv u-u_{\ast}$ as the difference between material velocity and mesh velocity in spatial configuration, when both material particle marker and mesh point marker are at the same spatial position at that moment. Finally, the transformation of $g$ and $g_{\ast}$ is

$Dg/Dt|_{\xi}=g_{\ast,t}|_c+b\cdot \nabla_x g$

All the conservation laws in ALE formulation can be derived from expression above.