#78 Brownian Carnot engine-II 微觀卡諾引擎-II

Hello, everyone. Today we'll continue our discussion about Brownian Carnot engine, by directly build a virtual one. To understand the code, you need to know how to do a computer simulation of a Langevin equation. Please refer to our previous post to check it out.
/*------------Divider--------------*/
So let's see our code written in MATLAB...

function  Result = Brownian_Carnot_Module
    %% permanent setting
    tic;
    global dt Dt Dt_obs kB R eta zeta T0 A perN_hit perN_obs N_multiple m N_record_init
    dt = 1E-7; Dt = 1E-5; Dt_obs = 1E-4;
    kB = 1.38E-23; R = 500E-9; eta = 8.90E-4;
    zeta = 6*pi*eta*R; rho = 1040; m = 4/3*pi*R^3*rho;
    T0 = 300; A = sqrt(2*kB*T0*zeta*dt);
    perN_hit = round(Dt/dt); perN_obs = round(Dt_obs/dt); N_multiple = 300;
    %initiation
    delT_init = 0.02;
    T_set_init = 300;
    kappa_init = 1E-6;
    IsoComp = 4;
    AdiaComp = 4;
    delT_isoComp = 0.03;
    delT_adiaComp = 0.03;
    delT_isoExp = 0.03;
    delT_adiaExp = 0.03;
    State.x_record = []; State.T_record= []; State.T_set = []; State.kappa_record = [];
    State.sign = []; State.time = 0; State.ProcessRecord = {}; State.QRecord = [];
    %% Carnot Engine
    State = Carnot_Initiation(State, delT_init, T_set_init , kappa_init);    toc;
    State = Carnot_Process(State, 'I', delT_isoComp , IsoComp);    toc;
    State = Carnot_Process(State, 'A', delT_adiaComp , AdiaComp);    toc;
    State = Carnot_Process(State, 'I', delT_isoExp , 1/IsoComp);    toc;
    State = Carnot_Process(State, 'A', delT_adiaExp , 1/AdiaComp);    toc;
    %% time course figure summary
    x_record = State.x_record;    kappa_summary = State.kappa_record;
    v_record = abs((x_record(:, 2:end) - x_record(:, 1:(end-1)))/Dt_obs).*repmat(sqrt(Lt(kappa_summary(1:(end-1)))), N_multiple ,1);
    Sx = EntropyCal(x_record);
    Sv = EntropyCal(v_record);
    S = Sx(1:(end-1)) + Sv; %calculate entropy
    Time = Dt_obs:Dt_obs:(State.time - Dt_obs);
    Temp = State.T_record; T_set = State.T_set;
    Result.x_record = x_record; Result.v_record = v_record; Result.kappa_summary = kappa_summary; Sign = State.sign;
    Result.Sx = Sx; Result.Sv = Sv; Result.S = S; Result.Time = Time; Result.Temp = Temp; Result.ProcessRecord = State.ProcessRecord;
    figure;
    [ax,h1,h2] = plotyy(Time, S, Time, Temp(1:(end-1))); hold on;
    set(h1,'linewidth',3); set(h2,'linewidth',2); set(ax, 'FontSize', 14);
    xlabel('Time','FontSize',16);
    ylabel('Entropy','FontSize',16);
    yLim = get(gca, 'YLim');
    [nProcess,~]=size(State.ProcessRecord);
    for i = 1:1:nProcess
        line([State.ProcessRecord{i,2}, State.ProcessRecord{i,2}], yLim,'LineWidth',2);
        text(State.ProcessRecord{i,2},yLim(2),State.ProcessRecord{i,1},'FontSize',20);
    end
    toc;
    %% Other Plots
    %T-kappa
    figure; scatter(kappa_summary, Temp,[],Sign); hold on;
    scatter(kappa_summary, T_set ,[] ,Sign,'filled');
    xlabel('\kappa (N/m)','FontSize',16);
    ylabel('Temperature(K)','FontSize',16);
    set(gca, 'FontSize', 14);
    %F-kappa
    figure; F_kappa_exp = mean(x_record.^2, 1)/2;
    F_kappa_theo = 1/2*kB*T_set./kappa_summary;
    scatter(kappa_summary, F_kappa_exp,[],Sign); hold on;
    scatter(kappa_summary, F_kappa_theo,[],Sign,'filled');
    ylabel('F_{\kappa} (m^2)','FontSize',16);
    xlabel('\kappa (N/m)','FontSize',16);    set(gca, 'FontSize', 14);
    %T-S
    figure; S_baseline = mean(S(1:N_record_init));
    S_theo = S_baseline + log(T_set/T_set_init)-1/2*log(kappa_summary/kappa_init);
    scatter(S, Temp(1:end-1),[],Sign(1:end-1,:)); hold on; scatter(S_theo, T_set ,[] ,Sign,'filled');
    xlabel('S/k_B (1/K)','FontSize',16);
    ylabel('Temperature(K)','FontSize',16);
    set(gca, 'FontSize', 14);
    S_cal =  0.5*log(2*pi*kB*T0/kappa_init)+1/(2*pi)+0.5*log(2*pi*kB*T0/m)+1/(2*pi);
    display([S_baseline, S_cal]);
    %efficiency
    dW = F_kappa_exp(1:end-1).*(kappa_summary(2:end)-kappa_summary(1:end-1));
    dQ = kB*(Temp(2:end)-Temp(1:end-1))- dW;
    recordOrNot = logical(State.QRecord(1:end-1));
    efficiency = - sum(dW)/sum(dQ(recordOrNot));
    efficiency_theo = 1- 1/sqrt(AdiaComp);
    display([efficiency, efficiency_theo]);
    Result.efficiency = efficiency; Result.efficiency_theo = efficiency_theo;
    Result.dQ =dQ; Result.dW =dW; Result.recordOrNot = recordOrNot;
    %% end
    beep;
end

function correction = Lt(kappa)
    global Dt_obs zeta m
    f = 1/Dt_obs; fp = zeta/m; fk = kappa/(2*pi*zeta); f0 = sqrt(fp*fk); f1 = sqrt(fp^2/4 - f0.^2);
    correction = 1./(2.*f.^2).*(1./f0.^2 + 1./f1.*(exp(-f1/f - fp/(2*f))./(fp+2*f1) - exp(f1/f - fp/(2*f))./(fp-2*f1))).^(-1);
end

function S = EntropyCal(record)
    [~, nLength] = size(record);
    S = zeros(1, nLength);
    for i = 1:1:nLength
        h = histogram(record(:,i), 25,'Normalization','pdf');
        pdf = h.Values; d = h.BinWidth; close all;
        S(i) = - nansum(pdf.*log(pdf))*d;
    end
end

function State = Carnot_Initiation(State, delT, TempSet, kappa_init)
    %%
    global dt perN_obs T0 kB N_multiple A zeta perN_hit Dt N_record_init
    display('initiation...');
    State.ProcessRecord{1,1}='initiation';State.ProcessRecord{1,2}=State.time;
    N_init = delT/dt;    N_record_init = round(N_init/perN_obs);
    if TempSet < T0
        display('Error!');display('The setted temperature should always be higher than 300K.')
        return
    end
    T_set_init = TempSet - T0;
    A1_init = sqrt(2*kB*T_set_init*zeta*Dt);
    %%
    x_record_init = zeros(N_multiple, N_record_init);
    for j = 1:1:N_multiple
        x_temp = zeros(1, N_init);
        x_temp(1) = sqrt(kB*TempSet/kappa_init)*randn(1,1);
        for i = 2:1:N_init
            x_old = x_temp(i-1);
            x_temp(i) = x_old + 1/zeta*(-dt*(kappa_init*x_old)+A*randn(1,1));
            if mod(i, perN_hit) == 0
                x_temp(i) = x_temp(i) + 1/zeta*(A1_init*randn(1,1));
            end
            if mod(i, perN_obs) == 0
                x_record_init(j, i/perN_obs) = x_temp(i);
            end
        end
    end
    %%
    T_eff_init = kappa_init*mean(x_record_init.^2,1)/kB;
    State.x_record = [State.x_record x_record_init];
    State.T_record = [State.T_record T_eff_init];
    State.T_set = [State.T_set TempSet*ones(1, N_record_init)];
    State.kappa_record = [State.kappa_record kappa_init*ones(1, N_record_init)];
    State.sign = vertcat(State.sign, repmat([0 0 0],N_record_init,1));
    State.time = State.time + delT; State.QRecord = [State.QRecord zeros(1, N_record_init)];
End

function State = Carnot_Process(State, IsoAdia, delT, multiply)
    %%
    global dt perN_obs T0 kB N_multiple A zeta perN_hit Dt
    N = delT/dt;    N_record_Process = round(N/perN_obs);
    kappa_init = State.kappa_record(end);    kappa_endProcess = multiply*kappa_init;
    kappa_Process = kappa_init:((kappa_endProcess-kappa_init)/(N-1)):kappa_endProcess;
    T_init = State.T_set(end);
    %%
    if strcmp(IsoAdia, 'I')
        if multiply >1; 
            display('isothermal compression ...'); sign = [0 0 1];
            State.ProcessRecord{end+1,1}='isothermal compression';
            State.ProcessRecord{end,2}=State.time; 
            State.QRecord = [State.QRecord zeros(1, N_record_Process)];
        else 
            display('isothermal expansion ...'); sign = [1 0 0];
            State.ProcessRecord{end+1,1}='isothermal expansion';
            State.ProcessRecord{end,2}=State.time;
            State.QRecord = [State.QRecord ones(1, N_record_Process)];
        end
        T_set_Process = (T_init - T0)*ones(1, N); A1_Process = sqrt(2*kB*T_set_Process*zeta*Dt);
    elseif strcmp(IsoAdia, 'A')
        if multiply >1; 
            display('adiabatic compression ...'); sign = [0 1 0];
            State.ProcessRecord{end+1,1}='adiabatic compression';
            State.ProcessRecord{end,2}=State.time;
            State.QRecord = [State.QRecord ones(1, N_record_Process)];
        else 
            display('adiabatic expansion ...'); sign = [1 0 1];
            State.ProcessRecord{end+1,1}='adiabatic expansion';
            State.ProcessRecord{end,2}=State.time;
            State.QRecord = [State.QRecord ones(1, N_record_Process)];
        end
        Const = T_init^2/kappa_init;
        T_set_Process = round(sqrt(Const*kappa_Process) - T0, 5); 
        A1_Process = sqrt(2*kB*T_set_Process*zeta*Dt);
    else display('incomprehensible input');
    end
    x_record_Process = zeros(N_multiple, N_record_Process);
    for j = 1:1:N_multiple
        x_temp = zeros(1, N);
        x_old = State.x_record(j, end);
        x_temp(1) = x_old + 1/zeta*(-dt*(kappa_Process(1)*x_old)+A*randn(1,1));
        for i = 2:1:N
            x_old = x_temp(i-1);
            x_temp(i) = x_old + 1/zeta*(-dt*(kappa_Process(i)*x_old)+A*randn(1,1));
            if mod(i, perN_hit) == 0
                x_temp(i) = x_temp(i) + 1/zeta*(A1_Process(i)*randn(1,1));
            end
            if mod(i, perN_obs) == 0
                x_record_Process(j, i/perN_obs) = x_temp(i);
            end
        end
    end
    kappa_Process_abre = kappa_Process(perN_obs: perN_obs: N);
    T_eff_Process = kappa_Process_abre.*mean(x_record_Process.^2,1)/kB;
    T_set_Process_abre = T_set_Process(perN_obs: perN_obs: N);
    %%
    State.x_record = [State.x_record x_record_Process];
    State.T_record = [State.T_record T_eff_Process];
    State.kappa_record = [State.kappa_record kappa_Process_abre];
    State.T_set = [State.T_set T_set_Process_abre+T0];
    State.sign = vertcat(State.sign, repmat(sign,N_record_Process,1));
    State.time = State.time + delT;
end

It looks a little complicated. Our simulation is composed of an initial waiting period and a full Carnot cycle. This is run for 300 times to take average and calculate entropy. The entropy is calculated from the probability density function of the location and speed of the microsphere using Shannon formula (that is, ). After these simulation, we made some plot to display our results.
In fact, it would be great if our readers try to modify these codes. For example, by adjusting the rate of compression and expansion, or change the order of each process.
Before we delve to our results, I have to briefly raise 3 issues. First, in the above code, we estimate the baseline entropy by taking the average of the entropy measured in initial waiting period. This is in fact not necessary. The entropy of the system could be calculated explicitly as

This analytical solution reads -20.3568, which is basically the same as our average in initial waiting period, -20.2964. This is left for the readers to prove that.

Secondly, the evaluation of efficiency is a little bit different from classical Carnot engine. In classical Carnot engine, the efficiency is calculated as:

However, in its stochastic counterpart, the heat flow during the namely "adiabatic" process must be taken into consideration. Thus the efficiency of a Brownian Carnot engine should be calculated as:

Thirdly, since we observe the microsphere and estimate its speed by a rather long time interval, it is necessary to correct our average speed to approximate the instantaneous speed. The authors have proved the validity of their correction formula.

And finally, let's see our results:

This is the summary of our experiment time course. As expected, the entropy decreases during isothermal compression, and nearly holds the same during 2 adiabatic processes.

The filled color lines indicate what a classical Carnot engine would behave in our experimental setting. There are some fluctuations but our stochastic engine acts like his classical counterpart. However, this would not be the case if we hasten the rate of compression and expansion. What do you expect to change? What is the source of irreversibility in Brownian Carnot engine?

This is the F-κ plot, which you could interpret as a P-V plot. The area enclosed in the cycle is the work our engine done.

And this is a T-S plot. In our simulation here, the efficiency is 0.4427, which is close to the efficiency of its classical counterpart, 0.50, because we compress and expand the system rather slowly and take average from 300 cycles. Efficiency above classical limit is possible, however, due to the stochastic nature of our system, and it is in fact a big issue to investigate the distribution of efficiency of our system. But we will stop here for it takes time to explore.