Residual computation under finite displacement

[Enginit]

Dear all,

I have a question related to the computation of a residual under a finite displacement.
In FEAP, the residual terms (for body force and traction) are computed at the current or initial configuration?

Thank you very much in advance.

[Prof. R.L. Taylor]

You can compute at either configuration provided the quantities you use are defined properly in that configuration.  For example, commonly the mass matrix is computed in the reference configuration since one knows the density there.  On the other hand if you are computing a follower force it is often easier to compute in the current configuration.

Similarly, you can use the second PK stress to compute the stress divergence in the reference configuration or the Cauchy stress if you compute in the current configuration.  Some of this is described in the FEA books (e.g., the Zienkiewicz, Taylor, Fox volume   
               "The Finite Element Method for Solid and Structural Mechanics"

[Prof. S. Govindjee]

If you are asking about FEAP's built-in elements, then it depends on which one and with what time integrator if your problem is dynamic.  If it is a static problem most likely it will be current.

[Enginit]

Dear Prof. Taylor and Prof. Govindjee,

Thank you very much for response and suggestion. I would like to discuss further as below. All problems are in FEAP v8.4.

A. Problem 1: (brick 8-node, 8 gauss points, finite deformation, static)
I established the weak governing equation at the initial configuration (using 2 PK stress) and from that, I developed the residual and stiffness matrix. I employed the SHP3D to get the derivative N,x; N,y and N,z (in figure 1) at each integration point. For numerical integration, I employed INT3D to get sw(4,l) (weight) and SHP3D to get xjac at each integration point.
So my question is:
When the shp3d computes the values (N,x...xjac), which configuration the shp3d is using? Whether it uses the initial coordinates (initial configuration) or the updated coordinate (current configuration).

B. Problem 2: (brick 8-node, 8 gauss points, finite deformation, static)
FEAP brick element behaves very well under patch test and another example with forced displacement. Meanwhile, my element converges slower.
However, I encountered a strange behavior of this element in a simple cantilever beam with external force at the far right end. Please view figure 2 (FEAP element) and figure 3 (my element, it still works but converges slower, energy decreases half after one iteration). The input file is attached.
So my question are:
1. Could you check the input file? Is there any mistake there?
2. Is it possible for me to code the residual terms (for external loads) at the initial configuration in FEAP. If it is possible, how can I make it? (since FEAP computes the residual terms (for external loads automatically internally)).

[Prof. S. Govindjee]

A.  If you follow the source code you will see which configurations are being used.  I suggest also making a 1 element mesh and examining the contents of the arrays that are being used to get a concrete picture of how the code works.

B. If you want to apply your own loads, one option is to write an element that does the loading.  This in fact is the way we implement follower and other pressure loads; see the pressure elements for an example.

[Prof. S. Govindjee]

FYI.  B.  You are solving a non-linear problem; you can not put so much load on the system in one step.  Try putting the load on in 10 increments and you will see that FEAP's element works just fine.

Code: [Select]

       1    5   0  1.000E+00  1.00E+00  5.00E+03  3.78E-07  1.02E+04  3.04E-20      0.16 ER
       2    5   0  2.000E+00  1.00E+00  5.00E+03  2.69E-04  1.00E+04  4.99E-13      0.27 EN
       3    6   0  3.000E+00  1.00E+00  5.00E+03  9.25E-09  9.50E+03  7.10E-24      0.41 ER
       4    6   0  4.000E+00  1.00E+00  5.00E+03  1.15E-08  8.78E+03  1.13E-23      0.55 ER
       5    6   0  5.000E+00  1.00E+00  5.00E+03  1.25E-08  7.93E+03  1.40E-23      0.68 ER
       6    6   0  6.000E+00  1.00E+00  5.00E+03  1.61E-08  7.07E+03  2.15E-23      0.81 ER
       7    6   0  7.000E+00  1.00E+00  5.00E+03  1.69E-08  6.24E+03  2.53E-23      0.94 ER
       8    5   0  8.000E+00  1.00E+00  5.00E+03  8.64E-05  5.48E+03  1.53E-13      1.05 ER
       9    5   0  9.000E+00  1.00E+00  5.00E+03  3.17E-05  4.81E+03  1.92E-14      1.17 ER
      10    5   0  1.000E+01  1.00E+00  5.00E+03  1.07E-05  4.23E+03  2.01E-15      1.27 ER

Code: [Select]

   Residual norm =     5.0000000E+03    1.0000000E+00      t=     0.03     0.02
   Residual norm =     2.4663043E+05    4.9326086E+01      t=     0.05     0.02
   Residual norm =     1.7402965E+03    3.4805931E-01      t=     0.08     0.02
   Residual norm =     2.4654287E-01    4.9308574E-05      t=     0.10     0.02
   Residual norm =     3.7763763E-07    7.5527527E-11      t=     0.12     0.02
   Residual norm =     5.0000000E+03    1.0000000E+00      t=     0.14     0.02
   Residual norm =     2.3995298E+05    4.7990596E+01      t=     0.17     0.02
   Residual norm =     1.6507200E+03    3.3014400E-01      t=     0.19     0.02
   Residual norm =     2.2126800E+01    4.4253599E-03      t=     0.21     0.02
   Residual norm =     2.6870677E-04    5.3741354E-08      t=     0.23     0.02
   Residual norm =     5.0000000E+03    1.0000000E+00      t=     0.26     0.02
   Residual norm =     2.2159458E+05    4.4318917E+01      t=     0.28     0.02
   Residual norm =     1.4189887E+03    2.8379774E-01      t=     0.30     0.02
   Residual norm =     6.4145758E+01    1.2829152E-02      t=     0.32     0.02
   Residual norm =     7.2117766E-04    1.4423553E-07      t=     0.35     0.02
   Residual norm =     9.2459832E-09    1.8491966E-12      t=     0.37     0.02
   Residual norm =     5.0000000E+03    1.0000000E+00      t=     0.39     0.02
   Residual norm =     1.9579750E+05    3.9159500E+01      t=     0.41     0.02
   Residual norm =     1.1273820E+03    2.2547641E-01      t=     0.44     0.02
   Residual norm =     9.2890347E+01    1.8578069E-02      t=     0.46     0.02
   Residual norm =     8.9451634E-04    1.7890327E-07      t=     0.48     0.03
   Residual norm =     1.1464053E-08    2.2928105E-12      t=     0.51     0.03
   Residual norm =     5.0000000E+03    1.0000000E+00      t=     0.53     0.03
   Residual norm =     1.6727268E+05    3.3454535E+01      t=     0.55     0.03
   Residual norm =     8.5023769E+02    1.7004754E-01      t=     0.57     0.03
   Residual norm =     9.5154257E+01    1.9030851E-02      t=     0.59     0.03
   Residual norm =     7.2184678E-04    1.4436936E-07      t=     0.61     0.03
   Residual norm =     1.2455248E-08    2.4910495E-12      t=     0.64     0.03
   Residual norm =     5.0000000E+03    1.0000000E+00      t=     0.66     0.03
   Residual norm =     1.3963478E+05    2.7926956E+01      t=     0.68     0.03
   Residual norm =     6.2534398E+02    1.2506880E-01      t=     0.70     0.03
   Residual norm =     7.8272258E+01    1.5654452E-02      t=     0.72     0.03
   Residual norm =     4.3353287E-04    8.6706573E-08      t=     0.74     0.03
   Residual norm =     1.6144237E-08    3.2288475E-12      t=     0.77     0.03
   Residual norm =     5.0000000E+03    1.0000000E+00      t=     0.79     0.03
   Residual norm =     1.1493354E+05    2.2986708E+01      t=     0.81     0.03
   Residual norm =     4.5818930E+02    9.1637861E-02      t=     0.83     0.03
   Residual norm =     5.5338220E+01    1.1067644E-02      t=     0.85     0.03
   Residual norm =     2.0976604E-04    4.1953208E-08      t=     0.88     0.03
   Residual norm =     1.6895044E-08    3.3790087E-12      t=     0.90     0.03
   Residual norm =     5.0000000E+03    1.0000000E+00      t=     0.92     0.03
   Residual norm =     9.3935452E+04    1.8787090E+01      t=     0.94     0.03
   Residual norm =     3.3855300E+02    6.7710600E-02      t=     0.97     0.03
   Residual norm =     3.5137845E+01    7.0275691E-03      t=     0.99     0.03
   Residual norm =     8.6447440E-05    1.7289488E-08      t=     1.01     0.03
   Residual norm =     5.0000000E+03    1.0000000E+00      t=     1.03     0.03
   Residual norm =     7.6618800E+04    1.5323760E+01      t=     1.05     0.03
   Residual norm =     2.5325586E+02    5.0651171E-02      t=     1.07     0.03
   Residual norm =     2.0680132E+01    4.1360263E-03      t=     1.10     0.03
   Residual norm =     3.1656162E-05    6.3312323E-09      t=     1.12     0.03
   Residual norm =     5.0000000E+03    1.0000000E+00      t=     1.14     0.03
   Residual norm =     6.2583833E+04    1.2516767E+01      t=     1.16     0.03
   Residual norm =     1.9160580E+02    3.8321160E-02      t=     1.18     0.03
   Residual norm =     1.1553867E+01    2.3107733E-03      t=     1.21     0.03
   Residual norm =     1.0654846E-05    2.1309692E-09      t=     1.23     0.03


[Prof. R.L. Taylor]

If you are coding your own element you should follow the routines you are using to see how the computations are performed, for example when you call the shape function routine you pass coordinates.  If you follow how these are use you will learn how each routine uses them and understand what the result is.

These are things your faculty supervisor can help you.  Can you tell us who is the license holder for the version you are using that can help us guide you if you need help in other areas.

[Enginit]

Dear Prof. Govindjee and Prof. Taylor,

Thank you very much for your advice and responses. It helps me a lot.
I read the suggested book and theory manual of FEAP so I know that FEAP computes Cauchy stress (the governing equation is pushed forward the current|deformed configuration).
Regarding the license holder, please check your message box

Now, everything is clearer to me. I would like to discuss as follows.
1. "A.  If you follow the source code you will see which configurations are being used.  I suggest also making a 1 element mesh and examining the contents of the arrays that are being used to get a concrete picture of how the code works."
" If it is a static problem most likely it will be current."
Thank you for your suggestion, I checked the subroutines fld3d1, kine3d so I understand that FEAP is coded the same as the theory manual (at current|deformed configuration).

2. "B. If you want to apply your own loads, one option is to write an element that does the loading.  This in fact is the way we implement follower and other pressure loads; see the pressure elements for an example."
Could you tell me the name of pressure element subroutine? I don't know where it is in FEAP folder.

3. Is that true that user can not code the residual terms (related to external loads) so that it is computed in the reference|initial configuration in user element elmt**.f because it's only could handled by FEAP?
If the user can do that, could you give me instruction?
If the user can not do that, I will probably need to formulate everything at the current|deformed configuration (since FEAP already computes the residual terms (related to external loads) at the current|deformed configuration).

Best regards

[Prof. R.L. Taylor]

2.  Pressure element is in ./elements/frame/presld.f

3. What type of load do you want to apply?  That can help answer your question.  Also you should talk to your Professor Ibrahimbegovic as he can help you too -- or another member of the staff working on finite element methods.

[Enginit]

Dear Professor,
Thank you very much for your guidance and support. I will try to think about it.
Best regards,

[Enginit]

Dear Prof. Govindjee and Prof. Taylor,
I have just made my element working properly. I did a mistake on my geometry tangent matrix since I misunderstand that 2nd Piola Kirchhoff stress can be decomposed to S=Cd and then I put d outside to be the unknown, which is totally wrong, the displacement matrix d should be taken out from the variation strain instead.
Now the results from my element match with FEAP built-in element (large deformation by Saint Vernant Kirchhoff model).  In FEAP, the residual for finite deformation probably computed at the current configuration, as mentioned in "The Finite Element Method for Solid and Structural Mechanics" that it's a good way to set governing equations at the current configuration. Because this would keep B matrix the same as small deformation theory and also for other conveniences.
Once again, I would like to thank you for all of your comments and guidance. I would like to close this topic also.
Best regards,