What is inside FEAP implementation of the B-Bar method?

[aortizb]

Dear FEAP community:

I would like to understand what is really going on inside the B-bar implementation of FEAP. According to the book “The Finite Element Method for Solids and Structural Mechanics”, 7th Edition, which I understand is mostly where FEAP implementation relies, it is stated that when the pressure and theta fields are locally approximated (hence, static condensation can be performed) the three-field u-p-theta mixed formulation can be treated as per the B-bar method. Hence, in Section 2.6.3 the B-bar method is used to construct the tangent stiffness matrix. It seems very reasonable. So far, so good. I gave a try to implement in MATLAB what is given there.

I mention that I am using the J2 model for plasticity with kinematic isotropic and kinematic hardening that is given in Section 4.6.2 of book stated above. I am working on a 1-element test. The element is a 9-node quadrilateral FE. My first try was to compare the standard displacement-based formulation. I looked at the stresses at the nine interior Gauss points. The results of my implementation in MATLAB perfectly match the results of FEAP. The same thing happens with nodal displacements. So far, so good.

However, when I do the same test using the mixed B-Bar formulation (9-node quadrilateral with 3 interior pressure/theta constraints, i.e., Q9/3 element as the name given in the book) the stresses and displacements do not match. I invested a lot of time in finding the source of this difference and debugging the code. I developed the equations in the book by myself to check whether there is a typo; but nothing; everything seems to be OK in the code and in the equations. Suddenly, I thought of doing the same exercise for linear elasticity. Then, I go to the book “The Finite Element Method Its Basis and Fundamentals”, 7th Edition. There, in Section 10.4, the u-p-theta approach is developed. In Section, 10.4.1 it is detailed the B-Bar method. Here, I found the source of difference between my code and what is implemented in FEAP… at least in linear elasticity: implementation of Equation (10.30):

\bar{A} = \int_Omega B’*D_d*B d\Omega + W’*H*W (IMPLEMENTATION 1)

gives different stresses at the nine interior Gauss points, and different nodal displacements, than the B-Bar implementation (Equation (10.31)):

\bar{A} = \int_Omega \bar{B}’*D*\bar{B} d\Omega,    (IMPLEMENTATION 2)

where,

\bar{B} = I_d*B + 1/3*m*N_v*W,

D = D_d + K*mm’.

I found that IMPLEMENTATION 1 perfectly matches the FEAP results, but not IMPLEMENTATION 2. Then, I plotted the entries of the stiffness matrices of IMPLEMENTATION 1 and 2 and compared them entry by entry. I found that there are differences in the entry values. The book, in the second paragraph below Equation (10.33), states that Equations (10.30) and (10.31) are equivalent. I find this to be, in theory ok, but numerically is apparently not as I get different entries in the stiffness matrices.

In conclusion, for the linear elastic case, what is implemented in FEAP matches with IMPLEMENTATION 1 and not with the B-Bar method (IMPLEMENTATION 2). So, could it be possible that something similar is going on with the elastoplastic case? That is, the B-bar method for the elastoplastic case as given in Section 2.6.3 of “The Finite Element Method for Solids and Structural Mechanics”, 7th Edition, is not exactly what is implemented in FEAP? I tried to figure out in files sldXdX.f, where X is a number, the implementation, but it seems to me that the implementation has many more things than the simple B-bar method described in Section 2.6.3. So, what is exactly implemented in FEAP when using the mixed 9-node element? I don’t find this to be critical as I tested the B-bar implementation of Section 2.6.3 on another problem (the classical elastoplastic Cook´s membrane) and found that as the mesh is refined, my B-Bar implementation based on Section 2.6.3 converges to the expected FEAP solution. The only thing is that my B-bar implementation needs about twice the number of 9-node elements than the elements in FEAP.

Any comments/corrections/additions would be appreciated.

Best regards,
Alejandro.

[Prof. S. Govindjee]

I have not looked at this in awhile.  But I do recall that there are some subtleties in the implementations in FEAP so that the constitutive models can be used without needing to have explicitly split pressure-volume forms.

Notwithstanding, have you had a look at the FEAP theory manual (chapter 7): http://feap.berkeley.edu/wiki/index.php?title=Manuals
to see how that compares to the Z&T volumes and what you see in your tests.

[aortizb]

Thank you so much Prof. Govindjee. I had not looked at the theory manual. But now I had it and found that the implementation for the linear elastic case is exactly the one I named IMPLEMENTATION 1 in my original post. So, there is no B-Bar implementation in FEAP... at least in the form that is described in the book. Regarding the nonlinear case, I found that the theory manual differs a lot from the book in what respects to mixed formulation and so I tend to think that FEAP's implementation is not the B-bar that is described in the book. I will try to implement what is in the theory manual and see if that resolves my doubts. Once again, thank you so much.

Alejandro.