[midpoint scheme]

[Xuannamdo-danang]

Dear Professor Taylor,

Currently, I intended to use mid-point scheme for visco plasticity softening model and am reading your paper (please see attached file). However, there are some problems I can not understand and would like to ask you:

1. How dose FEAP update displacement, velocity, strain... at tn+1/2 for mid-point scheme in dynamic problems?

2. Obviously, in comparison with original formulation, stiffness matrix is changed. So, what/how is the change of mass matrix when using the mid-point scheme for computations?

3. In formulations, you only computed at tn+alpha (alpha=1/2 in case of mid-point scheme) for parameters involving viscoplastic strain rates while others are still computed at tn+1. Is that right? Why? Because in midpoint scheme for dynamic problems, we need to solve this equations Ku(tn+1/2) = fint(tn+1/2) - fext(tn+1/2) - Ma(tn+1/2).

Thank you so much and I am waiting for your reply

Yours sincerely,
DO XUAN NAM

[Prof. R.L. Taylor]

1. Generally, feap computes solutions at t_n+1 from values at t_n; for example the HHT algorithm.  This means  we need all terms from the increments of the n+1 values -- if you have constitution then one has two options: (1) interpolate the displacement to the n+alpha point or (2) interpolate the stress to the n+alpha point (from the n and n+1 values).

Currently, there may be some implementation for only (1) - but it is recommended to do (2) as it uses constitution always at n+1. It requires knowing the effect from the t_n solution too (e.g., store the element residual without inertia as history variable)

In either case stored solutions can be at t_+1

2.  The algorithm will tell how to do each of the terms: inertia, internal and external forces.

3. We always interpolated the alpha from the n and n+1 points as in the table in the paper.

[Xuannamdo-danang]

Dear Professor Taylor,

Thank you so much for your quick reply. However, it has not yet been obvious for me in several points as follows:

1. You mean that we should choose option 2 which uses constitution always at n+1 and then the stress is only one field which should be interpolated to the n+alpha point (from the n and n+1 values). However, from my point of view, in addition to the stress, the strain should also be interpolated, leading to the need of interpolation of displacement field at n+alpha point because computational procedure at element level needs to compute "viscoplastic multiplier-gamma" at n+alpha point to determine whether current step is elastic or viscoplastic (please see attached "beta" file). Is that right?

2. From the algorithm shown in Table 2, I see that delta alpha should be zero. Is that right? (please see attached "delta alpha" file)

Could you please explain to me

Thank you so much

Yours sincerely,
DO XUAN NAM

[Prof. R.L. Taylor]

The best way for you to learn algorithms is to implement a few and see how they work.  Once you have done the implementation you will understand all the steps needed to get the solution values.  Discuss each one with your supervisor, professor or fellow worker and you will gain a lot of knowledge.

[Prof. S. Govindjee]

Have a look at this message and the linked report:

http://feap.berkeley.edu/forum/index.php?topic=1121.msg3524;topicseen#msg3524

[Xuannamdo-danang]

Thank you so much for your replies, Professor Taylor and Professor Govindjee

[Xuannamdo-danang]

Dear Admins,

1. I read the document you sent me and already obtained a good result for Wilson-theta scheme. However, I do not exactly understand what FEAP has done. How can FEAP recognize ud(n,1) = displacement at t+theta*dt,  ud(n,2) = velocity at t+theta*dt, ud(n,3) = acceleration at t+theta*dt, ud(n,4) = velocity at t+dt, ud(n,5) = acceleration at t+dt? If I change the number "*" in ud(n,*) what will happen?

2. Can we control so that FEAP can take external force at t(n+1/2) or t(n+1) or generally at t(n+i) according to our demand for computing the following residual:

residual(t_n+1) = external force(t_n+1/2) - internal force(t_n+1) - inertia force(t_n+1)...?

Thank you so much

Yours sincerely,
DO XUAN NAM

[Prof. R.L. Taylor]

I am not sure this is possible, or a consistent algorithm?  The applied nodal forces are added to the residual in pload.f and interpolated using theta(3), a parameter form the time integrator.  But this is set to be 1.0 to compute the internal forces so there is no separate parameter to do what you want.  Furthermore, some of the external loading is computed by residuals (e.g., body forces, etc.) so they would be done at n+1

Not a good idea to modify the code for this in my opinion.

[Xuannamdo-danang]

Thank you so much for your reply, Professor Taylor