Contact between matrix and inclusions (Should I specify contact pairs here?)

[MBrown]

Hello all,

I've created a mesh representing a circular cross section of one material inset with inclusions of another material (first attachment). Upon running, the mesh blows up, and doesn't make it past the first solution. I've also created a similar mesh, but where the inclusions are left as voids (second attachment). This model, although slow, seems to solve just fine. Nothing changes in the solution between the two models.

What could be the root of my problem? My working theory is that I might need to define contact between the two materials at the inclusion interface, but this hasn't been a problem in the past when modeling a quarter-symmetric inclusion (third attachment).

Thanks,
Michael

[Prof. R.L. Taylor]

1. Feap does not like meshes with missing elements (it is something we will try to fix).  If you insert a *ELEM -1564 before the
elements it element description, it ran for me.

2. The slow nature is basically because the profile is bad due to the wan the nodes have been numbered by the mesh generator.
If you add an OPTI in the batch solution it both of the in.S* meshes run much faster.

[MBrown]

Thank you! Both problems have been giving me trouble. What is causing the missing elements? Is it a problem with how I'm entering them in the input file?

[MBrown]

Actually, I'm sorry, but I've realized that I generated the StrandBuildMixedMatrix file from the other, and I accidentally saved the first. This is the real in.StrandBuildMixed file that keeps blowing up... I've implemented the changes that you suggested though, and I have to say it has definitely sped up. Anyway, do you know why it might be acting this way?

Thanks again,
Michael

[Prof. R.L. Taylor]

I do not know how you generate your meshes, but the result created errors in how the elements are numbered.

If you turn off contact and just put a point load at the top it does not work.

If you issue the command CHECk the programs imples some elements are incorrectly numbered. I do not know if all the diagnostics are meaningful since you have some quads and some triangles, however, after the check (feap tries to correct the incorrectly numbered elements). the  problem runs -- so does the contact.

However, you will not be able to solve the problem with plasticity when you use 3-node triangles.  These elements cannot deform consistently with a near volumetric response (this is a well known locking phenomena).  Even 4-node quads using a MIXEd form may not work well in highly confined problems.  In general you need to use higher order elements, however, this also has difficulties as the current contact in feap is not intended for higher order elements.  You may be able to still use the node to surface form if the SLAVE body uses 4-node elements and the disk higher order elements -- but test some before doing too many analyses.

[MBrown]

Thank you, this has been incredibly helpful. Especially using OPTImize and CHECk.

I'm using gmsh to generate meshes and I've switched to using 4 node elements now to stop the locking phenomena. If anyone is interested, I've attached another, more complicated model here, meshed in quad 4 node elements.

Interestingly (and annoyingly) I've found that it crashes fairly early when the FINIte MIXEd option is in effect, but when that option is turned off, it seems to run well. Could this be an indication of what you were talking about, that my elements are too low order? The error I get is *WARNING* no convergence in PLASFD. I assume this is referring to the routine in elements/material/finite.

[Prof. R.L. Taylor]

There are a couple of things I observe.

1. you seem to be using small deformation solid elements.

2. you have a material with nu = 0.5 -- FEAP is not coded to allow exactly incompressible.

3. the mixed elements are intended for cases which are nearly incompressible.

4. 4-node quads are not fully BB-consistent (for small strain); 9-node are for small strain -- what works in finite strain depends on being able to converge the solution.  Higher order elements can go unstable and result in a "tangled" mesh.  Nevertheless, for nearly incompressible behavior the 4-node elements may work -- but not for nu=0.5 exactly.

5. There is a bug in the routine 'ckiosp.f' in the /program directory.  The value of xsj returned from the call to shp2 should be assigned to jac(ineg) in the if check below the call to read
       if(xsj.le.0.0d0) then
         jac(ineg) = xsj
          .....

This will fix the fact that all jacobians are output as zeros.

This is a hard problem!  gmsh does a good job in meshing though.

[Prof. R.L. Taylor]

I have some success in solving with the finite mixed element if I use TANG LINE 1 for the solution instead of TANG,,1

You will find the deformed solution is quite different showing a lot of distortion

[MBrown]

Prof. Taylor,

I believe these pointers have, for the most part, done the trick. The last thing to add is plastic hardening, but it may not be necessary as the new models I am getting after correcting the error in the program file and using the tang line 1 command are approaching real word data. I will try to do a bit of parametric work and hone in on a more accurate model.

I will play around more with the poisson ratio as well. The 0.5 value was really a filler value. I'm trying to figure out how to approximate the behavior of a powder for that material in particular.

GMSH is also free! and from what I understand, it particularly shines for 3D meshes, and when using its post processing functions. But I've not had to do that yet.

Again, thanks,
Michael

[Prof. R.L. Taylor]

Michael,

Regarding your comments.

You could make the piece rigid by fixing all the nodes and not just the top row.  Then define displacements for all the nodes too.  In this way the material properties are really never used and you can put anything there.  It may help the conditioning when large deformations exist in the circular disk parts.

There is another plasticity model that exists in feap (It is the HILL model but may be made isotropic).  One should note that the kinematic hardening coded in the isotropic model is WRONG and we have removed it from future releases.  The kinematic hardening in the anisotropic model is good and may be used.  Isotropic hardening is acceptable in both models.

Your problem is one of the most challenging that I have observed!