Dear Jörgen,
I want to share a code with the polymer modeling community and figured your website is probably a good medium to try that out. Please review my message below and share if you agree.
Regards,
Travis
In 2017, I ran into an issue that required me to rely on crosshead (ie. grip) displacement to compute strain in rubber dumbbells in uniaxial tension. To prevent the grips from fracturing the specimen at very low strain (a big problem at high temperatures), I had to grip the dumbbell at a section that was much wider than the gage section. Of course, this results in large errors in the strain calculation. To circumvent the problem, I came up with a method to correct stress-strain data when a dumbbell is gripped at a wide section and direct measurement of strain is not possible. Since then, a surprising number of people have come to me with the same problem. Therefore, I am sharing my code here. The code is, hopefully, commented sufficiently well for general use, but I will follow up later with a full publication on the concept.
Disclaimer: Nobody has peer-reviewed the code for accuracy, but it has worked for 4 case studies. Use discretion when using the code. The author appreciates any feedback on bugs or issues encountered with the code.
PROGRAM nBlock10 ! v10.0
!
! *****************************************************************
! Author: Travis Hohenberger
! Date: 2019-11-04
! E-mail: twhohen@gmail.com
! System: Windows 10
! Software: Netbeans IDE Dev (201804200002)
! Java 1.8.0_101; Java Hotspot 64 Bit Server VM 25.101-b13
! Java Runtime Environment 1.8.0_101-b13
! **************************************************************
!
! Dumbbells are often stretched in uniaxial extension to collect stre-
! ss-strain data. When a dumbbell is gripped at a section wider than
! its gage section, and strain is computed from displacement of the g-
! rips, the calculated strain is lower than the actual strain in the
! gage section. This program corrects stress-strain data when direct
! measurement of strain (for instance, with a video extensometer) in
! the gage section is not possible, and therefore must be estimated f-
! rom grip (ie. crosshead) displacement. The program decomposes a dum-
! bbell into a series of blocks as illustrated in SCHEMATIC 1. As cur-
! rently written, the code creates 200 'n-blocks' between the upper (
! ie. tab) and lower (ie. gage) sections.
! ____________________
! &&&&&&&&&&&&&&&&&& | |
! &&&&&&&&&&&&&&&&&& | |
! &&&&&&&&&&&&&&&&&& | |
! &&&&&&&&&&&&&&&&&& | UPPER |
! &&&&&&&&&&&&&&&&&& ==> Tab section | BLOCK |
! &&&&&& TAB &&&&&&& | index = 202 |
! &&&&&&&&&&&&&&&&&& | in code |
! &&&&&&&&&&&&&&&&&& |__________________|
! &&&&&&&&&&&&&&&&& |_________________|
! &&&&&&&&&&&&&&&& ==> Smooth radius |________________|
! &&&&&&&&&&&&&&& |_______________|
! &&&&&&&&&&&&& ==> Tangent of radii |_____________| Seven
! &&&&&&&&&&&& |___________| ==> n-Block
! &&&&&&&&&&& ==> Smooth radius |_________| transitions,
! &&&&&&&&& |________| index = 2, 3, .
! &&&&&&&& | | .., 201 in code
! &&&&&&&& | |
! &&&&&&&& | |
! &&&&&&&& | LOWER |
! & GAGE & | BLOCK |
! &&&&&&&& | |
! &&&&&&&& ==> Gage section | index |
! &&&&&&&& | = 1 |
! &&&&&&&& | in |
! &&&&&&&& | code |
! &&&&&&&& |_______|
!
! (a) (b)
!
! SCHEMATIC 1: (a) Quarter-dumbbell with "smooth" geometry;
! (b) dumbbell decomposed into "block" geometry
!
! *****************************************************************
!
! The code uses the incompressible form of the generalized Yeoh strai-
! n-energy function:
!
! W = K1*(I1-3)^m + K2*(I1-3)^p + K3*(I1-3)^q
!
! which can be replaced by the user by changing the force equation.
!
! Reference:
!
! Hohenberger et al (2019). A constitutive model for both low and hig-
! h strain nonlinearities in highly filled elastomers and implementat-
! ion with user-defined material subroutines in Abaqus. Rubber Chemis-
! try and Technology, 92 (4), 653-686.
!
! Disclaimer: This code is not computationally optimal, but it has wo-
! ked for multiple cases tested by the author. More detai-
! ls will be shared in a journal publication.
!
! *****************************************************************
!
IMPLICIT NONE
!
REAL K1 , K2 , K3 , deL(202), w(202) , t , Delta ,&
Fo , Lam0 , F(202), Lam(202), ErrF(202), delNew, Fout ,&
deLtot, ErrDel, R1 , R2 , w0 , wTab , SigOut,&
alpha , h1 , h2 , L0 , dh1 , dh2 , y(202),&
Lgrip , Lbell , em , pe , qu , MxErrF, DelEnd,&
SigStore(200) , LamStore(200) , LamFact , ForcFac
!
INTEGER i, j, k, int1, int2
!
! ************************************************************
! - MATERIAL MODEL PARAMETERS FOR GEN-YEOH STRAIN-ENERGY FUNCTION --
! ************************************************************
!
K1 = 5.38 ! MPa
K2 = -2.85 ! MPa
K3 = 0.40 ! MPa
em = 0.89
pe = 1.08
qu = 1.85
!
! ************************************************************
! ------ CONVERGENCE FACTORS FOR STRETCH AND FORCE ITERATIONS ------
! ************************************************************
!
LamFact = 10. ! Increasing values speed up convergence but can
ForcFac = 0.1 ! result in unacceptable error. Monitor 'ErrDel'
! ! and 'MaxErrF' in the output window to ensure
! ! they remain << 1.
!
! ************************************************************
! ----------- DUMBBELL DIMENSIONS (mm), SEE SCHEMATIC 2 ------------
! ************************************************************
!
Lbell = 57.5 ! Dumbbell length (half-length, die C = 57.5)
Lgrip = 27.899 ! Grip separation (half-length)
L0 = 16.5 ! Gage length (half-length, die C = 16.5)
w0 = 3.12305 ! Width of gage section (half-width)
wTab = 12.5 ! Width of tab (half-width, die C = 12.5)
t = 2.3335 ! Dumbbell thickness (not illustrated)
R1 = 14. ! Dumbbell radius at gage (die C = 14)
R2 = 25. ! Dumbbell radius at tab (die C = 25)
Delta = 0.078 ! Initial displacement for starting iterations
DelEnd = 32. ! Target end displacement (half-displacement)
!
! | |
! |<--- wTab ---->| Delta, DelEnd
! | | ^
! -------------------- &&&&&&&&&&&&&&&&& |
! ^ &&&&&&&&&&&&&&&&& |
! | &&&&&&&&&&&&&&&&& |
! | &&&&&&&&&&&&&&&&& |
! | &&&&&&&&&&&&&&&&& |
! | ----------- &*********& ------------
! | ^ &&&&&&&&&&&&&&&&&
! | | &&&&&&&&&&&&&&&&&
! | | &&&&&&&&&&&&&&&& <---- R2
! | | &&&&&&&&&&&&&&&
! | | &&&&&&&&&&&&&
! Lbell | &&&&&&&&&&&&
! | | &&&&&&&&&& <---------- R1
! | | &&&&&&&&&
! | Lgrip &&&&&&&& -----------
! | | &&&&&&&& ^
! | | &&&&&&&& |
! | | &&&&&&&& |
! | | &&&&&&&& |
! | | --->& w0 &<--- L0
! | | &&&&&&&& |
! | | &&&&&&&& |
! v v &&&&&&&& v
! -------------------- &&&&&&&& -----------
!
! SCHEMATIC 2: Illustration of dimensions to be input by user
!
! ************************************************************
! ------------------- INITIALIZE SOME VARIABLES --------------------
! ************************************************************
!
! NOTE: The guess values for 'Lam0' and 'Fo' affect convergence and
! in addition to 'LamFact' & 'ForcFac', may result in unaccep-
! table error if they are too large. If large errors are enco-
! untered even with very small 'LamFact' & 'ForcFac', try red-
! ucing the initial guess values for 'Lam0' and 'Fo'.
!
Lam0 = 1.000001 ! Initial guess value for stretch
Fo = 0.1 ! Initial load increment (units = Newtons)
Lam = Lam0 ! Initialize all Lam(i)'s to Lam0
F = 0. ! Initialize pulling force to 0
ErrF = 500. ! Initialize force error to a large value
!
! ************************************************************
! ---------- CALCULATION OF HEIGHTS & WIDTHS OF n-BLOCKS -----------
! ***************************************************************
int1 = 100 ! Number of intervals in which to decompose dumbbell
int2 = 100 ! radius transitions. This code must be updated if
! '100' is not specified for 'int1' or 'int2'. When
! ! testing the code, 100 intervals was found to give
! ! less than 0.1% discretization error.
!
! The geometric calculations below do not require modification if t-
! he default dimensioning scheme for an ASTM dumbbell applies.
!
alpha = ACOS((R1+R2+w0-wTab)/(R1+R2)) ! Angle b/t horizontal and
! the line connecting the
! centers of R1 & R2.
!
h1 = L0 + R1*SIN(alpha) ! Height to tangent of R1 and R2
h2 = h1 + R2*SIN(alpha) ! Height to transition b/t R2 & wTab
dh1 = (h1-L0) / int1 ! n-block interval height in R1 zone
dh2 = (h2-h1) / int2 ! n-block interval height in R2 zone
y(1) = L0 ! Set y(1) equal to gage length
w(1) = w0 ! Set w(1) equal to gage width
w(202) = wTab ! Set w(202) equal to tab width
!
DO i=2,101 ! Loop to calculate the height to
IF (i.EQ.2) THEN ! the n-blocks and respective wid-
y(i) = y(i-1) + dh1/2. ! ths in the R1 transition zone.
ELSE
y(i) = y(i-1) + dh1
END IF
w(i) = -((R1**2 - (y(i)-L0)**2))**0.5 + w0 + R1
END DO
!
DO i=102,201 ! Loop to calculate the height to
IF (i.EQ.102) THEN ! the n-blocks and respective wid-
y(i) = h1 + dh2/2. ! ths in the R2 transition zone.
ELSE
y(i) = y(i-1)+dh2
END IF
w(i) = (R2**2 - (y(i) - (L0 + (R1+R2)*SIN(alpha)))**2)**0.5&
+ wTab - R2
END DO
!
IF (Lgrip.GT.h2) THEN ! Compute the height of the upper tab se-
y(202) = Lgrip - h2 ! ction. Neglect this section if the grip
ELSE ! location is below the upper tab region.
y(202) = 0.
END IF
!
! ***************************************************************
! --------- ITERATIVE CALCULATIONS FOR FORCE AND STRETCH -----------
! ***************************************************************
!
k = 1 ! Index k is used for data storage / output.
!
220 j = 1 ! Restart location if the the current target displ-
! ! acement, 'Delta', has not been reached after for-
! ! ce and stretch iterations.
!
210 DO i = j,202 ! Loop for force calculations in each n-block.
!
! Force from generalized Yeoh strain-energy function:
!
F(i) = 2*w(i)*t*(Lam(i)-1/(Lam(i)**2))*&
(em * K1 * (Lam(i)**2 + 2/Lam(i) - 3)**(em-1.) + &
pe * K2 * (Lam(i)**2 + 2/Lam(i) - 3)**(pe-1.) + &
qu * K3 * (Lam(i)**2 + 2/Lam(i) - 3)**(qu-1.))
!
ErrF(i) = Fo - F(i) ! Residual force error
!
IF (ErrF(i).GT.500.0) THEN ! Convergence check
WRITE(*,*) 'ErrF = ' , ErrF(i)
WRITE(*,*)
WRITE(*,*) 'Force error > 500. Try adjusting Fo.'
STOP
END IF
!
IF (ErrF(i).GT.0.) THEN ! If error in force
IF (ErrF(i).LT.0.05) THEN ! is positive, incr-
Lam(i) = Lam(i) + 0.000001*LamFact ! ement stretch and
ELSE IF (ErrF(i).LT.0.5) THEN ! use the GOTO stat-
Lam(i) = Lam(i) + 0.00001*LamFact ! ement to recalcul-
ELSE IF (ErrF(i).LT.5.) THEN ! ate the force unt-
Lam(i) = Lam(i) + 0.0001*LamFact ! il the force resi-
ELSE IF (ErrF(i).LT.50.) THEN ! dual, 'ErrF(i)',
Lam(i) = Lam(i) + 0.001*LamFact ! becomes negative.
ELSE IF (ErrF(i).LT.500.) THEN
Lam(i) = Lam(i) + 0.01*LamFact
END IF
GO TO 210
END IF
!
j = j+1 ! Increment j so that later, if necessary, the outer
END DO ! loop will restart until the final target displacem-
! ! ent, 'DelEnd', is reached or exceeded.
!
deL(1) = L0*(Lam(1)-1) ! Displacement of gage section
deL(202) = y(202)*(Lam(202)-1) ! Displacement of tab section
!
DO i=2,101 ! Displacement of each n-block
deL(i) = dh1*(Lam(i)-1) ! in R1 section
END DO
!
DO i=102,201 ! Displacement of each n-block
deL(i) = dh2*(Lam(i)-1) ! in R2 section
END DO
!
DO i = 2,201 ! Set each n-block displacement
IF (y(i).GT.Lgrip) THEN ! to 0 if the n-block is above
deL(i) = 0. ! the grip location.
END IF
END DO
!
deLtot = 0. ! Initialize total displacement.
!
DO i=1,202 ! Compute total displacement.
deLtot = deLtot + deL(i)
END DO
!
ErrDel = Delta - deLtot ! Compare the sum of all n-block
! ! displacements to the target d-
IF (ErrDel.GT.0.) THEN ! isplacement, 'Delta'. If error
IF (ErrDel.LT.0.1) THEN ! is positive, increment force &
Fo = Fo + 0.01*ForcFac ! restart the force and stretch
ELSE IF (ErrDel.LT.1.) THEN ! iteration loops with the GOTO
Fo = Fo + 0.1*ForcFac ! statement.
ELSE IF (ErrDel.LT.10.) THEN
Fo = Fo + 1.0*ForcFac
ELSE IF (ErrDel.LT.100.) THEN
Fo = Fo + 10.*ForcFac
ELSE IF (ErrDel.LT.1000.) THEN
Fo = Fo + 100.*ForcFac
END IF
GO TO 220
END IF
!
! ***************************************************************
! --------- END OF ITERATION LOOPS FOR FORCE AND STRETCH -----------
! ***************************************************************
!
DO i = 1,202 ! Store the maximum force er-
IF (ErrF(i+1).LT.ErrF(i)) THEN ! ror of all the n-blocks.
MxErrF = ErrF(i+1)
ELSE
MxErrF = ErrF(i)
END IF
END DO
!
Fout = F(1) ! Store the force in the gage section.
SigOut = F(1)/w(1)/t ! Store the stress in the gage section.
!
LamStore(k) = Lam(1) ! Store stress and stretch as output that
SigStore(k) = SigOut ! the user can copy from output window.
!
! This next sequence of commands may seem confusing. The desired st-
! ress & stretch outputs are stored in arrays LamStore(k) & SigStor-
! e(k), but a stress corresponding to an n-block corrected stretch
! is stored one increment behind the current increment due to the m-
! anner in which this code increments displacement from 'Delta' to
! 'DelEnd'. Specifically, stress and stretch during the first itera-
! tion are not n-block corrected. The code forces the crosshead str-
! etch at iteration 'k+1' to equal the gage stretch at iteration 'k'
! , and then maps the updated gage stretch (at 'k+1') to the stress
! at 'k', creating an n-block corrected (stretch, stress) pair.
!
IF (k.EQ.1) THEN ! Write headings for stress & stretch outputs.
WRITE(*,*)
WRITE(*,*) ' Gage Stretch' , ' Nominal Stress / MPa',&
' ErrDel ' , ' MaxErrF '
WRITE(*,*) ' ------------ ' , '--------------------',&
' ------------ ' , '-------------'
ELSE
WRITE(*,360) LamStore(k),SigStore(k-1),& ! Write stress, str-
ErrDel,MxErrF ! etch, and errors.
END IF
!
k = k+1 ! Increment k for the next target di-
! ! displacement, 'Delta'.
!
delNew = Lgrip*(Lam(1)-1.) ! Set the next displacement incremen-
! ! t, 'Delta', so that the stretch ba-
IF (delNew.LT.delEnd) THEN ! sed on crosshead displacement in i-
Delta = delNew ! ncrement 'k+1' is equal to the str-
Fo = Fout ! etch in the gage section at increm-
GO TO 220 ! ent 'k'. If the total target displ-
END IF ! acement, 'DelEnd', has not been re-
! ! ached, conduct another iteration of
! ! the entire sequence.
WRITE(*,*)
WRITE(*,*) 'Final errors, force, displacement, & average stretch'
WRITE(*,*) '----------------------------------------------------'
WRITE(*,310) 'ErrDel = ' , ErrDel ,' mm Max. displacement error'
WRITE(*,320) 'MxErrF = ' , MxErrF ,'N' , 'Max. force error'
WRITE(*,330) 'F = ' , Fout ,'N' , 'Pulling force'
WRITE(*,340) 'Delta = ', Delta ,'mm','Crosshead displacement'
WRITE(*,350) 'LamT = ' , Delta/Lgrip+1.,'Crosshead (avg) stretch'
WRITE(*,*)
!
310 FORMAT(A10,F12.6,A29)
320 FORMAT(A10,F12.6,A3,A19)
330 FORMAT(A10,F12.3,A3,A16)
340 FORMAT(A10,F12.3,A3,A25)
350 FORMAT(A10,F12.4,A29)
360 FORMAT(F12.4,F18.3,F26.6,F17.6)
!
END PROGRAM nBlock10
!
! ******************************************************************
! --------------------- USER GUIDE / INSTRUCTIONS ---------------------
! ******************************************************************
!
! 1. Plot uniaxial tension stress-strain data with strain based on grip
! (ie. crosshead) displacement.
!
! 2. Fit a strain-energy function to the stress-strain data and record
! model parameters.
!
! 3. Input model parameters in the n-block code, and update the force
! equation if the generalized Yeoh strain-energy function is not us-
! ed.
!
! 4. Run the program to generate corrected stress-strain data for the
! material. If more data points are desired, adjust inputs 'Delta'
! and 'DelEnd' to result in different starting and ending points in
! the calculator. See SCHEMATIC 3 for an illustration of n-block co-
! rrected stress-strain data.
!
!
! | / /
! | Hyperelastic model _/ __/
! | fit to stress-strain ___ __/ ____/
! | data based on cross- | ___/ ____/
! | head displacement | ___/ ____/
! | __V___/ _______/
! S | ___/ _____/ ^
! T | __/ ___/ |
! R | __/ ___/ |__ n-block corrected
! E | __/ ___/ stress-strain output
! S | _/ __/
! S | _/ __/
! | _/ _/
! | / /
! | / _/
! | / _/
! | / /
! | /_/
! |/_______________________________________________________________
! STRAIN
!
! SCHEMATIC 3: Illustration of n-block corrected stress-strain data
! which, via a hyperelastic material model, shifts st-
! ress-strain data based on grip (ie. crosshead) disp-
! lacement.
