Code covered by the BSD License  

Highlights from
Constitutive Models of Linear Viscoelasticity with Laplace Transform

Constitutive Models of Linear Viscoelasticity with Laplace Transform

by

 

The program computes compliance and relaxation function.

viscoelasticity(B)
function [J,R] = viscoelasticity(B)
%function that calculate stress relaxation function and creep
%% input
syms s
%% MATRIX SETTINGS_V2
[N,n]=size(B);
B1=double(B(:,[1,2]));
P=max(max(B1));
%setting A1
A1=zeros(N,P-1);
for i=1:N
    A1(i,[B1(i,1):B1(i,2)-1])=-1;
end
%setting A2
A2=sym(zeros(N,N+2));
for i=1:N
    A2(i,i)=1/B(i,3)+1/(B(i,4)*s);
end
%setting A3
A3=zeros(P,P-1);
%setting A4
A4=zeros(P,N+2);
A4(1,:)=[-kronecker(ones(1,N),B1(:,1)'),1,0];
for i=2:P-1
    A4(i,1:N)=kronecker(i*ones(1,N),B1(:,2)')...
        -kronecker(i*ones(1,N),B1(:,1)');
end
A4(P,:)=[-kronecker(P*ones(1,N),B1(:,2)'),1,0];
%setting A5
A5=[-ones(1,P-1),zeros(1,N),0,1];
%setting A
A=[A1,A2;A3,A4;A5];
clear A1 A2 A3 A4 B1 i m n
%% Gauss elimination
[G]=gauss(A);
%% inverse Laplace transform
C=-G(P+N,P+N)/G(P+N,P+N+1);
if (isequal(class(C),'sym')==1)&&(isequal(findsym(C,1),'s')==1)
    syms t
    J=ilaplace(C/s,s,t);%compliance function
    R=ilaplace(1/(s*C),s,t);%relaxation function
    J2=ilaplace(C/s^2,s,t);
    R2=ilaplace(1/(s^2*C),s,t);
end
end

Contact us