引言

第二类超导体在混合态下呈现出独特的磁通涡旋结构，这些涡旋以周期性晶格形式排列，被称为 Abrikosov 晶格。涡旋结构的形成机制与超导体的几何形状密切相关，这一现象在凝聚态物理研究中具有重要地位。本文通过 COMSOL Multiphysics 平台，系统模拟分析不同几何构型对 Abrikosov 涡旋晶格形态的影响规律。

一、理论框架

Abrikosov 晶格的形成源于 Ginzburg-Landau（GL）理论，该理论是描述第二类超导体混合态的基本框架。在 GL 理论中，超导态由复序参量 ψ(r) 和磁矢势 A(r) 共同描述，两者满足耦合的 GL 方程组。

系统的 GL 自由能泛可表示为：

其中 α 和 β 为与温度相关的 GL 系数，m* 为 Cooper 电子对的有効质量，ħ 为约化普朗克常数，e 为基本电荷，μ₀ 为真空磁导率。当外加磁场介于下临界场 Hc1 和上临界场 Hc2 之间时，系统能量极小化导致形成量子化的磁通涡旋，这些涡旋在弹性和磁相互作用的共同驱动下自组织形成 Abrikosov 晶格，通常呈现三角格子结构。

二、COMSOL 数值建模方法

2.1 物理场配置

在 COMSOL 环境中，我们采用 AC/DC 模块的超导子模块进行模拟。首先需要构建描述序参量和电磁场的耦合方程组。

% 创建电磁准静态分析模型
model = mphmodel('create');
model.label('Abrikosov Vortex Simulation');

% 添加磁场物理场接口
fe = addPhysics(model, 'mf', 'geometricalOrder', 'linear');

% 初始化序参量场（复数标量场）
psiField = addVariable(model, 'psi', 'domain', 'geometricDomain', 'unit', '1');
psiField.shallowDescription = 'Superconducting order parameter';
psiField.value = {'1', '0'};  % 初始幅值为1，相位为0

% 初始化矢势场（仅z分量用于二维模型）
aField = addVariable(model, 'Az', 'domain', 'geometricDomain', 'unit', 'Wb/m');
aField.shallowDescription = 'Magnetic vector potential';
aField.value = {'0'};

上述代码构建了包含复序参量和磁矢势的物理场模型框架，为后续的 GL 方程求解奠定基础。

2.2 几何参数化建模

几何形状对涡旋分布具有显著影响。我们通过参数化方式定义不同的几何构型：

% 定义矩形超导区域参数
rectParams.length = 2e-6;    % 长度 2 μm
rectParams.width = 2e-6;     % 宽度 2 μm

% 创建矩形几何
g1 = model.geom.create('rect1', 2);
g1.set('lx', rectParams.length);
g1.set('ly', rectParams.width);
g1.set('pos', [0 0]);

% 定义圆形超导区域参数
circleParams.radius = 1e-6;  % 半径 1 μm

% 创建圆形几何
g2 = model.geom.create('circle1', 2);
g2.set('r', circleParams.radius);
g2.set('pos', [0 0]);

% 定义六边形超导区域参数
hexParams.radius = 1e-6;     % 外接圆半径

% 创建六边形几何
g3 = model.geom.create('hex1', 2);
g3.set('r', hexParams.radius);
g3.set('rot', 0);

本模型采用二维平面近似，以便高效计算涡旋结构的细节分布。实际模拟中可根据需要扩展至三维构型。

2.3 边界条件与材料属性

超导体表面的边界条件对涡旋行为具有决定性影响。我们设置以下边界条件：

% 设置超导区域材料参数
mat = model.material.create('superconductor', 'domain', 'geom1');
mat.set('relativePermeability', '1');
mat.set('electricalConductivity', '0');  % 超导体电导率为零

% 序参量边界条件：在边界处设为零（超导相干长度外）
bc1 = addPhysicsFeature(fe, 'ZeroScalarBoundary', 'selection', 'boundary');
bc1.set('psi', '0');

% 磁矢势边界条件：磁绝缘边界
bc2 = addPhysicsFeature(fe, 'MagneticInsulation', 'selection', 'boundary');
bc2.set('Az', '0');

% 设置外部磁场激励（施加z方向均匀磁场）
extField = addPhysicsFeature(fe, 'MagneticField', 'selection', 'boundary');
extField.set('Hx', '0');
extField.set('Hy', '0');
extField.set('Hz', '0.01');  % 外加0.01 T磁场

上述边界条件模拟了超导体置于外部均匀磁场中的物理情景，序参量在边界处衰减至零反映了超导态在表面的 suppression 效应。

2.4 网格离散化策略

涡旋核心区域具有极高的序参量梯度，需要精细的网格划分以确保数值精度：

% 生成自适应网格
mesh = model.mesh.create('mesh1', 'geom1');

% 全局网格尺寸控制
mesh.set('hauto', 1);  % 自动选择基础网格级别

% 局部加密（涡旋核心区域）
refine1 = mesh.create('refineCore', 'edge');
refine1.set('hmax', 5e-8);  % 核心区域最大单元尺寸50 nm
refine1.set('selection', 'domain');

% 执行网格生成
mesh.run;

对于包含多个涡旋的复杂构型，建议采用自适应网格加密策略，以在保证精度的同时控制计算规模。

2.5 求解器配置与结果提取

GL 方程组的非线性特性要求采用迭代求解方法：

% 配置非线性求解器
sol = model.sol.create('sol1');
sol.set('stol', 1e-6);  % 相对残差容差
sol.set('maxiter', 500);  % 最大迭代次数

% 使用牛顿迭代法
newton = sol.create('n1', 'Newton');
newton.set('damping', 0.8);  % 阻尼因子
newton.set('maxiter', 100);

% 执行求解
sol.run;

% 提取序参量分布结果
psiResult = mphinterp(model, 'psi', 'dataset', 'sol1');
azResult = mphinterp(model, 'Az', 'dataset', 'sol1');

% 计算磁感应强度分布
bzResult = mpheval(model, 'mf.Bz', 'dataset', 'sol1');

通过上述求解过程，我们可以获得序参量 ψ 和磁感应强度 B 的空间分布，进而分析涡旋结构的形态特征。

三、几何效应的模拟结果分析

3.1 矩形边界条件下的涡旋排列

在矩形超导薄膜中，涡旋排列呈现出显著的方向性特征。数值模拟表明：

• 涡旋倾向于在矩形长边方向形成准直线排列
• 边界附近涡旋间距呈现梯度分布，角区涡旋密度最低
• 涡旋数量随外加磁场增加呈阶梯式增长

这种各向异性源于矩形边界对涡旋间相互作用力的调制效应。

3.2 圆形对称构型的涡旋环结构

圆形几何具有旋转对称性，涡旋在此构型下呈现独特的同心环状分布模式：

• 低磁场下，涡旋聚集于圆心附近形成单核心结构
• 中间磁场强度时，涡旋向外迁移形成环形壳层
• 高磁场下，涡旋分布趋于均匀，边界效应占主导

圆形构型的对称性使其成为研究涡旋-涡旋相互作用的理想模型系统。

3.3 多边形边数效应

多边形超导区域的几何对称性随边数变化，涡旋结构呈现明显的边数依赖性：

• 四边形（正方形）构型产生十字对称的涡旋阵列
• 六边形构型呈现六重旋转对称的涡旋模式
• 边数增加时，涡旋分布趋近圆形极限

多边形构型中的棱角处存在局域化的应力集中，导致涡旋优先在平直边附近排列。

结论

本文基于 COMSOL 平台的数值模拟，系统研究了不同几何形状对 Abrikosov 涡旋晶格的影响机制。研究结果表明：几何边界条件通过调制涡旋间的相互作用势能，深刻影响涡旋的平衡分布模式。矩形构型产生各向异性的涡旋链，圆形构型形成同心涡旋环，多边形构型则呈现与边数相关的对称模式。这些发现为理解超导磁通动力学以及设计新型超导电子器件提供了重要的理论依据。未来研究可进一步探索缺陷、pinning 中心以及多层异质结构对涡旋晶格的调控作用。