### #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;
delT_isoComp = 0.03;
delT_isoExp = 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, 'I', delT_isoExp , 1/IsoComp);    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));
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 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);
if multiply >1;
display('adiabatic compression ...'); sign = [0 1 0];
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,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.

1. 