finish
/clear,start

/sys,mkdir C:\ANSYS\double_edge_crack\3D\elastic
/cwd,C:\ANSYS\double_edge_crack\3D\elastic
*set,JName,'3D_elastic'
/filname,JName,0

/title, Außenrisse rechts und links, elastisch, lastgesteuert, K_J

!/unit,si

! Einheiten
! Masse        [kg]
! Längen       [mm]
! E-Modul      [N/mm²]
! Flächen      [mm²]
! Spannungen   [N/mm²]
! Querdehnzahl [-]

/VIEW,  1, -0.217341817295    , -0.872958677897    ,  0.436698615912
/ANG,   1,   35.1695878049
/REPLO

*ASK,t,Blechdicke (Default = 20 mm),20
*ASK,w,Blechbreite (Default = 260 mm),260
*ASK,l,Blechlänge (Default = 600 mm),600
*ASK,a,Anfangsrisslänge (Default = 10 mm),10
*ASK,F,Aufbebrachte Spannung in N/mm² (Default = 150),150

/Prep7

SMRTSIZE,off

!Elementauswahl
ET,1,solid45
ET,2,solid45
ET,3,solid95

!Material Eigenschaften
MP,EX,1,2.1e5
MP,ALPX,1,1.2E-5
MP,KXX,1,25.0
MP,DENS,1,7850.0
MP,PRXY,1,0.3
MP,C,1,600
MP,REFT,1,50
MP,HF,1,10

!Geometrie der Platte

!Generierung eines feinen Halbkreisförmigen Rissmeshes
*AFUN,deg
CSYS,0
N,1
CSYS,1
NGEN,9,20,1
N,11,0.8*a,0,0
N,171,0.8*a,180,0
FILL,11,171,7,31,20
CSYS,0
FILL,1,11,9,2,1,9,20,3

!Generiert eine bestimmte Anzahl an Knotenebenen, um die Elemente nicht
!zu sehr zu verzerren
NGEN,9,1000,1,171,1,,,(t/2)/8
E,2,22,1,1,1002,1022,1001,1001
EGEN,8,20,-1
EGEN,8,1000,1,8,1
E,2,3,23,22,1002,1003,1023,1022
EGEN,8,20,-1
EGEN,8,1000,65,72,1
E,3,4,24,23,1003,1004,1024,1023
EGEN,8,20,-1
EGEN,8,1000,129,136,1
E,4,5,25,24,1004,1005,1025,1024
EGEN,8,20,-1
EGEN,8,1000,193,200,1
E,5,6,26,25,1005,1006,1026,1025
EGEN,8,20,-1
EGEN,8,1000,257,264,1
E,6,7,27,26,1006,1007,1027,1026
EGEN,8,20,-1
EGEN,8,1000,321,328,1
E,7,8,28,27,1007,1008,1028,1027
EGEN,8,20,-1
EGEN,8,1000,385,392,1
E,8,9,29,28,1008,1009,1029,1028
EGEN,8,20,-1
EGEN,8,1000,449,456,1
E,9,10,30,29,1009,1010,1030,1029
EGEN,8,20,-1
EGEN,8,1000,513,520,1
E,10,11,31,30,1010,1011,1031,1030
EGEN,8,20,-1
EGEN,8,1000,577,584,1

!Speichert die Knoten in der Rissspitze
NSEL,S,LOC,X,0
NSEL,R,LOC,Y,0
CM,CRACKTIP,NODE
CMSEL,S,CRACKTIP
ESLN
TYPE,2
EMODIF,ALL
NUMMRG,NODE

!Solid95-Elemente in die Rissspitze legen
!Verschieben der midside-nodes in die Viertelpunkte
/NOPR
NSEL,ALL
*GET,N,NODE,,NUM,MAX
CMSEL,S,CRACKTIP                  !ruft die Knoten in der Rissspitze auf
ESLN                                   !wählt die dazugehörigen Elemente
*GET,ELMAX,ELEM,,NUM,MAX
*DO,IEL,1,ELMAX,1
     ELMI=IEL
     *IF,ELMI,LE,0,EXIT
     *GET,ELTYPE,ELEM,ELMI,ATTR,TYPE
     *IF,ELTYPE,NE,2,CYCLE
     N3 = NELEM(ELMI,3)
     *IF,NSEL(N3),LE,0,CYCLE
     N7 = NELEM(ELMI,7)
     *IF,NSEL(N7),LE,0,CYCLE
     N1 = NELEM(ELMI,1)
     N2 = NELEM(ELMI,2)
     N5 = NELEM(ELMI,5)
     N6 = NELEM(ELMI,6)

     X3 = 0.75*NX(N3)
     Y3 = 0.75*NY(N3)
     Z3 = 0.75*NZ(N3)
     X = 0.25*NX(N2) + X3
     Y = 0.25*NY(N2) + Y3
     Z = 0.25*NZ(N2) + Z3
     N = N + 1
     N10 = N
     N,N10,X,Y,Z
     X = 0.25*NX(N1) + X3
     Y = 0.25*NY(N1) + Y3
     Z = 0.25*NZ(N1) + Z3
     N = N + 1
     N12 = N
     N,N12,X,Y,Z
     X7 = 0.75*NX(N7)
     Y7 = 0.75*NY(N7)
     Z7 = 0.75*NZ(N7)
     X = 0.25*NX(N6) + X7
     Y = 0.25*NY(N6) + Y7
     Z = 0.25*NZ(N6) + Z7
     N = N + 1
     N14 = N
     N,N14,X,Y,Z
     X = 0.25*NX(N5) + X7
     Y = 0.25*NY(N5) + Y7
     Z = 0.25*NZ(N5) + Z7
     N = N + 1
     N16 = N
     N,N16,X,Y,Z
     N4 = N3
     N8 = N7
     NSEL,ALL
     TYPE,3
     EN,ELMI,N1,N2,N3,N4,N5,N6,N7,N8         !Redefinition des Elementes
     EMORE,0,N10,0,N12,0,N14,0,N16
     EMORE,
*ENDDO

CMSEL,U,CRACKTIP
NUMMRG,NODE
NSEL,ALL
ESEL,ALL
/GOPR

!Generierung der restlichen Platte
!Knotenkoordinaten
P1X=0
P2X=0.8*a
P3X=w/2-a
P4X=P3X
P5X=-a
P6X=P5X
P7X=-0.8*a
P8X=-cos(22.5)*0.8*a
P9X=-cos(22.5*2)*0.8*a
P10X=-cos(22.5*3)*0.8*a
P11X=0
P12X=cos(22.5*3)*0.8*a
P13X=cos(22.5*2)*0.8*a
P14X=cos(22.5)*0.8*a
P15X=P4X
P16X=P5X

P1Y=0
P2Y=0
P3Y=0
P4Y=l/4
P5Y=l/4
P6Y=0
P7Y=0
P8Y=sin(22.5)*0.8*a
P9Y=sin(22.5*2)*0.8*a
P10Y=sin(22.5*3)*0.8*a
P11Y=0.8*a
P12Y=sin(22.5*3)*0.8*a
P13Y=sin(22.5*2)*0.8*a
P14Y=sin(22.5)*0.8*a
P15Y=l/2
P16Y=l/2

K,1,P1X,P1Y,0
K,2,P2X,P2Y,0
K,3,P3X,P3Y,0
K,4,P4X,P4Y,0
K,5,P5X,P5Y,0
K,6,P6X,P6Y,0
K,7,P7X,P7Y,0
K,8,P8X,P8Y,0
K,9,P9X,P9Y,0
K,10,P10X,P10Y,0
K,11,P11X,P11Y,0
K,12,P12X,P12Y,0
K,13,P13X,P13Y,0
K,14,P14X,P14Y,0
K,15,P15X,P15Y,0
K,16,P16X,P16Y,0

K,8001,P1X,P1Y,t/2
K,8002,P2X,P2Y,t/2
K,8003,P3X,P3Y,t/2
K,8004,P4X,P4Y,t/2
K,8005,P5X,P5Y,t/2
K,8006,P6X,P6Y,t/2
K,8007,P7X,P7Y,t/2
K,8008,P8X,P8Y,t/2
K,8009,P9X,P9Y,t/2
K,8010,P10X,P10Y,t/2
K,8011,P11X,P11Y,t/2
K,8012,P12X,P12Y,t/2
K,8013,P13X,P13Y,t/2
K,8014,P14X,P14Y,t/2
K,8015,P15X,P15Y,t/2
K,8016,P16X,P16Y,t/2

!Linien
L,1,2
L,2,3
L,3,4
L,4,5
L,5,6
L,6,7
L,7,8
L,8,9
L,9,10
L,10,11
L,11,12
L,12,13
L,13,14
L,14,2
L,4,15
L,15,16
L,16,5
L,7,1
L,1,11

L,8001,8002
L,8002,8003
L,8003,8004
L,8004,8005
L,8005,8006
L,8006,8007
L,8007,8008
L,8008,8009
L,8009,8010
L,8010,8011
L,8011,8012
L,8012,8013
L,8013,8014
L,8014,8002
L,8004,8015
L,8015,8016
L,8016,8005
L,8007,8001
L,8001,8011

l,1,8001
l,2,8002
l,3,8003
l,4,8004
l,5,8005
l,6,8006
l,7,8007
l,8,8008
l,9,8009
l,10,8010
l,11,8011
l,12,8012
l,13,8013
l,14,8014
l,15,8015
l,16,8016

!Flächen
AL,7,8,9,10,18,19
AL,7,26,45,46
AL,8,27,46,47
AL,9,28,47,48
AL,10,29,48,49
AL,19,38,39,49
AL,18,37,39,45
AL,26,27,28,29,37,38

AL,1,11,12,13,14,19
AL,11,30,49,50
AL,12,31,50,51
AL,13,32,51,52
AL,14,33,52,40
AL,1,20,39,40
AL,20,30,31,32,33,38

A,6,3,15,16
A,6,3,8003,8006
A,3,15,8015,8003
A,15,16,8016,8015
A,16,6,8006,8016
A,8006,8003,8015,8016

!Volumen
VA,1,2,3,4,5,6,7,8
VA,9,6,11,12,13,14,10,15
VA,16,17,18,19,20,21

VSBV,3,1
VSBV,4,2

!Zuweisung der Eigenschaften zu den Volumen + Meshen
VSEL,ALL
VATT,1,1,1

LDELE,3
LDELE,4
LDELE,5
LDELE,15
LDELE,17
LDELE,22
LDELE,23
LDELE,24
LDELE,34
LDELE,36

LESIZE,ALL,,,1
LESIZE,2,,,w/10,8
LESIZE,21,,,w/10,8
LESIZE,16,,,w/10
LESIZE,35,,,w/10
LESIZE,56,10,,
LESIZE,57,10,,,1/3
LESIZE,59,10,,
LESIZE,60,10,,,3
LESIZE,40,,,-0.75*a+9.5
LESIZE,41,,,-0.75*a+9.5
LESIZE,42,,,-0.75*a+9.5
LESIZE,43,,,-0.75*a+9.5
LESIZE,44,,,-0.75*a+9.5
LESIZE,45,,,-0.75*a+9.5
LESIZE,46,,,-0.75*a+9.5
LESIZE,47,,,-0.75*a+9.5
LESIZE,48,,,-0.75*a+9.5
LESIZE,49,,,-0.75*a+9.5
LESIZE,50,,,-0.75*a+9.5
LESIZE,51,,,-0.75*a+9.5
LESIZE,52,,,-0.75*a+9.5
LESIZE,53,,,-0.75*a+9.5
LESIZE,54,,,-0.75*a+9.5

VSWEEP,ALL

!Verbinden der Knoten, die auf selber Position liegen
CPINTF,UX
CPINTF,UY
CPINTF,UZ
CEINTF,,UX
CEINTF,,UY
CEINTF,,UZ

!Lagerung
NSEL,S,LOC,X,w/2-a
DSYM,SYMM,X
NSEL,ALL
NSEL,S,LOC,X,0,w/2-a
NSEL,R,LOC,Y,0
DSYM,SYMM,Y
NSEL,ALL
NSEL,S,LOC,Z,0,0
DSYM,SYMM,Z
NSEL,ALL
NSEL,S,LOC,Z,t/2,t/2
D,ALL,UZ,0
NSEL,ALL

!Aufbringen der Spannung
NSEL,S,LOC,Y,l/2,l/2
SF,ALL,PRES,-F
NSEL,ALL
ALLS
FINISH

!Starten der Berechnung

/COM                  !Places a comment in the output
/OUTPUT,SCRATCH       !Redirects text output to a file or to the screen.


!Berechnung
/Solu
/nerr,5,99999999        !Limitiert die Anzahl von Warnungen pro Kommando
antype,static
NEQIT,50            !Specifies the max. number of equilibrium iterations
                    !for nonlinear analyses.

*ASK,N,Anzahl der Substeps (Voreingestellt sind 50 substeps),50
NSUBST,N                                        !Anzahl der Lastschritte
OUTPR,ALL
OUTRES,ALL,ALL
ALLSEL
SOLVE
SAVE
FINISH


!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!------------- Ausgabe der Dehnungen über die Breite -------------------
!------------- Lage 1: y = L/2 -----------------------------------------
!------------- Lage 2: y = L/4 -----------------------------------------
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------

/POST1
*DIM,strain,array,50,4,N+1,,,,
*DIM,f_res,array,N,1,,,,
*DIM,stress_strain,array,N,8,9,,,,

SET,First

CSYS,0

*DO,q1,1,N,1

       *DO,q2,1,21,1
       !Pfad 1(y=0 und z=0)                             TOTALE DEHNUNGEN
       PATH,Pfad1,2,30,20
       PPATH,1,,0,0,0
       PPATH,2,,w/2-a,0,0
       PDEF,Pfad1,epto,1
       *GET,Eps1,PATH,0,ITEM,Pfad1,PATHPT,q2
       *VFILL,strain(q2,1,q1),data,Eps1
       PADELE,Pfad1
       !Pfad 2 (y=L/2 und z=0)
       PATH,Pfad2,2,30,20
       PPATH,1,,-a,l/4,0
       PPATH,2,,w/2-a,l/4,0
       PDEF,Pfad2,epto,1
       *GET,Eps2,PATH,0,ITEM,Pfad2,PATHPT,q2
       *VFILL,strain(q2,2,q1),data,Eps2
       PADELE,Pfad2
       !Pfad 3(y=0 und z=0)                         PLASTISCHE DEHNUNGEN
       PATH,Pfad3,2,30,20
       PPATH,1,,0,0,0
       PPATH,2,,w/2-a,0,0
       PDEF,Pfad3,epel,1
       *GET,Eps3,PATH,0,ITEM,Pfad3,PATHPT,q2
       *VFILL,strain(q2,3,q1),data,Eps3
       PADELE,Pfad3
       !Pfad 4 (y=L/2 und z=0)
       PATH,Pfad4,2,30,20
       PPATH,1,,-a,l/4,0
       PPATH,2,,w/2-a,l/4,0
       PDEF,Pfad4,epel,1
       *GET,Eps4,PATH,0,ITEM,Pfad4,PATHPT,q2
       *VFILL,strain(q2,4,q1),data,Eps4
       PADELE,Pfad4
       *ENDDO

       *DO,q3,1,9,1
       *GET,Stress_y,NODE,(q3*1000)-999,S,Y
       *GET,Strain_el_y,NODE,1,EPEL,Y
       *GET,Strain_pl_y,NODE,(q3*1000)-999,EPPL,Y
       *GET,Strain_tot_y,NODE,(q3*1000)-999,EPTO,Y
       *GET,Stress_1,NODE,(q3*1000)-999,S,1
       *GET,Strain_el_1,NODE,(q3*1000)-999,EPEL,1
       *GET,Strain_pl_1,NODE,(q3*1000)-999,EPPL,1
       *GET,Strain_tot_1,NODE,(q3*1000)-999,EPTO,1
       *VFILL,stress_strain(q1,1,q3),data,Stress_y
       *VFILL,stress_strain(q1,2,q3),data,Strain_el_y
       *VFILL,stress_strain(q1,3,q3),data,Strain_pl_y
       *VFILL,stress_strain(q1,4,q3),data,Strain_tot_y
       *VFILL,stress_strain(q1,5,q3),data,Stress_1
       *VFILL,stress_strain(q1,6,q3),data,Strain_el_1
       *VFILL,stress_strain(q1,7,q3),data,Strain_pl_1
       *VFILL,stress_strain(q1,8,q3),data,Strain_tot_1
       *ENDDO

       !Querschnitt in die äußere Knotenebene legen und
       !Spannungsausgabe auf diesem Querschnitt
       KWPLAN,1,16,15,8016
       SUCR,SURFACE,CPLANE
       SUSEL,S,SURFACE
       SUMAP,F_Y,S,Y
       SUGET,SURFACE,F_Y,FORCE,OFF
       SUEVAL,F_R,F_Y,AVG
       f_res(q1,1)=F_R
       WPCSYS,1,0
       SUSEL,NONE

       SET,NEXT

*ENDDO

PARSAV,ALL,C:\ANSYS\double_edge_crack\3D\elastic\eps_Standard_plane_strain,txt,

!------------------------------------K1-WERTE---------------------------
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!----Bestimmung von K1-Werten mit unterschiedlichen Pfadpunkten---------
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------

/POST1
*DIM,KTable,ARRAY,N,9,,,,,
CSYS,0
SET,First

*DO,q4,1,N,1

         !Definition Pfad
         PATH,K_1,3,,48
         PPATH,1,1
         PPATH,2,8200
         PPATH,3,162
         KCALC,0,1,1,0
         *GET,KI_1,KCALC,0,K,1
         *VFILL,KTable(q4,1),data,KI_1

         midnode=8202

         Pfade i zur Berechnung der K-Werte in den Ebenen
         *DO,q5,2,9,1
         PATH,K_i,3,,48
         PPATH,1,(q5*1000)-999
         PPATH,2,midnode
         PPATH,3,(q5*1000)-838
         KCALC,0,1,1,0
         *GET,KI_i,KCALC,0,K,1
         *VFILL,KTable(q4,q5),data,KI_i
         midnode=midnode+32
         *ENDDO
         
         SET,NEXT

*ENDDO

PARSAV,ALL,C:\ANSYS\double_edge_crack\3D\elastic\K_Standard_plane_strain,txt,




!-----------------------------------------------------------------------
!-----------------------------------J-Integral--------------------------
!--------------------------------- für 6 Radien ------------------------
!-----------------------------------------------------------------------

/Post1
*DEL,q1
SET,First

!-----Schleife zur Bestimmung des J-Integralverlaufs
*DIM,JIntTable,ARRAY,50,2,6,,,,
*DIM,radius,ARRAY,6,1,,,,,
radius(1) = a/(8*4)
radius(2) = a/8
radius(3) = a/6
radius(4) = a/4
radius(5) = a/3
radius(6) = a/2

*DO,q1,1,N,1

M=1

!Unterschiedliche Radien r1 bis r6 um die Rissspitze
*DO,J,1,6,1

!Laden des Volumens und der Dehnungsenergie
ETABLE,VOL1,VOLU,!Fills a table of element values for further processing
                 !Vol1 = Any unique user defined label for use in
                 !       subsequent commands and output headings
                 !VOLU = Element volume. Based on unit thickness for
                 !       2-D plane elements (unless the thickness option
                 !       is used) and on the full 360 degrees for 2-D
                 !       axisymmetric elements.

ETABLE,STRAEN,SENE,!SENE = Stiffness" energy or thermal heat dissipation
                   !(applies to all elements where meaningful).

!Berechnung der Dehnungsenergiedichte
SEXP,STRENDEN,STRAEN,VOL1,1,-1, !Forms an element table item by
                                !         exponentiating and multiplying
                                !STRENDEN = Label assigned to results
                                !STRAEN = First labeled result item in
                                !         operation (Lab1)
                                !VOL1 = Second labeled result item in
                                !       operation (Lab2)
                                !1 = Exponent applied to STRAEN
                                !-1 = Exponent applied to VOLU1
                                !LabR = (|Lab1|^EXP1) x (|Lab2|^EXP2)

!-----Erstellung der Pfade für die Integration

*afun,deg

L=1

PATH,J2,29,30,20  !J2 = Name for this path
                  !29 = The number of points used to define this path
                  !30 = The number of sets of data which you can map to
                  !     this path.
                  !20 = The number of divisions between adjacent points.

*DO,I,0,180,180/28

       AX=radius(J)*COS(I)
       BY=radius(J)*SIN(I)
       PPATH,L,,AX,BY,0
       L=L+1

*ENDDO

!Mappen der Dehnungsenergiedichte auf dem Pfad

PDEF,MSTRENDE,ETAB,STRENDEN,AVG   !Interpolates an item onto a path.
                                  !MSTRENDE = Label assigned to the
                                  !           resulting path item
                                  !ETAB = Label identifying the item for
                                  !       interpolation (
                                  !STRENDEN = Component of the item
                                  !           (Dehnungsenergiedichte)
                                  !AVG = Option to average across
                                  !      element boundaries

!Integration der Dehnungsenergiedichte über y

PCALC,INTG,INT1,MSTRENDE,YG,1,  !PCALC,INTG,LabR,Lab1,Lab2,FACT1
                                !INTG = Finds an integral.
                                !Int1 = Label assigned to the resulting
                                !       path item.
                                !MSTRENDE = First labeled path item in
                                !           operation
                                !YG = Second labeled path item in
                                !     operation
                                !1 = Factor applied to Lab1
                                !Formula LabR = FACT1 x Int(Lab1 d(Lab2)

!Zuweisen der Ergebnisse zu einem Parameter

*GET,Term1,PATH,,LAST,INT1  !Term1 = The name of the resulting parameter
                            !PATH = Entity keyword.
                            !, = The number or label for the entity. In
                            !    some cases, a zero (or blank) ENTNUM
                            !    represents all entities of the set.
                            !LAST = The name of a particular item for
                            !       the given entity. Last value of path
                            !       item Lab from last path operation
                            !INT1 = The number (or label) for the
                            !       specified Item1

!Mappen der Spannungen sx,sy undsxy auf den Pfad

PDEF,stressx,S,X,AVG!Interpolates an item onto a path.
                    !stressx = Label assigned to the resulting path item
                    !S = Component stress
                    !X = Component stress
                    !AVG = Option to average across element boundaries

PDEF,stressy,S,Y,AVG
PDEF,stressxy,S,XY,AVG

!Definition des Einheitsnomalenvektors zum Pfad

PVECT,NORM,normx,normy,normz   !Interpolates a set of items onto a path.
                               !NORM = Defines a unit normal vector at
                               !       each interpolation point in the
                               !       direction of the cross product of
                               !       the tangent to the path and the
                               !       active Z axis. Resulting vector
                               !       components are in the active
                               !       coordinate system
                               !normx = Label assigned to X-component of
                               !        the resulting vector

!Berechnung des Zugvektors tx

PCALC,MULT,TX1,STRESSX,NORMX,1
PCALC,MULT,TX2,STRESSXY,NORMY,1
PCALC,ADD,TX,TX1,TX2,1,1,,

!Berechnung des Zugvektors ty

PCALC,MULT,TY1,STRESSY,NORMY,1
PCALC,MULT,TY2,STRESSXY,NORMX,1
PCALC,ADD,TY,TY1,TY2,1,1,,

!Berechnung der Ableitung der Verschiebung

!Berechnung der totalen Pfadlänge

*GET,DX0,PATH,,LAST,S

*SET,DX,DX0/10
*SET,DX1,1/DX

!Verschiebung des Pfades um DX/2 in die negative X-Richtung + Mappen der
!Verschiebungen UX1 und UY1 auf den Pfad

PCALC,ADD,XG,XG,,,,-DX/2

PDEF,UX1,U,X,AVG      !Interpolates an item onto a path.
                      !UX1 = Label assigned to the resulting path item
                      !U = Label identifying the item for interpolation;
                      !    X structural displacement
                      !X = Component of the item
                      !AVG = Option to average across element boundaries
PDEF,UY1,U,Y,AVG

!Verschiebung des Pfades um + DX/2 in die positive X-Richtung + Mappen
!der Verschiebungen UX2 und UY2 auf den Pfad
PCALC,ADD,XG,XG,,,,DX

PDEF,UX2,U,X,AVG
PDEF,UY2,U,Y,AVG

!Verschiebung des Pfades auf die ursprüngliche Position
PCALC,ADD,XG,XG,,,,-DX/2

!Berechnung von (UX2-UX1)/DX und (UY2-UY1)/DX

PCALC,MULT,SUM1,UX2,TX,DX1,                       !Sum1 = UX2 x TX x DX1
PCALC,MULT,SUM2,UX1,TX,DX1,

PCALC,MULT,SUM3,UY1,TY,DX1,
PCALC,MULT,SUM4,UY2,TY,DX1,

PCALC,ADD,SUM5,SUM1,SUM2,1,-1,,     !SUM5 = (1 x SUM1) + (-1 x SUM2) + 0
PCALC,ADD,SUM6,SUM3,SUM4,-1,1,,
PCALC,ADD,SUM7,SUM5,SUM6,1,1,,

!Integration über s

PCALC,INTG,INT2,SUM7,S,1,             !Formula: INT2 = 1 x Int(SUM1 d(S)

!Zuweisen des Ergebnisses zu einem Parameter

*GET,Term2,PATH,,LAST,INT2

*SET,JINT,Term1-Term2
*SET,JIntTable(q1,2,M),JINT*2

*SET,JIntTable(M,1,1),radius(J)

M=M+1

!PDEF,CLEAR       !Clear all previous defined paths

*ENDDO

*DEL,M

SET,NEXT

*ENDDO

PARSAV,ALL,C:\ANSYS\double_edge_crack\3D\elastic\J_Standard_plane_strain,txt,