``` %% CNC自组装动力学仿真主程序（边界条件修正版） clear; clc; close all; %--------------------------- % 仿真参数设置 %--------------------------- Lx = 200; % x方向模拟区域尺寸(μm) Ly = 200; % y方向尺寸 Nx = 16; % x方向网格数 Ny = 128; % y方向网格数 dx = Lx/Nx; % 空间步长 dy = Ly/Ny; dt = 0.1; % 时间步长 total_time = 750; % 总模拟时间 steps = total_time/dt; % Landau-de Gennes参数 a = -0.5; % 温度相关线性系数 B = 1.0; % 三阶项系数 C = 1.0; % 四阶项系数 L1 = 1.0; % 弹性常数L1 L2 = 0.5; % 手性扭曲耦合系数 p0 = 15; % 手性螺距(μm) np = p0/(2*pi); % 螺距扭曲参数 %--------------------------- % 初始化场量 %--------------------------- Q = struct('xx', zeros(Ny,Nx), 'xy', zeros(Ny,Nx), 'yy', zeros(Ny,Nx)); % 初始条件：各向同性态添加随机扰动，并固定y边界 rng(2024); Q.xx = 0.01*(rand(Ny,Nx)-0.5); Q.xy = 0.01*(rand(Ny,Nx)-0.5); Q.yy = -Q.xx; % 固定y边界条件（取向沿x轴） Q.xx(1,:) = 0.01; % y=0边界 Q.xx(end,:) = 0.01; % y=Ly边界 Q.xy(1,:) = 0; Q.xy(end,:) = 0; Q.yy = -Q.xx; %--------------------------- % 预定义微分算子 %--------------------------- laplacian = @(M) (circshift(M,[0 1]) + circshift(M,[0 -1]) + ... circshift(M,[1 0]) + circshift(M,[-1 0]) -4*M)/(dx^2); gradientX = @(M) (circshift(M,[0 1]) - circshift(M,[0 -1]))/(2*dx); gradientY = @(M) (circshift(M,[1 0]) - circshift(M,[-1 0]))/(2*dy); %--------------------------- % 时间演化主循环 %--------------------------- for step = 1:steps %=== 自由能导数计算 === Qnorm = Q.xx.^2 + 2*Q.xy.^2 + Q.yy.^2; dF_dQ = a*(Q.xx+1i*Q.xy) - B*(Q.xx.^2-Q.xy.^2+1i*(Q.xx.*Q.xy+Q.xy.*Q.yy)) + C*Qnorm.*(Q.xx+1i*Q.xy); %=== 弹性项计算 === elastic_term = L1*complex(laplacian(Q.xx), laplacian(Q.xy)) - 1i*L2*np*(gradientX(Q.xy)-gradientY(Q.xx)); %=== 时间推进 === Q.xx = Q.xx + dt*(real(dF_dQ) + real(elastic_term)); Q.xy = Q.xy + dt*(imag(dF_dQ) + imag(elastic_term)); Q.yy = -Q.xx; %=== 边界条件处理 === Q = applyBoundaryConditions(Q, dx, dy); % 应用修正后的边界条件 %=== 可视化 === if mod(step,50)==0 plotCNCstructure(Q, dx, dy, step*dt); end end %% 边界条件处理函数 function Q = applyBoundaryConditions(Q, dx, dy) % x方向周期性边界 Q.xx = applyXPeriodic(Q.xx); Q.xy = applyXPeriodic(Q.xy); Q.yy = applyXPeriodic(Q.yy); % y方向固定边界（取向沿x轴） Q.xx(1,:) = 0.01; % y=0边界 Q.xx(end,:) = 0.01; % y=Ly边界 Q.xy(1,:) = 0; Q.xy(end,:) = 0; Q.yy = -Q.xx; end function M = applyXPeriodic(M) % x方向周期性边界处理 M(:,1) = M(:,end-1); M(:,end) = M(:,2); end %% 增强型可视化函数 function plotCNCstructure(Q, dx, dy, currentTime) [X, Y] = meshgrid(0:dx:(size(Q.xx,2)-1)*dx, 0:dy:(size(Q.xx,1)-1)*dy); % 计算序参数和取向角 S = sqrt(Q.xx.^2 + 2*Q.xy.^2 + Q.yy.^2); theta = 0.5 * atan2(Q.xy, 0.5*(Q.xx - Q.yy)); % 角度范围[-π/2, π/2] % 创建HSV颜色映射 hue = (theta + pi/2)/pi; % 映射到[0,1] value = S/max(S(:)); hsvImg = cat(3, hue, ones(size(S)), value); rgbImg = hsv2rgb(hsvImg); % 矢量场参数 stride = 6; [m,n] = size(S); xq = 1:stride:n; yq = 1:stride:m; [Xq,Yq] = meshgrid(xq*dx, yq*dy); U = cos(theta(yq,xq)); V = sin(theta(yq,xq)); % 绘制主图 figure(1); clf; imagesc([0, max(X(:))], [0, max(Y(:))], rgbImg); hold on; quiver(Xq, Yq, U, V, 0.5, 'k', 'LineWidth', 1.1); axis equal tight; % 坐标轴标注 xlabel('X (μm)'); ylabel('Y (μm)'); title(sprintf('CNC自组装过程 t = %.1f s', currentTime)); % 定制化颜色条 hcb = colorbar('Position',[0.92 0.1 0.02 0.8]); colormap(hsv); hcb.Ticks = 0:0.25:1; hcb.TickLabels = {'-90°','-45°','0°','45°','90°'}; hcb.AxisLocation = 'in'; drawnow; end```此代码平面内一直显黑色，请修改后重新输出