[求助] 塑性变形对铁磁材料磁化过程的影响物理模型

function [x,y]=lx15()ms=1.585e6;  %饱和磁化强度r1=7e-18;r11=-1e-25;r2=-3.3e-30;r22=2.1e-38;u0=4*pi*(1e-7);G1=0.81e-6;G2=3.8e-12;G3=G1;G4=G2;alphak=0.76;C11=2.863e11;C12=1.41e11;poisson=C12/(C11+C12);yield=C11-2*C12*poisson;G=yield/(2*(1+poisson));b=2.5e-10;d=20e-6;k0_u0=500;a0=600;zetad0=1e12;k0=k0_u0*4*pi*10^(-7);f=0.075;sigmay=235e6;epsilonp=0.0001;sigma=yield*epsilonp+sigmay;zetad=(f*yield*epsilonp/(alphak*G*b)+zetad0^0.5)^2;k=(G1+G2/d)*zetad^(1/2)*k0;a=(G3+G4/d)*zetad^(1/2)*a0;alpha=8.44e-6;%求解x12=linspace(0,7e3,100);   %H=0:7000A/mx13=linspace(7e3,-7e3,200);  %H=7000:-7000A/m;x14=linspace(-7e3,7e3,200);  %H=-7000:7000A/m;[x,y11]=ode45(@ode_fun11,x12,[0.0001,0]); %anhysteretic loop when H=0:7000A/m[x,y12]=ode45(@ode_fun12,x13,[y11(100,1),y11(100,2)]);  %hysteretic process when H=7000:-7000A/m[x,y13]=ode45(@ode_fun11,x14,[y12(200,1),y12(200,2)]);  %hysteretic process when H=-7000:7000A/mMan_coth1=ms*(coth((x13'+alpha*y12(:,1)+3*sigma/u0*((r1+r11*sigma)*y12(:,1)+2*(r2+r22*sigma)*y12(:,1).^3))/a)-a./(x13'+alpha*y12(:,1)+3*sigma/u0*((r1+r11*sigma)*y12(:,1)+2*(r2+r22*sigma)*y12(:,1).^3))); %Man公式计算结果Man_coth2=ms*(coth((x14'+alpha*y13(:,1)+3*sigma/u0*((r1+r11*sigma)*y13(:,1)+2*(r2+r22*sigma)*y13(:,1).^3))/a)-a./(x14'+alpha*y13(:,1)+3*sigma/u0*((r1+r11*sigma)*y13(:,1)+2*(r2+r22*sigma)*y13(:,1).^3)));  %Man公式计算结果figureplot(x13',Man_coth1,x14',Man_coth2)title('Man公式计算结果');figureplot(x13',y12(:,2),'r',x14',y13(:,2),'b');grid on;xlabel('H(A/m)');ylabel('M(A/m)');title('toatal magnetization M');figureplot(x13',y12(:,1),'r',x14',y13(:,1),'b');grid on;xlabel('H(A/m)');ylabel('Man(A/m)');title('anhystersis magnetization Man');He12=x13'+alpha*y12(:,2)+3*sigma/u0*((r1+r11*sigma)*y12(:,2)+2*(r2+r22*sigma)*y12(:,2).^2);He13=x14'+alpha*y13(:,2)+3*sigma/u0*((r1+r11*sigma)*y13(:,2)+2*(r2+r22*sigma)*y13(:,2).^2);figureplot(x13',He12,'r',x14',He13,'b');grid on;xlabel('H(A/m)');ylabel('He(A/m)');title('effective magnetic field He');figureBe1=u0*(x13'+y12(:,2)+3*sigma/u0*((r1+r11*sigma)*y12(:,2)+2*(r2+r22*sigma)*y12(:,2).^2));Be2=u0*(x14'+y13(:,2)+3*sigma/u0*((r1+r11*sigma)*y