Appendix B
#pragma rtGlobals=1 // Use modern global access method.
Function diff(L)
Variable L
Variable Jm, arat, lrat, ts, xp, dx, dt, csrat, modrat, phi
// Physical constants for the problem
Jm = 97.2 // max. value of first strain invariant related to max extension by J =
// lambda^2+2/lambda-3
arat = 3 // ratio of cross-sectional area of fibril to that of foot
lrat = 0.3333333 // ratio of intial length of foot to that of fibril
phi = 0.001 // ratio of diffusion time scale to viscoelastic relaxation time scale
modrat = 5 // ratio of shell modulus to core modulus
csrat = 0.01 // ratio of shell thickness to core thickness
Make /N=17 /O omeg // non-dimensional pulling velocity = V/L_fib*(L^(m+1)/D)^(1/m)
// range of omega values to explore
omeg[] = {0.01,0.025,0.05,0.075,0.1,0.25,0.5,0.75,1,2.5,5,7.5,10,25,50,75,100}
// Computational range and resolution
ts = 80000 // time steps
xp = 51 // number of points in x direction
Variable i, j, k, lam, bc, m, jmin, g, ct, ctmin, sq, m2, t, lamdot
Make /N=10 /O gridsL, dxs, dts
Make /N=(10,xp) /O xval
Make /N=(xp) /O utemp
Make /N=17 /O stopind, lamnommax
Make /N=(17,10) /O rmpt
Make /N=(17,ts) /O sig, stretch, stretchnom, lzero, ftlength, ttot, ftstravg
Make /N=(17,ts) /O lzerodot, txx, tzz, sigvl, siggent
Make /N=(17,xp,ts) /O uxt
// set up the waves of x values
// first specify lengths where remeshing occurs
gridsL[] = {1, 0.6, 0.36, 0.216, 0.1296, 0.07776, 0.046656, 0.0279936, 0.01679616, 0.010077696}
i = 0
do
dxs[i] = gridsL[i]/(xp-1) // spacial resolution
dts[i] = 0.9*dxs[i]*dxs[i]/2 // temporal resolution; for stability, we need dt < dx*dx/2
j = 0
do
xval[i][j] = 1-gridsL[i] + j*dxs[i]; j = j+1
while (j < xp)
i = i + 1
while(i < 10)
// Main loop to solve the governing equations as a function of Omega
k = 0
do
uxt[k][][0]= 0 // initial condition
// initial values of time dependent variables of interest
sig[k][0] = 0; stretch[k][0] = 1; stretchnom[k][0] = 1; lzero[k][0] = 0; lzerodot[k][0] = 0;
ttot[k][0] = 0; ftlength[k][0] = 1; txx[k][0] = 0; tzz[k][0] = 0; sigvl[k][0] = 0; siggent[k][0] = 0
i = 0; jmin = 0; g = 0; t = 0
dx = dxs[g]; dt = dts[g]
do
// current time
t = t + dt
ttot[k][i+1] = t
// actual overall stretch of fibril, nominal stretch, strain rate
stretch[k][i+1] = (1 + omeg[k]*t)/(1+lzero[k][i])
stretchnom[k][i+1] = 1 + omeg[k]*t
lam = stretch[k][i+1]
lamdot = omeg[k]/(1+lzero[k][i]) - lzerodot[k][i]*stretch[k][i]/(1+lzero[k][i])
// calculation of viscoelastic, elastic components as well as total stress in fibril
txx[k][i+1] = txx[k][i] - dt*(lamdot*(txx[k][i]+1/3)/lam+phi*txx[k][i])
tzz[k][i+1] = tzz[k][i] + dt*(2*lamdot*(tzz[k][i]+1/3)/lam-phi*tzz[k][i])
sigvl[k][i+1] = -txx[k][i+1] + tzz[k][i+1]
siggent[k][i+1] = (lam-1/lam/lam)/(1-(lam*lam+2/lam-3)/Jm)/3
sig[k][i+1] = ((1-csrat)^2*sigvl[k][i+1]+(1-(1-csrat)^2)*modrat*siggent[k][i+1])/((1-csrat)^2+(1-(1-csrat)^2)*modrat)
// solve for displacement in foot as function of x at current time step
j = jmin + 1
do
uxt[k][j][i+1] = uxt[k][j][i] + dt/dx/dx*(uxt[k][j+1][i]-2*uxt[k][j][i]+uxt[k][j-1][i])
j = j +1
while (j < xp)
// assign boundary conditions
bc = sig[k][i+1]*arat // stress at boundary with fibril
uxt[k][jmin][i+1] = uxt[k][jmin+2][i+1]-bc*2*dx // left bc
uxt[k][xp-1][i+1] = uxt[k][xp-3][i+1] // right bc
// finding lzero (right hand boundary position) & remeshing if necessary
j = jmin
do
if (abs(uxt[k][j][i+1]) < xval[g][j])
m = (abs(uxt[k][j][i+1])-abs(uxt[k][j-1][i+1]))/(xval[g][j]-xval[g][j-1])
// unstressed length of transferred material for next step
lzero[k][i+1]= (abs(uxt[k][j][i+1]) - m*xval[g][j])/(1-m)*lrat/arat
lzerodot[k][i+1] = (lzero[k][i+1]-lzero[k][i])/dt
jmin = j-1
// remeshing, uses a linear interpolation between grid points to get proper values for new grid
if (jmin == 20)
if (g<9)
printf "remeshing - Omega = %g, Mesh Number = %d, Time Step = %d, Time = %g\r", omeg[k], g+1, i+1, ttot[k][i+1]
rmpt[k][g+1] = i + 1; g = g + 1; ctmin = 0; sq = 0
jmin = 0; dx = dxs[g]; dt = dts[g]
do
ct = 0
do
if (xval[g][sq]<=xval[g-1][ct])
m2 = (uxt[k][ct][i+1]-uxt[k][ct-1][i+1])/(xval[g-1][ct]-xval[g-1][ct-1])
utemp[sq] = uxt[k][ct-1][i+1] + m2*(xval[g][sq]-xval[g-1][ct-1])
ctmin = ct-1
break
endif
ct = ct + 1
while(ct<xp)
sq = sq+1
while(sq<xp-1)
utemp[xp-1] = uxt[k][xp-1][i+1]
ct = 0
do
uxt[k][ct][i+1] = utemp[ct]
ct = ct+1
while(ct<xp)
endif
endif
break
endif
j = j +1
while (j < xp)
// failure criterion – when slip reaches 99% of original foot length
if (abs(uxt[k][xp-1][i+1]) > 0.99)
printf "Omega = %g, Critical Time Step = %d\r", omeg[k], i
stopind[k]=i
lamnommax[k] = 1 + omeg[k]*t -1
break
endif
ftlength[k][i+1] = 1 - abs(uxt[k][xp-1][i+1])
ftstravg[k][i+1] = (lzero[k][i+1]/lrat*arat-abs(uxt[k][xp-1][i+1]))/(1-lzero[k][i+1]/lrat*arat)
i = i + 1
while (i < ts)
k = k+1
while (k < 17)