Material parameter estimation using GPL Maxima

Maxima is a Computer Alegbra System which is freely available under the Gnu Public License. Its home page is http://maxima.sourceforge.net .

For many computer algebra tasks, including material law parameter estimation, it is a viable alternative to Mathematica.

There are several frontends available but I recommend wxMaxima http://wxmaxima.sourceforge.net .

Here is my Maxima script for estimating some parameters for the Costa law. It is not optimal but gives an idea of what can be done.

{{{ /* Optimise material parameters for non-linear finite strain elasticity using Levenberg-Marquardt algorithm. Glenn Ramsey November 2005. Based on Holger Schmid's Mathematica implementation. NOTE Requires maxima > 5.9.2, currently (24/11/05) that means using the cvs version. */

kill(all);

/* set up some matrices to represent the deformation gradient tensor */ Fnf: matrix(1,k,0, 0,1,0, 0,0,1)$ Ffn: transpose(Fnf)$ Fns: matrix(1,0,0, 0,1,0, 0,k,1)$ Fsn: transpose(Fns)$ Ffs: matrix(1,0,0, 0,1,0, k,0,1)$ Fsf: transpose(Ffs)$ Fss: matrix(1,0,0,0,1,0,0,0,1+k)$

Nsf: matrix(0,0,1)$ Nsn: matrix(0,0,1)$ Nss: matrix(0,0,1)$

/* F, NA are global lists of matrices / / Put the matrices in a list so we can refer to them by index */ F:Fsf, Fsn, Fss; NA:Nsf, Nsn, Nss;

dim:3$ (f:1, n:2, s:3)$ maxstrain:0.25$

/* Use named variables for indices so we know what mode we are dealing with. / / Data is available for 2 shear modes and 1 normal mode so here I'm setting things up so that the modes we don't have data for are not used. Doing it this way makes it easier if I change the material law, then I only have to rearrange the index numbers in one place. */ nparam:4; (ff:nparam+1, fn:nparam+1, fs:1, nn:nparam+1, ns:2, ss:3)$ (aa:4)$

/* Strain energy function - Costa et al, Ph.Tr.R.Soc.Lond.,359,2001 / / Q(b):=bff,1*(Ef,f)^2 + 2*bfn,1*((1/2)*(En,f+Ef,n))^2 + 2*bfs,1*((1/2)*(Ef,s+Es,f))^2 + bnn,1*(En,n)^2 + 2*bns,1*((1/2)*(En,s+Es,n))^2 + bss,1*(Es,s)^2 ;

W(b):=(1/2)*ba,1*(exp(Q(b))-1); */

/* The above form doesn't work, so remove unused terms and subst for Q(b) */

W(b):=(1/2)*baa,1*(exp( 2*bfs,1*((1/2)*(Ef,s+Es,f))^2 + 2*bns,1*((1/2)*(En,s+Es,n))^2 + bss,1*(Es,s)^2 )-1);

/* 2nd Piola-Kirchhoff stress / / derivative of strain energy function wrt components of strain tensor */ S:zeromatrix(dim,dim)$ for i thru dim do for j thru dim do S[i,j]:diff(W(b),Ei,j)$

/* function for t (model force) in terms of F, N, x */ tm(v,F,N,x):=block(b,k,E,b:v,k:x,Fev:ev(F),E:(1/2)*(transpose(Fev).Fev-ident(dim)),ev(Fev.S.N));

/* Experimental force is given by polynomials fitted to the data These cubics do not fit the data well for small strains, but this is what was published... Douglas et el,JEB,201,1998. */

fpd(x):=4970*x^3 + 6320*x^2 - 425*x + 26$ flm(x):=-1090*x^3 + 3800*x^2 - 300*x +21$ fnn(x):=1.39e5*x^3 + 230*x^2 - 522*x + 28$

Te:fpd,flm,fnn;

te(x,mode):=(Temode)(x);

ta(x,mode,v):=tm(v,Fmode,NAmode,x)mode,1;

plot2d((te(x,1)),(te(x,2)),(te(x,3)),x,0,maxstrain,y,0,200, gnuplot_term,ps, gnuplot_out_file, "experimental.ps");

/* Generate a list of strain steps for the objective function. */ u:block(y,y:[], for n:0 thru maxstrain step maxstrain/10 do y:append(y,n),y)$ uln:length(u)$

/* objective function in terms of material parameters*/

/* display the objective function in noun form */ 'sum('sum((ta(un,mode,v)-te(un,mode))^2,n,1,uln),mode,1,3);

/* now actually define it */ forcediff(mode,v):=sum((ta(un,mode,v)-te(un,mode))^2,n,1,uln); obj(v):=sum(forcediff(mode,v),mode,1,3);

/* have to do this substitute and evaluate thing there, not sure why */ m:obj(b)$ mobj(v):=block(b:v,(ev(m)))$

/* partial derivatives of the objective function wrt material parameters */ g:zeromatrix(nparam,1)$ for i:1 thru nparam do g[i,1]:diff(obj(b),bi,1); grad(v):=block(b:v,r:zeromatrix(nparam,1),for i:1 thru nparam do (r[i,1]:ev(gi,1)),r);

/* 2nd partial derivatives of the objective function wrt material parameters */ h:zeromatrix(nparam,nparam)$ for i:1 thru nparam do for j:1 thru nparam do h[i,j]:diff(gj,1,bi,1); hess(v):=block(b:v,r:zeromatrix(nparam,nparam),

for i:1 thru nparam do for j:1 thru nparam do (r[i,j]:ev(hi,j)),r);

/* Now the equations are set up */

/* Optimization loop Levenberg-Marquardt*/

/* Initial guess / / theta[0]:transpose(matrix(20,20,20,20)); */ theta[0]:matrix(17,9,17,90); lmda:1; aom:1; ath:1;

for i:0 while (aom > 10^-5 or ath > 10^-5) do (

omega[0]:float(mobj(theta0)), delta:-invert(float(hess(theta0)+lmda*ident(nparam))).float(grad(theta0)), display(delta), theta[1]:theta0 + delta, omega[1]:float(mobj(theta1)), display(omega1), display(omega0), display(theta1), if (omega1<omega0) then (

ath:(theta1-theta0).(theta1-theta0), aom:abs(omega1-omega0), theta[0]:theta1, display(ath), display(omega1), display(theta1), lmda:lmda/10 ) else ( lmda:lmda*10 ),

i:i+1, display(i)

);

/* r:zeromatrix(nparam+1,1); for i thru nparam do r[i,1]:theta1i,1; */

r:theta1;

/* display some results */ te(0.1,1); float(ta(0.1,1,r)); te(0.1,2); float(ta(0.1,2,r)); te(0.1,3); float(ta(0.1,3,r)); plot2d((te(x,1)),ta(x,1,r),(te(x,2)),ta(x,2,r),(te(x,3)),ta(x,3,r),x,0,maxstrain,y,0,2000, gnuplot_term,ps, gnuplot_out_file, "fitted.ps"); }}}