contact and its assembly

[leifamer]

Hello dear FEAP users,

I am looking for some possibility to define surface to surface contact on 8 node brick element. It results a stiffness matrix to couple two surface, each with 4 nodes, so that the dimension of the stiffness matrix is 8*3x8*3=24x24.

Now there are some issues that I can’t understand, thank you first of all for any help.

1.) In the Feap Example Manual Chapter 8, it calculates the contact between two lines with non-matching mesh. Why is "TANG,,1" used, instead of "UTAN"?

2.) In file pmacr1.f, the subroutine contact is called. If I assembly the stiffness matrix of contact element explicitly "after" the utan/tang command. Is that possible? for example

UTAN,,-1
NEWContact  Stiffness assembly
FORM
NEWContact  Force assembly
SOLV

(Note: use unfactored stiffness matrix, and contact is assembled after utan, both without factorization, NEWContact means the new defined user macro)


Thank you very much.
Best regards

[FEAP_Admin]

program/formfe.f is the routine that calls contact( ) for computing tangents/right-hand sides (depending on the value of isw).
Look for

Code: [Select]

c     Contact contributions

      call contact (isw)


[leifamer]

Thanks a lot, now I find the right position to insert some code. 

One further question to assembly is, if the utan,,1 is called, there is some operation on the stiffness matrix, and factored with 0.8 (User Manual Page 180). Why is it necessary? How can I prevent it?

(In file pptang.f, utan,,-1 is used to plot the matrix for matlab. I am looking for some similar procedure)

[Prof. R.L. Taylor]

There is no factor on assembly.  There may be factors of 0.8 for the case where line search is used with a TANG command (line search is not a useful when UTAN is employed).  FEAP does all the assembly steps before it does any factoring of the coefficient matrix.  First elements and then contact.

You should be cautious about surface to surface contact.  Facets can over lap and then there will be more than 8 nodes affected.  The only surface to surface routines included in the contact module are for 2-d ties of parts -- that is the parts are joined at the start and never change.  Once they change the complexity is much increased and we have not yet attempted to implement any of the algorithms in the literature ourselves. 

[leifamer]

Thank your for reply, Prof. Taylor.

For overlap case, there are only two or more separate 8-nodes-contact parts, because the overlapped areas have only influence on the corresponding 8 nodes. It means that the shape of the stiffness due to contact does not change.
But the pairing is unknown and thus very difficult. Till now I just implemented something about contact through using a full matrix, which is not efficient. Now I found that the main problem is how change the pointer mr(np(21)) and the same variable in called subroutine jp( for  assembly the stiffness matrix in profile format.

Now I have to use FEAP, and find a solution for it. It is great thankful, if you could just give me some tips, and I can try it out and see if it works. Could I may ask how to adaptively change the profiled stiffness matrix, according to the current contact elements? (Assuming that the nodes of slave-master pairs are known.)

Thanks for your time.
Best regards

[Prof. R.L. Taylor]

If you write your routine to use the same routines to set the equation structure as in the existing contact modules the equation structure should get reset.  These are "modprof" type routines and the assembly routines "constass".

When feap solves a contact problem it calls contact(103) which activates calls to the contact modules for csw = 103 and csw = 403.  These options are used to set the equation structure before the solve is performed.  If you run a 2-d node to surface problems with "contact debug" in the contact setup you will get some files that show you the sequence of calls to contact subroutines.  Perhaps this will help you understand how the routines are working.

[leifamer]

Hello Prof. Taylor,

thank you again for your great advice. I have changed the integer pointer position mr(np(21)) according to the connection of two surfaces, at least it works now for one simple example with two blocks.

For feedback, I’d like explain a little what I have done. The current contact algorithms are implemented independent of the available subroutine().

1.)   For the change new profile due to contacted node/surface, a subroutine “myContactUpdate” is placed after contact(103). First of all a contact search of connected surface is carried out. After that the profile of the stiffness matrix is reset (call rstprf.f) and updated (call conprof.f), after that the new numbering for the equation is generated by calling nwprof().
(I know it is not so correct, and I should not use the available contact element)

2.)   In file formfe.f, I added a subroutine “myContactAlgorithm” analog to the last step after contact(isw). In the new subroutine there is only one assembly process. The stiffness matrix for the connection is generated and the constass() is called.

The graphic output is still a problem. The history variables for contact are also very difficult to understand. There should be a way to do that. Thank you for the help so far! In attachment is some plot of the example, still now I calculated it only with 8 dofs, and for stick case.