Диссертация (1103090), страница 23
Текст из файла (страница 23)
Случайная сила ξ задается температурой в соответствии с флуктуационно-диссипационной теоремой:⟨ξ(t1 )ξ(t2 )⟩ = 2kB T m.γШаг интегрирования выбирался достаточно малым (∼ 0.01τ ).112500κoH = 0.7450∆, σLJ4001.43503002.92502002.42.52.62.7t, τ2.82.93.0×105Рисунок А.3.
Изменение полуширины профиля концентрации малых ионов, ∆, со временемв ходе измерения силы взаимодействия мембран. Компьютерное моделирование проводилось для отношения зарядов ионов Z̃ = −1, а безразмерная ширина щели варьировалась отκo H = 0.7 до κo H = 2.9.Сначала помещали ионы в пространство между мембранами, после чего система моделировалась некоторое время для достижения равновесного состояния.
Релаксацию системык равновесному состоянию контролировали с помощью временной зависимости полуширины профиля концентрации ионов ∆. Как видно из рис. А.1, кинетическая энергия системыпрактически не изменяется со временем в процессе уравновешивания и находится вблизизначения 3/2kB T , поэтому она не может служить характерным параметром для контролярелаксации системы к равновесию. На рис. А.2 мы показываем, что можно использовать полуширину (или среднеквадратичное отклонение) профиля концентрации малых ионов, ∆, вкачестве такого параметра. Действительно, ∆ медленно увеличивается со временем по меревыхода малых ионов из щели во внешний раствор. Измерения свойств системы, профилейконцентрации и потенциала, силы взаимодействия мембран производились в равновесномсостоянии.
При этом профиль концентрации и ∆ не изменялись со временем, см. рис. А.3.Далее приведем код для компьютерного моделирования для пакета ESPResSo на языкеTCL.А.5 Численное решение уравнения ПБ для мембранКод в вычислительной системе Mathematica:ZP=.;√ ([ ][ ])Xout[ϕ_,ϕs_]:= h2 + 2 Exp ϕ2 − Exp ϕs2[]1Xin[ϕ_,ϕm_]:=NIntegrate √,{x,ϕm,ϕ}112((Exp[−x]− ZPExp[−ZPx]−Exp[−ϕm]+ ZPExp[−ZPϕm]))“Xin[ϕs,ϕm]:=h/2”;ϕsurf[ϕm_]:= −1Log[−ZPExp[−ϕm] + Exp[−ZPϕm]]ZPZP = −1; h = 2;113Plot[Xin[ϕsurf[ϕm],ϕm],{ϕm,0.01,10}]ClearAll[φ]sol = FindRoot[{Xin[ϕsurf[y],y] == h/2},{y,0.1,0.2},AccuracyGoal → 4,PrecisionGoal → 4]NIntegrate::nlim : x = y is not a valid limit of integration. ⟩⟩{y → 0.337706@ − 9.187876825590759`*∧ -6i}2Xin[ϕsurf[Abs[y]],Abs[y]]/.%NIntegrate::nlim : x = Abs[y] is not a valid limit of integration.
⟩⟩2.00013inv[f_,s_]:=Function[{t},s/.FindRoot[f − t,{s,1}]]ϕm:=Abs[y/.sol]ϕs:=ϕsurf[ϕm]ClearAll[X]X = {}; npoints = 500;δϕ = ϕs − ϕm;ϕs;δϕ2[{= 50iDo ϕ = ϕm + npoints∗ δϕ; X = Append[ {[[[]]] ()}]}]iiX, Abs Evaluate Xin ϕm + npoints ∗ δϕ,ϕm, ϕm + npoints ∗ δϕ,{i,0,npoints} ;[{iDo ϕ = ϕs + npoints∗ δϕ2; X = Append[X,{Abs[Evaluate[Xout[ϕs + i ∗ δϕ2,ϕs]]],(ϕs + i ∗ δϕ2)}]},{i,0,npoints}];string = “path to save”;Export[string,X,“table”];Shooting method(*φ”[x]== − Exp[−φ[x]]*)eq1 = φ”[x] == −Exp[−φ[x]] + Exp[−Zφ[x]]; (*The inner PB equation to be solved*)Z = −15; h = 1;φh[s_?NumberQ]:=First[φ[h/2]/.NDSolve[{eq1,φ[0] == ϕm,φ′ [0] == s},φ,{x,0,h/2}]](*dependence of solution upon the derivative ofpotential atthe initial point. So we have initial valueproblem.
Wemust choose the value of initial slope so thatthe surface pot. become equal ϕs*)Plot[φh[s],{s, − 4,0.1}]114ϕin[x_] = First[φ[x]/.NDSolve[{eq1,φ[0] == ϕm,φ′ [0] == s/.FindRoot[φh[s] == ϕs,{s,0}]},φ,{x,0,h/2}]]InterpolatingFunction[{{0.,0.5}},<>][x] ][()1ϕout[x_] = 2Log Exp[ϕs/2] + √2 x − h2 ;ϕtotal[x_]:=Piecewise[{{ϕin[x],0 ≤ x < h/2},{ϕout[x],x ≥ h/2}}]Plot[ϕtotal[x],{x,0,10},PlotRange → {{0,2},{0,2}},PlotPoints → 1000]Abs[FindRoot[φh[s] == ϕs,{s, − 10}]]{s → 0.650309}φh[0.6503]0.185376Programforsurfaceandmiddlepotentialcalculation;Phim = {};midpot[i_]:=Abs[y/.FindRoot[{Xin[ϕsurf[y],y] == i/2},{y,0,0.00001},AccuracyGoal → 7,PrecisionGoal → 7]]For[i = 0,i ≤ 500,i++,Phim = Append[{[]}]]Phim, 10(i−400)/100 //N,midpot 10(i−400)/100;string2 = “path to save”;Export[string2,Phis,“table”];Phim[[All,2]]//MatrixFormPhis = Transpose[{Phim[[All,1]],ϕsurf[Phim[[All,2]]]}]//MatrixForm115Приложение БЧисленное решение НТПБ для заряженных мембран врастворе электролитаЧисленное решение ур-ий (2.27) и (2.28) производилось методом коллокаций [160],реализованном с помощью скриптового языка python и библиотеки “scikits.bvp1lg”.
Для решения использовали безразмерные переменные для координаты: x → κi x. Вместо решенияуравнений для внешней и внутренней области по отдельности запишем уравнение НТПБследующим образом:()∂ 2 ϕ(x)−ϕo−Z̃ϕo=−e−U(x)e− qδ(x − H/2)/d,∂x2(Б.1)где U = 0 для |x| < H/2 и 1 в остальных случаях. Вместо граничного условия на поверхности мембраны мы добавляем дополнительный член в правую часть уравнения. Данный членучитывает тот факт, что мембраны заряжены посредством равномерного распределения поверхностного заряда σ в тонком слое толщиной d.
Толщина мембраны выбиралась так, чтоd ≪ H.Перед решением данное ур. (Б.1) для ϕ преобразовывали в систему уравнений: y1 = ϕи y2 = dϕ/dx. Система решаемых уравнений теперь выглядит следующим образом:y1′= y2 ;y2′=−κ2i(−y1e− U (x)e−Z̃y1).(Б.2)Данное уравнение решалось на промежутке от 0 до H/2 + 16λD . Устанавливали граничноеусловие симметрии системы как y2 = 0 при x = 0 и y2 = 0 при x = H/2 + 16. Расчетыостанавливались, когда решение находилось с точностью до 10−10 .
В качестве выходныхпараметров выводились потенциал y1 , потенциал в центре щели y1 (0) и мембранный потенциал y1 (H/2).Код для расчета потенциала приведен ниже:# загрузка~основных~библиотекimport numpy as npimport scikits.bvp1lg as bvp5 import math,os#~~ширина~мембраныd=0.01# расстояние~между~узлами~сетки10 dx = 0.001116# отношение~зарядов~ионовZ= -1.0phis0= -1/Z*np.log(1-Z)15 ###### Surface~properties# dimensionless~surface~charge# безразмерный~заряд~мембраныq =-0.3520 # Ступенчатая~функцияdef delta(x,d,x0):return 0.5 * (np.sign(x-x0) + 1) - 0.5* (np.sign(-d+x-x0) +1.0)# Потенциал~мембраны,~полученный~ с ~использованием~ЛТПБ25 def phis_DH(q,H):return ((q)*H*np.cosh(H/2.)+H*np.sinh(H/2.))/(H*np.sinh(H/2.) + np.sqrt(2.)*H*np.cosh(H/2.))# Потенциал,~полученный~с~использованием~ЛТПБdef pot_DH(x,q,H):30return np.piecewise(x,[np.abs(x)<H/2.,np.abs(x)>=H/2.],[lambda x: (phis_DH(q,H)-1.)/np.cosh(H/2.)*np.cosh(x)+1.,lambda x: phis_DH(q,H)*np.exp(np.sqrt(2.)*(H/2.-np.abs(x)))])# Производнаяпотенциала , полученная~с~использованием~ЛТПБdef dpot_DH(x,q,H):return np.piecewise(x,[np.abs(x)<H/2.,np.abs(x)>=H/2.],[lambda x: (phis_DH(q,H)-1.)/np.cosh(H/2.)*np.sinh(x),lambda x: -np.sqrt(2.)*phis_DH(q,H)*np.exp(np.sqrt(2.)*(H/2.-np.abs(x)))])35def potential(H):# импорт~глобальных~переменныхglobal dxglobal d,q40def fsub(x,Y):"""The equations, in the form: Y' = f(x, Y)."""1174550y, dy = Yd_y= dyd_dy = ( -np.exp(-y)+ delta(x,L-H/2.,H/2.)*np.exp(y) )-q*delta(x,d,H/2.)/dreturn np.array([ d_y, d_dy ])def gsub(Y):Граничные"""~условия."""y, dy = Yreturn np.array([dy[0] - 0,dy[1] - 0])# Начальное~приближениеdef guess(x):55y = pot_DH(x,q,H)dy = dpot_DH(x,q,H)606570Y = np.array([y, dy ] )dm = fsub(x, Y)return Y, dmL = H/2.+16.0# channel~widthnx = (L)/dxnx = 32000.# количество_~точек# 2 times zero as we have two bcs at x=0 :boundary_points = [0, L]tol = 1e-10 * np.ones_like(boundary_points)# print toldegrees = [1,1]solution = bvp.colnew.solve(boundary_points,degrees, fsub,gsub,initial_guess=guess,tolerances=tol, vectorized=True,maximum_mesh_size=256000)x = np.linspace(dx, L, nx)y, dy = solution(x).transpose()75# вывод_~решенияreturn x,y, y[0], y[(np.abs(H/2.
- x)).argmin()]118Приложение ВВзаимодействие мембраны и заряженной поверхности:вывод формулы для расклинивающего давленияКак обсуждалось ранее, в равновесии электростатическое расклинивающее давлениесостоит из двух частей, а именно, давления, вызванного электростатической объемной силой (ρE), и осмотического давления [74]. Однако в ряде работ показано, что соотношениядля расчета давления, выведенные в рамках НТПБ, нуждаются в корректировке и не могутбыть непосредственно применены для расчетов в рамках ЛТПБ [67; 70]. Мы выведем соотношение для давления, руководствуясь тем, что теория должна быть самосогласованной.Условие механического равновесия предполагает, что электростатический потенциали объемный заряд удовлетворяют уравнению гидростастики0 = −∇p + ρE = ∇ · (T − I p) ≡ −∇ · Π,(В.1)где T – тензор электростатических напряжений Максвелла:[]εE2Tij =Ei Ej − δij.4π2(В.2)Тензор Π(x,y) = T(x,y) − I p(x,y) есть электростатическая составляющая расклинивающего давления.
Уравнение (В.1) показывает [9; 161], что для одномерной системыс однородными поверхностями в равновесии расклинивающее давление Π инвариантно ипринимает одно и то же значение в каждой точке x [137], но наиболее удобным выбором является точка симметрии, например, центр щели между двумя поверхностями. В этой точкеэлектростатические напряжения равны нулю, T = 0, а расклинивающее давление сводитсяк осмотическому [4; 5; 9; 78].Неоднородность поверхности приводит к тому, что давление Π в общем случае не инвариантно, что следует из ур. (В.1); и к тому, что симметрия нарушается и точку симметриивыше необходимо заменить другой подходящей.Мы предлагаем следующий подход к решению проблемы:1.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
