Enable quasi-static analysis

[mrp089]

Part of #32.

I'm trying to figure out how I can achieve a quasi-static analysis (i.e. no inertia/viscosity effects but time-dependent load) within the svFSI time integration. However, I'm unsure how to modify the gen-alpha time integration for this purpose.

For a quasi-static analysis, I think I would have to:

1. set mass = damping = 0 to ignore inertia effects
2. set velocity = acceleration = 0 to enforce a static analysis
3. set alpha_f = alpha_m = 0 to eliminate any dependence on the previous time step
   (beta, gamma, Delta t should now be irrelevant)

So far, I have done for each step:

1. set Density: 0.0 in the input file
2. set An, Ao, Ag, Yn, Yo, Yg to zero before and after GLOBALEQASSEM in MAIN.f

For step 3, I'm running into problems.

From my nonlinear FEM text book, this is what I expect the dynamic effective tangential stiffness matrix (lK in svFSI) to look like:

This is what happens in svFSI to calculate lK:

• calculate factor amd = am * rho + af * gamma * Delta t * damping
  

  svFSI/Code/Source/svFSI/STRUCT.f

  Line 170 in e258b9f

  amd = eq(cEq)%am*rho + eq(cEq)%af*eq(cEq)%gam*dt*dmp

• calculate factor afl = af * beta * Delta t**2
  

  svFSI/Code/Source/svFSI/STRUCT.f

  Line 171 in e258b9f

  afl = eq(cEq)%af*eq(cEq)%beta*dt*dt

• add mass matrix and geometrical stiffness matrix:
  

  svFSI/Code/Source/svFSI/STRUCT.f

  Line 270 in e258b9f

  T1 = amd*N(a)*N(b) + afl*NxSNx

• add material stiffness matrix (do for all components):
  

  svFSI/Code/Source/svFSI/STRUCT.f

  Line 278 in e258b9f

  lK(1,a,b) = lK(1,a,b) + w*(T1 + afl*BmDBm)

This looks different from the textbook. For one, it looks like everything has been multiplied by beta * Delta t**2. Also, there's af instead of 1-af and am instead of 1-am. Here, setting alpha_f = af = 0 results in afl = 0 and thus lK = 0, which makes the system singular.

@vvedula22 Do you have any idea where I went wrong? I don't think I fully understand the svFSI implementation of gen alpha. Also, I couldn't find the part in the code where velocities and accelerations are updated.

[mrp089]

I found the update of velocities and accelerations in PIC.f:

svFSI/Code/Source/svFSI/PIC.f

Lines 137 to 145 in e258b9f

	coef(1) = 1._RKIND - eq(i)%am
	coef(2) = eq(i)%am
	coef(3) = 1._RKIND - eq(i)%af
	coef(4) = eq(i)%af
	DO a=1, tnNo
	Ag(s:e,a) = Ao(s:e,a)*coef(1) + An(s:e,a)*coef(2)
	Yg(s:e,a) = Yo(s:e,a)*coef(3) + Yn(s:e,a)*coef(4)
	Dg(s:e,a) = Do(s:e,a)*coef(3) + Dn(s:e,a)*coef(4)

The parameters are initialized here:

svFSI/Code/Source/svFSI/INITIALIZE.f

Line 79 in e258b9f

am = (2._RKIND - eq(iEq)%roInf)/(1._RKIND + eq(iEq)%roInf)

svFSI/Code/Source/svFSI/INITIALIZE.f

Line 134 in e258b9f

eq(iEq)%af = 1._RKIND/(1._RKIND + eq(iEq)%roInf)

Looks like this is complementary to the definitions in my textbook:

I would thus require am=af=1 for a quasi-static analysis.

[vvedula22]

For quasi-static analysis, all you have to do is remove the contribution of inertia by setting Density: 0.0. I'd not recommend changing the code to zero out acceleration (An, Ao) and velocity (Yn, Yo) variables. Even if you were to leave them as is, they would not have any contribution to the residue if you set density to 0.

It looks like you already identified portions of the code where the generalized-alpha parameters are set and computed. Note that the definitions we use are somewhat different from the original paper by Chung and Hulbert. Our implementation follows the definitions for alpha_m, alpha_f, gamma, and beta from Bazilevs et al.. Also note that we are solving for increments in acceleration (\ddot{\Delta u}) and not for increments in displacement (\Delta u).

I am not sure why you'd need to set am=af=1 for a quasi-static analysis. Such a setting would make the scheme to be backward Euler and would lose second-order time accuracy. This is not required. Simply use \rho_\infty based one-parameter definitions for evaluating alpha_m and alpha_f that guarantee unconditional stability and second-order accuracy. Setting rho_infty to be 0.5 will give you optimal damping.

[mrp089]

Thanks for the quick reply @vvedula22 and the reference!

I think we have a different understanding of the quasi-static analysis. In my understanding, we don't want any time integration, but n individual solutions for n time steps. The only place where we'd use a previous solution would be as an initial guess for the new solution. So I think we don't want any damping.

I wasn't aware that we solve for increments in acceleration. (To me, it seemed like we're solving for beta * Delta t**2-times the displacements.) In that case, is it even possible to achieve a truly static solution? For a quasi-static solution the increments in acceleration should always remain zero, right?

[vvedula22]

@mrp089 You're right! we don't do a quasi-static analysis in svFSI but pose the problem in a dynamic setting and remove the contribution from the inertial force. This leads to a time integration error proportional to beta * Delta t**2 as we solve for the increments in acceleration and not for increments in the displacement. The code was originally developed for FSI problems that involve inertia and so we never really had to deal with this. We could get away by setting Density: 0.0 for testing quasi-static problems. A truly static solution may not be possible for large deformation problems as the loading is usually done in small increments.

However, a pure quasi-static problem could be done in svFSI with some minor modifications to the code. You may define a global variable, say quasiStatic in COMMOD, and read a user input Quasi-static problem: t only for struct and shell equations. If set to true, you may set rho_inf = 0 resulting in alpha_f = 1. You may also set afl = 1 and amd = 0 in STRUCT.f and only update Dn directly with the increments stored in R in PICC subroutine. If you would like, I can set this special case in a day or two.

[mrp089]

Thanks, again, for your help @vvedula22! I think we're on the same page now. What I was wondering:

We multiply the stiffness matrix contributions by beta * Delta t**2 and then also multiply the update by beta * Delta t**2. Thus, these two factors should cancel out. However, when testing my quasi-static analysis with a simple tube inflation example, I'm still seeing a dependence of the results on the time step Delta t. You also mentioned the time integration error being proportional to beta * Delta t**2. How is that? What part of the time integration am I missing?

[vvedula22]

@mrp089 Sorry for the delay.

I was referring to the final update step that involves multiplication with \beta \Delta t^2 as the time integration error as the increments are on acceleration.

The first \beta \Delta t^2 used in assembling the stiffness matrix results from applying the chain rule to compute the derivative of the residue with respect to the acceleration.

K := dRm / d \ddot{u}_{n+1} = (dRm / du_{n+ \alpha_f)) * (du_{n+\alpha_f} / d \ddot{u}_{n+1}

The last term of the above equation equals \alpha_f beta \Delta t^2. See equations 159 and 166 in Bazilevs et al. (2008). See Eq. 194-196 of Bazilves et al. for the final update step.

Regarding the time-step dependence, there could be two reasons I could think of:

1. Predictor step: As we use a predictor - multi-corrector scheme to advance in time, you are not using U_n as your starting value but a predicted U_{n+1} based on the choice of \beta. The predicted value also depends on the velocity but I believe if you are zeroing out both acceleration An and velocity Yn, then the predicted U_{n+1} will equal U_n. See Eq. 179 of Bazilevs et al. (2008).

2. Linear solver convergence: Because you are scaling the stiffness matrix with \alpha_f \beta \Delta t^2, it is possible that the linear solver is converging to a different solution compared to the one that unscaled case. Normally, we would expect these differences to vanish as we use relative tolerances and a rescaling later in the update step, but it may not be the case for iterative solvers. You may use a finer tolerance for both linear and non-linear solvers, and/or try to directly work on the displacement increments instead of the increments in acceleration.

Please try these and see if they make sense.

[mrp089]

Thanks, @vvedula22, for that comprehensive reply! I think that makes a lot of sense and I've figured out the static time integration for now. Once G&R #32 is working, we can figure out a way of how to best integrate "static time integration" into svFSI.