% Berechnung der monotonen Arbeitslinie

% Zerstauchungsenergie Gcl
Gcl=30
% Elementgrösse
lE=500
% charakteristische Länge lch
lch=(lE^3/8)^(1/3)
% Druckfestigkeit
fcm=40
% Dehnung bei Erreichen von fcm
ec1=0.0022
% E-Modul
Ec=40000
% Zerstauchungsenergie pro Länge
gcl=Gcl/lch
% Neigung des Entlastungsastes im Druckbereich
bc=0.75
% Längenparameter lc
lc=Gcl/(fcm*(ec1*(1-bc)+bc*fcm/Ec))
% Steuerung der Kurvenneigung im Nachbruchbereich
gc=3.14^2*fcm*ec1/(2*(gcl-0.5*fcm*(ec1*(1-bc)+bc*fcm/Ec))^2)
% ideeller E-Modul
Eci=2/(3*Ec)*(fcm/ec1)^2-4/3*fcm/ec1+5/3*Ec

% Bereich 1 der Arbeitslinie
x1=0:0.0001:0.4*fcm/Ec
y1=Ec*x1
% Bereich 2 der Arbeitslinie
x2=0.4*fcm/Ec:0.0001:0.0022
y2=(Eci*x2-(x2/ec1).^2*fcm)/(1+(Eci*ec1/fcm-2)*(x2/ec1))
% Bereich 3 der Arbeitslinie
x3=0.0022:0.0001:0.04
y3=((2+gc*fcm*ec1)/(2*fcm)-gc*x3+gc/(2*ec1)*x3.^2).^-1
plot(x1,y1,'b')
hold on
plot(x,y,'r')
hold on
plot(x3,y3,'g')