粘弾性解析の基礎①：弾性と粘性

%% 第1回：弾性と粘性 ― 材料点の最小数値実験
% MATLAB R2020a以降。追加Toolbox不要。
% このファイルを実行すると、figures/ にPNG・PDF、data/ にCSVを保存します。
% 応力: MPa、時間: s、ひずみ: 無次元、粘性係数: MPa s。
clear; close all; clc;
root = fileparts(mfilename('fullpath'));
figdir = fullfile(root,'figures'); datadir = fullfile(root,'data');
if ~exist(figdir,'dir'), mkdir(figdir); end
if ~exist(datadir,'dir'), mkdir(datadir); end
set(groot,'defaultAxesFontName','Yu Gothic','defaultTextFontName','Yu Gothic');
set(groot,'defaultAxesFontSize',11,'defaultLineLineWidth',1.8);

E = 1000; eta = 1000; eps0 = 0.01;
ramps = [1 10 100]; dt = 0.02; t = (0:dt:120)';
colors = [0.10 0.35 0.65; 0.85 0.38 0.10; 0.15 0.55 0.38];
labels = {'1 sで負荷','10 sで負荷','100 sで負荷'};
strain = zeros(numel(t),3); rate = strain; spring = strain; dashpot = strain;
workS = zeros(1,3); workD = workS;
for j = 1:numel(ramps)
    tr = ramps(j);
    strain(:,j) = eps0*min(t/tr,1);
    % 各区間の平均ひずみ速度。0番目の点は負荷開始直前の状態。
    rate(:,j) = [0; diff(strain(:,j))/dt];
    spring(:,j) = E*strain(:,j);
    dashpot(:,j) = eta*rate(:,j);
    de = diff(strain(:,j));
    % ばねは区間両端応力の平均、ダッシュポットは区間一定応力で仕事を計算。
    workS(j) = sum(0.5*(spring(1:end-1,j)+spring(2:end,j)).*de);
    workD(j) = sum(dashpot(2:end,j).*de);
    assert(max(abs(spring(:,j)-E*strain(:,j))) < 1e-10);
    assert(abs(workS(j)-0.5*E*eps0^2) < 1e-10);
    assert(abs(workD(j)-eta*eps0^2/tr) < 1e-10);
    inside = t > 0 & t < tr;
    assert(max(abs(dashpot(inside,j)-eta*eps0/tr)) < 1e-9);
    assert(max(abs(dashpot(t > tr+dt,j))) < 1e-9);
end

f=figure('Color','w','Position',[100 100 1000 920]);
tl=tiledlayout(3,1,'TileSpacing','compact','Padding','compact');
nexttile; hold on;
for j=1:3, plot(t,100*strain(:,j),'Color',colors(j,:)); end
ylabel('ひずみ [%]'); title('(a) 同じ最終ひずみを、異なる速さで与える');
legend(labels,'Location','southeast'); grid on; xlim([0 120]); ylim([0 1.15]);
nexttile; hold on;
for j=1:3, plot(t,spring(:,j),'Color',colors(j,:)); end
ylabel('ばねの応力 [MPa]'); title('(b) ばね：同じひずみなら同じ応力');
grid on; xlim([0 120]); ylim([0 11.5]);
nexttile; hold on;
for j=1:3, stairs(t,dashpot(:,j),'Color',colors(j,:)); end
xlabel('時間 [s]'); ylabel('ダッシュポットの応力 [MPa]');
title('(c) ダッシュポット：速く変形させるほど大きな応力');
grid on; xlim([0 120]); ylim([0 11.5]);
export_pair(f,figdir,'ramp_response');

f=figure('Color','w','Position',[100 100 1000 440]);
tiledlayout(1,2,'TileSpacing','compact','Padding','compact');
nexttile; hold on;
for j=1:3, plot(t,100*strain(:,j),'Color',colors(j,:)); end
title('(a) 負荷開始直後のひずみ'); xlabel('時間 [s]'); ylabel('ひずみ [%]');
xlim([0 2]); ylim([0 1.15]); grid on;
nexttile; hold on;
for j=1:3, stairs(t,dashpot(:,j),'Color',colors(j,:)); end
title('(b) 負荷開始直後の粘性応力'); xlabel('時間 [s]'); ylabel('応力 [MPa]');
xlim([0 2]); ylim([0 11.5]); grid on; legend(labels,'Location','northeast');
export_pair(f,figdir,'ramp_zoom');

% 一定応力を時刻0の直後から10秒間与え、その後に除荷する。
% ばねの瞬間変形は繰返し時刻で描く。慣性を無視した理想化。
sig0=1; thold=10;
tc=[0;0;(dt:dt:thold)';thold;(thold+dt:dt:20)'];
sc=[0;sig0;sig0*ones(round(thold/dt),1);0;zeros(round(thold/dt),1)];
ecS=sc/E; ecD=sig0/eta*min(tc,thold);
assert(abs(ecD(end)-0.01)<1e-12);
assert(ecS(end)==0);
f=figure('Color','w','Position',[100 100 1000 680]);
tiledlayout(2,1,'TileSpacing','compact','Padding','compact');
nexttile; plot(tc,sc,'Color',colors(1,:)); ylabel('応力 [MPa]');
title('(a) 1 MPaを10秒間与え、その後に除荷'); grid on; ylim([-0.1 1.2]);
nexttile; plot(tc,100*ecS,'Color',colors(1,:)); hold on;
plot(tc,100*ecD,'Color',colors(2,:)); xlabel('時間 [s]'); ylabel('ひずみ [%]');
title('(b) ばねは直ちに回復し、単独のダッシュポットは変形が残る');
legend('ばね','ダッシュポット','Location','southeast'); grid on; ylim([-0.05 1.15]);
export_pair(f,figdir,'creep_recovery');

writetable(table(ramps',eta*eps0./ramps',workS',workD', ...
    'VariableNames',{'ramp_time_s','viscous_stress_MPa','spring_work_MJ_m3','dashpot_work_MJ_m3'}), ...
    fullfile(datadir,'summary.csv'));
writetable(table(t,strain(:,1),strain(:,2),strain(:,3),spring(:,1),spring(:,2),spring(:,3), ...
    dashpot(:,1),dashpot(:,2),dashpot(:,3),'VariableNames', ...
    {'time_s','strain_1s','strain_10s','strain_100s','spring_1s_MPa','spring_10s_MPa', ...
    'spring_100s_MPa','dashpot_1s_MPa','dashpot_10s_MPa','dashpot_100s_MPa'}),fullfile(datadir,'ramp_history.csv'));
writetable(table(tc,sc,ecS,ecD,'VariableNames', ...
    {'time_s','stress_MPa','spring_strain','dashpot_strain'}),fullfile(datadir,'creep_history.csv'));
disp('PASS: 構成式・解析解・仕事・除荷後ひずみを確認しました。');
disp(readtable(fullfile(datadir,'summary.csv')));

function export_pair(f,folder,name)
    exportgraphics(f,fullfile(folder,[name '.png']),'Resolution',180,'BackgroundColor','white');
    exportgraphics(f,fullfile(folder,[name '.pdf']),'ContentType','vector','BackgroundColor','white');
end

function viscoelasticity_videos
% 記事第1回の2本のMP4を再生成。追加Toolbox不要。
% 描画の伸びは誇張。動画の再生時間とモデル時間を区別して表示。
root=fileparts(mfilename('fullpath')); out=fullfile(root,'videos');
if ~exist(out,'dir'),mkdir(out);end
set(groot,'defaultAxesFontName','Yu Gothic','defaultTextFontName','Yu Gothic');
set(groot,'defaultAxesFontSize',12);
E=1000;eta=1000;fps=20;
for mode=1:2
    f=figure('Color','w','Position',[60 60 1280 720],'MenuBar','none','ToolBar','none','Visible','off');
    annotation(f,'textbox',[.02 .92 .96 .07],'String', ...
        pick(mode==1,'動画1：同じ最終ひずみでも、速さによって粘性応力が変わる', ...
        '動画2：伸ばす → 長さを保つ → 力をゼロにする'), ...
        'FontName','Yu Gothic','FontSize',22,'FontWeight','bold','EdgeColor','none','HorizontalAlignment','center');
    clock=annotation(f,'textbox',[.02 .86 .96 .05],'String','','FontSize',15,'EdgeColor','none','HorizontalAlignment','center');
    annotation(f,'textbox',[.02 .005 .96 .05],'String', ...
        '青：ばね   橙：ダッシュポット   破線：変形前の端   ／   伸びは拡大表示。慣性を無視した理想要素。', ...
        'FontSize',12,'EdgeColor','none','HorizontalAlignment','center');
    colors={[.08 .34 .62],[.73 .29 .07]};
    models=cell(2,2); traces=cell(2,2); axesTrace=cell(2,2);
    for col=1:2
        a=axes(f,'Position',[.045+(col-1)*.48 .43 .43 .40]);
        axis(a,[0 340 0 200]);axis(a,'off');hold(a,'on');
        if mode==1
            title(a,pick(col==1,'速い負荷：1秒で1%','遅い負荷：10秒で1%'),'FontSize',17);
        else
            title(a,pick(col==1,'ばね：力を抜くと戻る','ダッシュポット：伸びたまま'),'FontSize',17);
        end
        for row=1:2
            if mode==1, kind=row; y=145-(row-1)*105;
            else, kind=col;y=110;if row==2,continue;end,end
            models{row,col}=make_element(a,kind,y,colors{kind});
        end
        for row=1:2
            at=axes(f,'Position',[.06+(col-1)*.48 .29-(row-1)*.16 .40 .12]);hold(at,'on');grid(at,'on');
            horizon=pick(mode==1,12,7);
            xlim(at,[0 horizon]);ylim(at,pick(row==1,[-.05 1.15],[-.5 11.5]));
            ylabel(at,pick(row==1,'ひずみ [%]','応力 [MPa]'));
            if row==2,xlabel(at,'モデル時間 [s]');end
            traces{row,col}=plot(at,nan,nan,'LineWidth',2,'Color',pick(mode==1,colors{1},colors{col}));
            axesTrace{row,col}=at;
        end
        if mode==1
            traces{3,col}=plot(axesTrace{2,col},nan,nan,'LineWidth',2,'Color',colors{2});
        end
    end
    if mode==1
        name='video01_loading_speed';ts=linspace(0,12,321);ramps=[1 10];
    else
        name='video02_hold_release';ts=linspace(0,7,281);ramps=[];
    end
    v=VideoWriter(fullfile(out,[name '.mp4']),'MPEG-4');v.FrameRate=fps;v.Quality=92;open(v);
    saved=false(1,3);
    for k=1:numel(ts)
        t=ts(k);
        if mode==1
            clock.String=sprintf('モデル時間 %.2f s   ／   両条件とも最終ひずみ1%%。時間は同じ速さで進みます。',t);
            for col=1:2
                e=.01*min(t/ramps(col),1);rate=.01/ramps(col)*(t<ramps(col));
                update_element(models{1,col},e,E*e,rate,'ばね');
                update_element(models{2,col},e,eta*rate,rate,'ダッシュポット');
                et=.01*min(ts(1:k)/ramps(col),1);
                set(traces{1,col},'XData',ts(1:k),'YData',et*100);
                set(traces{2,col},'XData',ts(1:k),'YData',E*et);
                set(traces{3,col},'XData',ts(1:k),'YData',eta*.01/ramps(col)*(ts(1:k)<ramps(col)));
            end
        else
            stage=1+(t>=1)+(t>=4);
            names={'① 1秒で1%まで伸ばす','② 長さを固定：動かさずに保つ','③ 外力ゼロ：ばねの戻りは理想的に瞬間'};
            clock.String=sprintf('モデル時間 %.2f s   ／   %s',t,names{stage});
            eh=.01*min(ts(1:k),1);springe=eh.*(ts(1:k)<4);
            current=[springe(end),eh(end)];stresses=[E*current(1),eta*.01*(t<1)];
            for col=1:2
                update_element(models{1,col},current(col),stresses(col),.01*(t<1),pick(col==1,'ばね','ダッシュポット'));
                set(traces{1,col},'XData',ts(1:k),'YData',pick(col==1,springe,eh)*100);
                set(traces{2,col},'XData',ts(1:k),'YData',pick(col==1,E*springe,eta*.01*(ts(1:k)<1)));
            end
        end
        drawnow;
        frame=getframe(f);writeVideo(v,frame);
        for j=1:3
            mark=pick(mode==1,[.5 4 11],[.8 2.5 5]);
            if ~saved(j) && t>=mark(j)
                imwrite(frame.cdata,fullfile(out,sprintf('%s_check%02d.png',name,j)));saved(j)=true;
            end
        end
        if mode==1 && k==round(numel(ts)*.35) || mode==2 && k==round(numel(ts)*.37)
            imwrite(frame.cdata,fullfile(out,[name '_poster.png']));
        end
    end
    close(v);close(f);
    r=VideoReader(fullfile(out,[name '.mp4']));assert(r.Width==1280 && r.Height==720);
    assert(r.Duration>10);fprintf('PASS: %s, %.2f seconds, %dx%d\n',name,r.Duration,r.Width,r.Height);
end
end

function m=make_element(a,kind,y,color)
plot(a,[25 25],[y-18 y+18],'Color',[.5 .55 .6],'LineWidth',4);
plot(a,[180 180],[y-24 y+24],'--','Color',[.6 .65 .7]);
m.a=a;m.kind=kind;m.y=y;m.color=color;
m.line=plot(a,nan,nan,'Color',color,'LineWidth',2.7);
m.end=plot(a,nan,nan,'Color',color,'LineWidth',5);
m.force=quiver(a,0,0,0,0,0,'Color',color,'LineWidth',2,'MaxHeadSize',.8);
m.label=text(a,170,y-35,'','FontName','Yu Gothic','FontSize',15,'Color',color,'HorizontalAlignment','center');
m.motion=text(a,170,y+31,'','FontName','Yu Gothic','FontSize',12,'Color',[.1 .45 .35],'HorizontalAlignment','center');
end

function update_element(m,e,s,rate,label)
xend=180+e/.01*65;y=m.y;
if m.kind==1
    xx=[25 45 linspace(48,xend-18,17) xend];yy=[y y y+10*(-1).^(0:16) y];
else
    piston=xend-75;
    xx=[25 65 65 185 nan 65 65 185 nan piston piston nan piston xend];
    yy=[y y y+18 y+18 nan y y-18 y-18 nan y-15 y+15 nan y y];
end
set(m.line,'XData',xx,'YData',yy);set(m.end,'XData',[xend xend],'YData',[y-12 y+12]);
set(m.force,'XData',xend+7,'YData',y,'UData',35*(s>1e-9),'VData',0);
m.label.String=sprintf('%s：ひずみ %.2f%% ／ 応力 %.2f MPa',label,100*e,s);
m.motion.String=pick(rate>0,'右端を動かす →','右端は静止');
end

function y=pick(c,a,b)
if c,y=a;else,y=b;end
end