How to output contact force for those nodes in contact?

[Tao Jin]

Hi All,

I am doing some contact formulation and I want to output contact force on every mesh node which belongs to the slave surface and is in contact with the master surface.

I tried to use command 'reac' to output contact forces, but because of every node was in balance, the contact force can not be shown.

Does anyone know how to output contact force? By the way, I can get the node number of the node which is in contact.

Thank you

Tao

[Prof. R.L. Taylor]

it is currently not easy to just output the forces.  If you are in interactive mode you can see the contact state for every contact condition using SHOW CH -- in interactive mode -- then enter R and it will run through all the points.  You can scan the output file to see the values for all the contact variables -- one will be the normal force, others are gap, etc.

YOu can also try to get the reactions by first entering the solution command CONT OFF followed by reactions you want to see.  Be sure to activate contact again before continuing any solution (use CONT ON).

[Tao Jin]

it is currently not easy to just output the forces.  If you are in interactive mode you can see the contact state for every contact condition using SHOW CH -- in interactive mode -- then enter R and it will run through all the points.  You can scan the output file to see the values for all the contact variables -- one will be the normal force, others are gap, etc.

YOu can also try to get the reactions by first entering the solution command CONT OFF followed by reactions you want to see.  Be sure to activate contact again before continuing any solution (use CONT ON).

Thank you, Professor.

Here is the output of history variables for one contact node:

  Dump of set  #                1 Internal set #    60
          History set CH1         History set CH2  # in CH1   Absol #  Name
    0.000000000000000E+00   1.290000000000000E+02         1         1  masts   
    0.000000000000000E+00   1.472000000000000E+03         2        25  nachbar
    0.000000000000000E+00   1.000000000000000E+00         3         4  istgn   
    0.000000000000000E+00   4.906749999999960E-02          4        19  a1     
    0.000000000000000E+00   0.000000000000000E+00         5
    0.000000000000000E+00   1.204500000001385E-03          6
    0.000000000000000E+00  -9.790027019122823E-01          7        24  xi     
    0.000000000000000E+00   6.666567048859136E-01          8

                                  History set CH3  # in CH3   Absol #  Name
                            9.375184868097138E-02         1        16  flaeche
                            0.000000000000000E+00         2         2  knflg   
                           -4.610765386400466E-04         3         9  gn     
                           -2.454042392442261E-02         4        18  nvec   
                           -1.479851557514163E-15         5
                            9.996988384476645E-01         6
                            0.000000000000000E+00         7        20  a2     
                            3.750000000000000E-01         8
                            5.551115123125783E-16         9
                           -4.610765386400466E-06        10        14  fn     
                            0.000000000000000E+00        11        13  faug


Here is the explanation for the meaning of those symbols such as fn, nvec in the feap code:
c       Load dictionary of history variables

c       masts   - MASTter Segment
c       istgt   - Index of STatus for GT
c       istgn   - Index of STatus for GN
c       gn      - Normal projection Gap  (+ if open)
c       nvec    - normal vector to master segment
c       xi      - ...
c       a1      - vector tangential to surface
c       a2      - vector tangential to surface
c       a11,a12,a22 - metric tensor
c       diffxi1 - difference of xi(1) in this timestep to previous timestep
c       diffxi2 - difference of xi(2) ...

c       var. and lin:    slave
c                        1.master
c                        2.master
c                        3.master
c                        4.master

Although the dictionary of history variables does not mention, I guess that 'fn' is the normal force of contact node, right?


In the aspect of your second suggestion, my understand is that the solution batch should be like following:

BATCh   !Non-Augmented method
  LOOP,time,100
    TIME
      CONTact off
      REAC all
      CONTact on
      LOOP,newton,20
        CONTact,CHECk
        TANG,,1
      NEXT,newton
    REACtion all
    Show CH
  NEXT
END

Thank you for your time again.

Tao

[fheinrich]

Dear all,

although this topic is already more than 2 years old, it is still very relevant to me. I therefore hope
someone could share his/her thoughts.

The goal:
In fact, I want to achieve exactly the same thing Tao Jin asked in the first post:

... I want to output contact force on every mesh node which belongs to the slave surface and is in contact with the master surface. ...

The approach:
I followed Prof. Taylors hint and wrote the contact variables to my output-file, using "show ch", "R".
This yields a similar result to what Tao Jin posted before. I can see the normal forces of the slave nodes
in contact in the output-file:

Code: [Select]

                                  History set CH3  # in CH3   Absol #  Name
                            3.110620609763960E-05         1        16  flaeche
                            0.000000000000000E+00         2         2  knflg   
                           -1.671406229323534E-04         3         9  gn     
                            2.428244390954400E-01         4        18  nvec   
                            2.170097731499787E-02         5
                            9.698274894854026E-01         6
                           -3.390727629047547E-04         7        20  a2     
                            3.751198684065231E-03         8
                            9.594241532918913E-07         9
                           [b]-4.186473237250383E+06        10        14  fn[/b]
                            0.000000000000000E+00        11        15  ft     
                            0.000000000000000E+00        12         5  ft3 
The problem:
As I need to read (and maybe modify) the normal forces during run-time, I need to access
the forces in the routine "cnts3d". According to the contact manual, chapter 2.4 "History data
management and assignment" I should be able to read the normal forces with the command
"CH3(P3(10),NELM".
Unfortunately, this does not work for a couple of reasons:

• To me it is not clear, why the second parameter is "NELM" (the element number).
  I am looking for the force-value at a certain node(!), therefore I would have expected something like
  a global node-number as an input.
• To find the correct values, I added the following write-command to the loop of CSW=103 and CSW=403
  (i.e., the ones where it is looped over all slave nodes):

Code: [Select]

        loop_counter = 1
        do loop_counter = 1,160
          write(iow,*) 'loop_counter in csw403 is', loop_counter
          write(iow,*) 'ch3(p3(10,loop_counter) is ', ch3(p3(10),loop_counter)
        end do
Note: The loop_counter runs to 160 as I figured out that this is the maximum
possible loop-counter for my particular problem.

The output-file for this looks like:

Code: [Select]

loop_counter in csw403 is          88
 ch3(p3(10,loop_counter) is   0.000000000000000E+000
 loop_counter in csw403 is          89
 ch3(p3(10,loop_counter) is   4.456191161027072E-312
 loop_counter in csw403 is          90
 ch3(p3(10,loop_counter) is   0.000000000000000E+000
 loop_counter in csw403 is          91
 ch3(p3(10,loop_counter) is   0.000000000000000E+000
 loop_counter in csw403 is          92

So in short, CH3(P3(10),loop_counter) does not give the same results as from
"SHOW CH", "R"

As always, any help is highly appreciated.
Many thanks in advance!

Best regards

[fheinrich]

Additional note:

I found that, when using "SHOW CH", "R", the data is written to the output-file through the routine "dumpmem.f". In this routine I wanted to learn how the forces are being extracted from "CH3". Here, CH3 is defined as ch3(*), i.e. a vector. This is different to how it is defined in the contact manual: ch3(lh3,*), i.e. a matrix.
This fact, together with the above mentioned question why a history variable should be accessed in a way ch3(p3(...),NELM) instead of ch3(p3(...),NODE) is a bit confusing.

[fheinrich]

Additional note #2:

In this thread http://feap.berkeley.edu/forum/index.php?topic=1008.0, Prof. Taylor explains that the parameter for the normal force in ch3 is "p3(14)" - in contrast to what I wrote above ("p3(10)").

However, accessing the normal forces from the routine cnts3d still does not work for me.
I will post an update as soon as I have new insights.

[Prof. R.L. Taylor]

 The table is setting w3(14) to be the normal contact force  'fn', this shows up as the 10th item output for your problem -- but if you changed models it may move.  What do you see iwhen you look at ch3(p3(14), ...)??

The dumpmem routine is treating the array as single subscripted, however, it is using offsets to access the data it wants.  I would not suggest you try to do this, it involves extracting the data from multiple tables which you could also output.  It would seem you could do the same thing that the 2-d contact does and get output through the csw.eq.204 option

[fheinrich]

Dear Professor,

many thanks for your fast reply, it already helped a lot! Below, I provide you with some more background-information
and in addition I would like to ask another question.

The background:
My current simulation are two horizontal plates with clamped edges, a small gap between the plates and a vertical force pushing downwards in the middle of the upper plate. Please find the image "contact_problem_FEAP_15102015.JPG" attached.

I wrote the following input-file to achieve the above simulation:

Code: [Select]

FEAP * *
0 0 0 3 3 8

! Author:   fheinrich
! Date:     15th of October 2015
! Purpose:  Testing of modifications made in cnts3d

!---------------------------------------------------------------------------
! Input parameters
!---------------------------------------------------------------------------
PARAmeter
l = 0.08
h = 0.002
b = 0.015
g = 0.001
h1 = h + g
h2 = h1 + h
Q3 = 30000
Q2 = 100000
m = 20.
n = 2.
o = 2.


!---------------------------------------------------------------------------
! Mesh generation
!---------------------------------------------------------------------------
MATErial 1
SOLId
ELAStic ISOTropic 165e9

MATErial 2
SOLId
ELAStic ISOTropic 165e9

BLOCk 1
CARTesian,m,n,o,1,1,1
1 0. 0. 0.
2 l 0. 0.
3 l b 0.
4 0. b 0.
5 0. 0. h
6 l 0. h
7 l b h
8 0. b h

BLOCk 2
CARTesian,m,n,o,0,0,2
1 0. 0. h1
2 l 0. h1
3 l b h1
4 0. b h1
5 0. 0. h2
6 l 0. h2
7 l b h2
8 0. b h2


!---------------------------------------------------------------------------
! Boundary conditions
!---------------------------------------------------------------------------
EBOU
1 l 1 1 1
1 0 1 1 1


!---------------------------------------------------------------------------
! External loads
!---------------------------------------------------------------------------
CFORCE
NODE l/2 b/2 h2 0. 0. -Q3

END


!---------------------------------------------------------------------------
! Contact settings
!---------------------------------------------------------------------------
CONTact

SURFace 1
  QUADrilateral
  BLOCK SEGMent
  1 0. 0. h1
  2 0. b h1
  3 l b h1
  4 l 0. h1

SURFace 2
  QUADrilateral
  BLOCK SEGMent
  1 0. 0. h
  2 l 0. h
  3 l b h
  4 0. b h

MATE 1
  STAND
    FRIC COUL 1E15



PAIR 1
  CEL1 1 2 ! Author's note to the reader: Please be aware that here, CELMT01 is basically a mirrored version of the cnts3d-routine. This was done for the reason that the routine shall be modified in the future
 
  SOLM LAGM 2e9


END contact


!---------------------------------------------------------------------------
! Execution
!---------------------------------------------------------------------------
BATCh
    LOOP,augment,1
      LOOP,newton,30
        CONTACT CHECK
        TANG,,1
      NEXT
      AUGMent
    NEXT
DISPlacement COORdinate 2 0
END

INTEractive
STOP
After running the executable, I typed "SHOW CH", "R" to write history-data to the output-file.
In this particular example, I get an output-like this: (only the most relevant part is shown)

Code: [Select]

[...]

    0 - Not in Contact           =    60
    1 - Number to a Facet        =     3
    2 - Number to a Concave Edge =     0
    3 - Number to a Convex  Edge =     0
    4 - Number to a Point        =     0

[...]

Dump of set  #                1 Internal set #    19
          History set CH1         History set CH2  # in CH1   Absol #  Name
    0.000000000000000E+00   9.000000000000000E+00         1         1  masts   
    0.000000000000000E+00   1.570000000000000E+02         2        25  nachbar
    0.000000000000000E+00   1.000000000000000E+00         3         4  istgn   
    0.000000000000000E+00   1.961951509095453E-03         4        19  a1     
    0.000000000000000E+00   1.939704956858794E-07         5
    0.000000000000000E+00  -2.946218998362117E-04         6
    0.000000000000000E+00   8.798281349398949E-01         7        24  xi     
    0.000000000000000E+00   9.973037430893493E-01         8
    0.000000000000000E+00   0.000000000000000E+00         9         3  istgt   
    0.000000000000000E+00   0.000000000000000E+00        10        10  diffxi1
    0.000000000000000E+00   0.000000000000000E+00        11        11  diffxi2
    0.000000000000000E+00   3.587287106820496E-02        12        26  gt     
    0.000000000000000E+00   7.494598191551141E-03        13
    0.000000000000000E+00  -5.382467150958952E-03        14
    0.000000000000000E+00  -4.222053078158567E+07        15        31  lagmu   

                                  History set CH3  # in CH3   Absol #  Name
                            3.103208359189389E-05         1        16  flaeche
                            0.000000000000000E+00         2         2  knflg   
                           -8.076007169683797E-08         3         9  gn     
                            1.485013399755753E-01         4        18  nvec   
                            3.706140486058205E-03         5
                            9.889052616647138E-01         6
                            4.766808248855988E-06         7        20  a2     
                            3.746835089438059E-03         8
                           -1.475790977928862E-05         9
                           -4.282165938900373E+07        10        14  fn     
                            0.000000000000000E+00        11        15  ft     
                            0.000000000000000E+00        12         5  ft3     

  [...]

  Dump of set  #                1 Internal set #    21
          History set CH1         History set CH2  # in CH1   Absol #  Name
    0.000000000000000E+00   1.000000000000000E+01         1         1  masts   
    0.000000000000000E+00   1.580000000000000E+02         2        25  nachbar
    0.000000000000000E+00   1.000000000000000E+00         3         4  istgn   
    0.000000000000000E+00   1.947935301587423E-03         4        19  a1     
    0.000000000000000E+00   2.153316070647483E-07         5
    0.000000000000000E+00  -1.200356589293877E-04         6
    0.000000000000000E+00   9.952859509758349E-01         7        24  xi     
    0.000000000000000E+00   9.970681676903157E-01         8
    0.000000000000000E+00   0.000000000000000E+00         9         3  istgt   
    0.000000000000000E+00   0.000000000000000E+00        10        10  diffxi1
    0.000000000000000E+00   0.000000000000000E+00        11        11  diffxi2
    0.000000000000000E+00   3.999532920708605E-02        12        26  gt     
    0.000000000000000E+00   7.494168524660495E-03        13
    0.000000000000000E+00  -5.657374204101396E-03        14
    0.000000000000000E+00  -3.733081816449547E+08        15        31  lagmu   

                                  History set CH3  # in CH3   Absol #  Name
                            3.098482557413283E-05         1        16  flaeche
                            0.000000000000000E+00         2         2  knflg   
                           -1.392519785017587E-06         3         9  gn     
                            6.150027655609480E-02         4        18  nvec   
                            1.115175863662407E-02         5
                            9.980447656607565E-01         6
                            1.121379559845283E-06         7        20  a2     
                            3.742907783539232E-03         8
                           -4.189087583315588E-05         9
                           -3.733109666845284E+08        10        14  fn     
                            0.000000000000000E+00        11        15  ft     
                            0.000000000000000E+00        12         5  ft3     

[...]

  Dump of set  #                1 Internal set #    23
          History set CH1         History set CH2  # in CH1   Absol #  Name
    0.000000000000000E+00   1.200000000000000E+01         1         1  masts   
    0.000000000000000E+00   1.590000000000000E+02         2        25  nachbar
    0.000000000000000E+00   1.000000000000000E+00         3         4  istgn   
    0.000000000000000E+00   1.961416308644787E-03         4        19  a1     
    0.000000000000000E+00  -1.933369751104703E-07         5
    0.000000000000000E+00   2.939095858227143E-04         6
    0.000000000000000E+00  -8.859770673460258E-01         7        24  xi     
    0.000000000000000E+00   9.973032219071830E-01         8
    0.000000000000000E+00   0.000000000000000E+00         9         3  istgt   
    0.000000000000000E+00   0.000000000000000E+00        10        10  diffxi1
    0.000000000000000E+00   0.000000000000000E+00        11        11  diffxi2
    0.000000000000000E+00   4.412088522839708E-02        12        26  gt     
    0.000000000000000E+00   7.494597623479893E-03        13
    0.000000000000000E+00  -5.388225139907011E-03        14
    0.000000000000000E+00  -4.282149786886012E+07        15        31  lagmu   

                                  History set CH3  # in CH3   Absol #  Name
                            3.103706815062365E-05         1        16  flaeche
                            0.000000000000000E+00         2         2  knflg   
                           -7.783684611665069E-08         3         9  gn     
                           -1.481896386613505E-01         4        18  nvec   
                            3.848751581424443E-03         5
                            9.889514740899489E-01         6
                           -4.607721649441815E-06         7        20  a2     
                            3.746931225051313E-03         8
                           -1.527256339658783E-05         9
                           -4.222068645527987E+07        10        14  fn     
                            0.000000000000000E+00        11        15  ft     
                            0.000000000000000E+00        12         5  ft3
This means, after the final iteration, there are only 3 slave-nodes that are in contact with the master surface.
This seems to correct to me, as the entry "fn" of the "internal sets" 19,21 and 23 are non-zero. All other
internal sets have a zero "fn" value.

Using this output as a benchmark, I tried to extract the forces myself from my CELMT01-routine (i.e. the cnts3d-algorithm).
To do so, I added two write-commands to the section where csw.eq.403:

Code: [Select]

     elseif (csw.eq.403) then

c       Loop over all surface 1 elements = Loop over all slave nodes

        fp(1)  = np(191)+mr(np(192)+nsurf1-1)
        do ke = 1,mr(fp(1))
          kset = ke

          write(iow,*) 'kset in csw403 is', kset
          write(iow,*) 'ch3(p3(14),kset) is ', ch3(p3(14),kset)

[...]

The output in the output-file for the above write-commands looks like this:

[code] kset in csw403 is          19
 ch3(p3(14),kset) is   -42220686.4552739     
...
 kset in csw403 is          21
 ch3(p3(14),kset) is   -373310966.684531     
...
 kset in csw403 is          23
 ch3(p3(14),kset) is   -42821659.3890012 [/code]

This means, your hint works and I can extract the same forces by using ch3(p3(14),...). Very good!

The other question: Finding the global node numbers
The other question I have now is: How can I find the global node numbers of the (in this case)
3 nodes that are in contact with the master surface? Tao Jin mentioned in the first post that he
solved this, but did not provide the information how he did it.
I tried to find the nodes by using:

ix1(1,19)
ix1(1,21)
ix1(1,23)

, but this yields clearly wrong global node numbers.

Again, many thanks in advance!

Best regards

[Prof. R.L. Taylor]

I inserted the following code into cnst3d.f for the csw.eq.204 option.  Does this do what you want?  It does not use the element numbers but does loop over the nodes of the pair.

!       Get printout flag and range

        call setcprt (ifprt,fel,lel)

!       Print title
        if (ifprt) then
          write (iow,2001) npair
          fp(1)  = np(191)+mr(np(192)+nsurf1-1)
          do ke = 1,mr(fp(1))

            kset  = ke
            ns    = mr(fp(1)+kset)
            xs(1) = x(1,ns) + u(1,ns)
            xs(2) = x(2,ns) + u(2,ns)
            xs(3) = x(3,ns) + u(3,ns)

            istgn = nint(ch2(p1(4),kset))
            if(istgn.gt.0) then
              write(*,'(a,i5,1p,4e12.4)') 'NS=',ns,xs,ch3(p3(14),kset)
            endif

          end do ! ke
         

[fheinrich]

Dear Professor,

many thanks for the fast reply! Again, it helped a lot and the problem is almost solved, but:
It seems that the output of the nodes is working - however, in my particular problem it gives
wrong node numbers, because the contact algorithm itself describes a non-plausible deformation
of the structure.

I did the following:

1. I added your code-snippet of csw.eq.204 to the cnts3d routine at csw 403.
(In fact, I added it to my CELMT01 routine, which is a mirrored version of cnts3d for
modification purposes, as I stated earlier). The relevant part of the code then looks like
this:

Code: [Select]

c-----[--.----+----.----+----.-----------------------------------------]
cc--> csw 403 <--> csw <--> 403 <--> csw <--> 403 <--> csw <--> 403     /
c     Reset profile for contacts                                       V
c-----[--.----+----.----+----.-----------------------------------------]

      elseif (csw.eq.403) then

c       Loop over all surface 1 elements = Loop over all slave nodes

!       Get printout flag and range
        call setcprt (ifprt,fel,lel)

!       Print title
        if (ifprt) then
          write (iow,2001) npair
          fp(1)  = np(191)+mr(np(192)+nsurf1-1)
          do ke = 1,mr(fp(1))

            kset  = ke
            ns    = mr(fp(1)+kset)
            xs(1) = x(1,ns) + u(1,ns)
            xs(2) = x(2,ns) + u(2,ns)
            xs(3) = x(3,ns) + u(3,ns)

            istgn = nint(ch2(p1(4),kset))
            if(istgn.gt.0) then
              write(iow,'(a,i5,1p,4e12.4)') 'NS=',ns,xs,ch3(p3(14),kset)
            endif

          end do ! ke
        end if
2. I run this input-file:

Code: [Select]

FEAP * *
0 0 0 3 3 8

! Author:   fheinrich
! Date:     15th of October 2015
! Purpose:  Testing of modifications made in cnts3d

!---------------------------------------------------------------------------
! Input parameters
!---------------------------------------------------------------------------
PARAmeter
l = 0.08
h = 0.002
b = 0.015
g = 0.001
h1 = h + g
h2 = h1 + h
Q3 = 50000
m = 20.
n = 2.
o = 2.


!---------------------------------------------------------------------------
! Mesh generation
!---------------------------------------------------------------------------
MATErial 1
SOLId
ELAStic ISOTropic 165e9

MATErial 2
SOLId
ELAStic ISOTropic 165e9

BLOCk 1
CARTesian,m,n,o,1,1,1
1 0. 0. 0.
2 l 0. 0.
3 l b 0.
4 0. b 0.
5 0. 0. h
6 l 0. h
7 l b h
8 0. b h

BLOCk 2
CARTesian,m,n,o,0,0,2
1 0. 0. h1
2 l 0. h1
3 l b h1
4 0. b h1
5 0. 0. h2
6 l 0. h2
7 l b h2
8 0. b h2


!---------------------------------------------------------------------------
! Boundary conditions
!---------------------------------------------------------------------------
EBOU
1 l 1 1 1
1 0 1 1 1


!---------------------------------------------------------------------------
! External loads
!---------------------------------------------------------------------------
CFORCE
NODE l/2 b/2 h2 0. 0. -Q3

END


!---------------------------------------------------------------------------
! Contact settings
!---------------------------------------------------------------------------
CONTact

SURFace 1
  QUADrilateral
  BLOCK SEGMent
  1 0. 0. h1
  2 0. b h1
  3 l b h1
  4 l 0. h1

SURFace 2
  QUADrilateral
  BLOCK SEGMent
  1 0. 0. h
  2 l 0. h
  3 l b h
  4 0. b h

PAIR 1
  CEL1 1 2 ! Author's note to the reader: Please be aware that here, CELMT01 is basically a mirrored version of the cnts3d-routine. This was done for the reason that the routine shall be modified in the future
  SOLM LAGM 2e9

END contact


!---------------------------------------------------------------------------
! Execution
!---------------------------------------------------------------------------
BATCh
    LOOP,augment,1
      LOOP,newton,30
        CONTACT CHECK
        TANG,,1
      NEXT
      AUGMent
    NEXT
DISPlacement COORdinate 2 0
END

INTEractive
STOP
After this, I checked the "nodes in contact" from the output-file:

Code: [Select]

     C o n t a c t   O u t p u t   f o r   P a i r     1
NS=  197  2.7557E-02  1.1590E-05 -7.8833E-03  0.0000E+00
NS=  198  3.1660E-02  1.3599E-05 -9.4701E-03  0.0000E+00
NS=  220  3.5788E-02  7.5137E-03 -1.0578E-02 -1.5843E+08
NS=  199  3.5815E-02  2.5157E-06 -1.0533E-02  0.0000E+00
NS=  221  3.9990E-02  7.5146E-03 -1.1056E-02 -4.6052E+08
NS=  222  4.4200E-02  7.5137E-03 -1.0590E-02 -1.6103E+08
NS=  201  4.4173E-02  2.6437E-06 -1.0542E-02  0.0000E+00
NS=  202  4.8330E-02  1.3810E-05 -9.4857E-03  0.0000E+00
NS=  203  5.2435E-02  1.1751E-05 -7.9013E-03  0.0000E+00
, and plotted them:

see attachment "contact_problem_FEAP_20102015_wrong_nodes.jpg".

This shows nodes in contact with a presumably wrong antisymmetry.

To validate this, I run the above input-file again with a higher force
(Q3 = 200000). The result-plot (see attachment "contact_problem_FEAP_20102015_higher_force.jpg")
then shows indeed that the beams are twisting in a way they are not expected to.

The question:
The major question is now, where does this twisting of the plates come from?
If the twisting was not present, I assume that the above code would give me the
correct "nodes in contact".

Some notes:
- I figured out that it has to do with the contact. Leaving the contact option out, does
not twist the upper plate
- It seems like a numerical effect: When I applied a force of Q3 = 50000 N, the antisymmetry
was locating the nodes mainly at the y-coordinate (2-coordinate) = 0. When I use the higher
force, the plates twist in the other direction.
- If you want to use my above input-file, you have to replace "CEL1 1 2" by "NTOS 1 2"

Many thanks for any input! Solving this last question will surely solve the whole question!

Best regards

[Prof. R.L. Taylor]

I have run your input file using the NTOS element.  To get a solution, I had to reduce the size of the load and add an incremental solution  so that there was not too much in contact and not too much penetration.  This allowed for convergence.  Basically,
I used the set of commands

prop
dt,,1
loop,,3
  time
  loop,check,3
     contact  check
     loop,,10
        tang,,1
     next
  next
   .... outputs and plots
next

I reduced Q3 to 2000 for this analysis and your other data.

There are still issues that point to an error in the implementation of the contact checking.  I do get the right list of nodes in contact using the STRE CONT 1 command, but the edge forces are positive which is not correct they should release.  I will have to do further checking to see why this is not happening.  But try the above to see if you are getting the same result first.

At present you are doing contact, but the solid elements are set to be linear.  I would suggest you make them non-linear elastic (use NEOHook instead of ISOT as a simple way to do this).

Let me know your results.

[fheinrich]

Dear Professor,

many thanks for the helpful reply!

I can already tell that, with your proposed modifications, the simulation converges and that the "nodes in contact" now seem reasonable. I am going to do more simulations in the next days and will post more detailed results in the course of next week.

Best regards and have a nice weekend!