surface load behaves different for user element and standard feap element

[crepes]

Hi,

some time ago I created a (modified-) neo-hookean model, mostly based on

Code: [Select]

$FEAPHOME8_5/ver85/elements/solid2d/*
$FEAPHOME8_5/ver85/elements/material/finite/*

The convergence is quadratic. For displacement boundary conditions, the user element calculates the same resulting displacements and stresses as the built-in standard feap element.

But for traction-based loads like in example12 from the inputs85.zip, e.g. with

Code: [Select]

csurf
  surface
    1 0 0 a -p
    2 a 0 a -p
    3 a a a  p
    4 0 a a  p

the calculated displacements are not the same in the user element compared with the standard element.

As far as I know, I do not have to implement something for traction boundary conditions in the user-element code, but I start getting sceptical.

Am I missing something?

Kind regards

[Prof. R.L. Taylor]

It is hard to know what the problem is based on your description.

1. What test did you perform with displacement b.c. -- be sure it was similar to the one you do with traction, not a homogenous tension or other uniform deformation.

2. It is odd you used the linear element to make the development.  Is your model based on Green strain?  If it is referential then maybe the traction loading will be different.  One thing to check is if the reactions from the displacement loading are the same for both elements tested.

3. It would help if you described how you did the development, which feap element you tested against. Example 12 seems to be 3-d in the test suite -- post the input file you used and the two output files.

[crepes]

Thank you very much for your reply.

My explanation was a little misleading. The model is indeed 3D, even though I mentioned 2D-routines. I also used 3D-routines. My first version of this model was in 2D plane strain and then later I changed it to 3D. And I think I only used finite deformation routines like fld2d1.f, fld3d1.f, modlfd.f, neoh3f.f and fengy3.f.

Your first suggestion helped to transform the problem into another one. I actually only had checked displacement b.c. by using a constant edge displacement. Now using a similar constant traction load reveals that it also gives the correct results for this kind of loads.

So the problem is a different one that I originally assumed.

The problem seems to only occur in the presence of non-homogeneous stress fields. I am testing against the MNEO element with volume function 2. I also use this volume function in my user element which is also modified neo-hookean law.

Since the convergence is still approximately quadratic, I might have messed up something that affects the residual and the tangent in the same way. I will try to track the problem based on these new informations.

I am not sure, if the input file and output files are still relevant. I attach them anyway. For this I used some kind of simple shear deformation. Here the results differ pretty much.

[Prof. R.L. Taylor]

I would debug with a much smaller mesh  2 x 2 x 2 or 1 x 1 x 1.

Probably in how you set up the shear behavior if uniform axial stress works.

[crepes]

Hello,

after spending several days on debugging, I finally found the problem. Under "isw.eq.1" I defined "d(182) = 2". According to pmanual, this is "Nodal quadrature parameters". I have thought that this represents the number of gauss points for each mesh dimension, but now I found out that it does not. According to inmate.f this is something else and it seems like it should be zero for gauss quadrature. Somehow I had set this to d(182)=2 several months ago and did not notice that it is a problem in some situations. While debugging, I found that shp3(...) is calculated wrong with d(182) = 2 in my case.

And it seems like quadr3d.f automatically sets 2x2x2 gauss points for brick elements, so I guess there is no need for defining the amount of gauss points in my element routine.

Without "d(182) = 2" it now works....

[Prof. S. Govindjee]

If you are writing your own elements with their own material models, then I highly recommend NOT using inmate.  There are simply too many things that can go wrong.

[crepes]

Thank you very much for your reply.

I, again, was too fast and careless/lazy here. I think I do not use inmate.f. Yesterday I thought that it was used automatically when creating a user element, but apparently, I was wrong. I should try to track the actually used subroutines and find out how the quadrature order is set in user elements.

[Prof. R.L. Taylor]

Feap does not set any quadrature or shape functions for user elements.  It is your option to set these as input data in the material set and choose which shape functions you want.

You can still use the quadrature routines that are in the solid2d, etc. directories but you need to set the d(*) parameters in your input data.  Download some of the user elements at the web site to see you they set the parameters and shape functions.  Most do not use inmate.f