Tag: coupling

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.

Fortran 77 essentially does not have global variables for subroutines, instead, common blocks are used to passing data across subroutines when arguments are not preferred. This is conceptually unsafe, though widely used in some legacy F77 codes. The code I’ve been working on adopts this “feature” extensively. The problem is, in Dirichlet-Neumann-like domain decomposition implementation, the solver will be called twice for each iteration, each time with different stiffness matrix setup(due to different subdomains) and boundary conditions. On the other hand, the common blocks make those two calls undistinguishable.

To reuse the legacy solver, my first try was to construct two copies of the solver by putting it into two modules, hoping the explicit interface of the modules would identify the common blocks in different copies. This turns out to be wrong: the subroutines from the two modules still share the same common blocks.

What I finally decided to do is to build a copy of the original solver by renaming all the common blocks. Sed & AWK make this not to difficult. But still I don’t think that’s the most efficient way to tackle this. Consider this is a lesson that not all sequential solvers are easy to parallelize.

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.