#include <math.h>
void modl13(double d[],
double eps[6], double deps[6], double dt,
double hold[6], double hnew[6],
double sig[6], double ds[6][6]) {
/*
c-------------------------------------------------------------------
c * * F E A P *
* A Finite Element Analysis Program
c Copyright (c)
1984-2000: Robert L. Taylor
c
c C Language Example
Element for FEAP
c Copyright (c)
2000 : Sanjay Govindjee
c
c Simple viscoelastic
model
c with deviatoric
relaxation in one exponential mode
c uses mid-point
rule with recursion for history integral evaluation
c
c Input:
c d[0]
- lam
c d[1]
- mu
c d[2]
- relaxation time
c d[3]
- viscosity strength
c eps[6]
- strains
c deps[6]
- strain increments
c hold[6]
- old history
c dt
- time step
c
c Output:
c hnew[6]
- new history
c sig[6]
- stresses
c ds[6][6]
- tangent
c-------------------------------------------------------------------
*/
double lam,mu,tre,trs,p,muv,tau,k;
double trde,xfact;
int i,j;
tre = eps[0]+ eps[1]+ eps[2];
trde = deps[0]+deps[1]+deps[2];
lam = d[0];
mu = d[1];
tau = d[2];
muv = d[3];
xfact = exp(-dt/(2.0*tau));
for(i=0;i<6;i++) {
for(j=0;j<6;j++)
ds[i][j] = 0.0;
}
/* Elastic volumetric
response */
k = 2.0*mu/3.0 + lam;
p = k*tre;
sig[0] = p;
sig[1] = p;
sig[2] = p;
sig[3] = 0.0;
sig[4] = 0.0;
sig[5] = 0.0;
/* Viscous Deviatoric
Response */
for(i=0;i<3;i++) {
eps[i] = eps[i]
- tre/3.0;
deps[i] = deps[i] - trde/3.0;
}
/* Compute and
update new history */
for(i=0;i<3;i++) {
hnew[i] = xfact*xfact*hold[i]
+ xfact*deps[i];
sig[i] = sig[i]
+ 2.0*mu*eps[i] + 2.0*muv*hnew[i];
}
for(i=3;i<6;i++) {
hnew[i] = hold[i]
+ xfact*deps[i];
sig[i] = sig[i]
+ mu*eps[i] + muv*hnew[i];
}
/* Compute consistent
tangent */
xfact = muv*xfact + mu;
for(i=0; i<3; i++){
for(j=0; j<3; j++) ds[i][j]
= k-xfact*2.0/3.0;
}
ds[0][0] += 2.0*xfact;
ds[1][1] += 2.0*xfact;
ds[2][2] += 2.0*xfact;
ds[3][3] += xfact;
ds[4][4] += xfact;
ds[5][5] += xfact;
}