COMSOL——LiveLink for MATLAB学习2——随机几何的构建

上一篇学习了一下如何打开与使用Livelink for MATLAB。
https://blog.csdn.net/jessica0307/article/details/105444407
这一篇学习一下具体如何使用MATLAB脚本构建COMSOL模型。COMSOL官方给出了几个学习案例，对于自己的模型，也可以在COMSOL中另存为.m文件然后学习各命令行的具体含义（注意导出之前先压缩历史记录）。不需要一次记住所有，边用边查就行。

先从构建模型几何开始，在实践中学习。
模型描述：构建一个随机生成的颗粒填充几何，2D，粒径满足正态分布。

MATLAB 代码

function out = random_circle
% Model created on Apr 11 2020,  by COMSOL 5.5.

通用模型语句，相当于COMSOL desktop的模型引导部分，确定几何维度、物理场等

import com.comsol.model.*
import com.comsol.model.util.*

model = ModelUtil.create('Model');
model.component.create('comp1', true); % 生成组件1
model.component('comp1').geom.create('geom1', 2); % 生成2D几何
model.component('comp1').mesh.create('mesh1'); % 生成网格
model.component("comp1").geom("geom1").lengthUnit("mm");

设置全局参数,目前在字符串和数值这里有点混乱，如果直接model.param.set(‘Length’,‘1’),那么后面无法用直接调用Length。不太清楚，先这么着吧。

Length = 1;
model.param.set('Length', num2str(Length)); % 参数名称，值
model.param.descr('Length', 'The Length of a cube');% 参数描述
vf = 0.5;
model.param.set('vf', num2str(vf));
model.param.descr('vf', 'The volume fraction of fillers');
Vsq = Length^2;
model.param.set('Vsq', num2str(Vsq));
miu = 0.05;
model.param.set('miu', num2str(miu)); % 用于控制circle大小的正态分布参数
model.param.descr('miu', 'Average radius');
sigma = 0.02;
model.param.set('sigma', num2str(sigma));
model.param.descr('sigma', 'standard deviation');

接下来我们去生成几何, 先画一个正方形

model.component('comp1').geom('geom1').create('sq1', 'Square');
model.component('comp1').geom('geom1').feature('sq1').set('size', 'Length');
model.component('comp1').geom('geom1').feature('sq1').set('pos', [0 0]); % 坐标原点为默认值，此句可以省略
model.component('comp1').geom('geom1').run('sq1');

此时可以用mphgeom命令在MATLAB里看一下效果

 mphgeom(model,'geom1');

接下来去生成填充的随机圆的坐标与半径

n = 10000;
Vsum = 0;
Pos = zeros(n,2);
R = zeros(n,1);
idx = 1; % index for circle 
flag = 0;
while (Vsum < Vsq * vf)
    r = abs( normrnd(miu,sigma) ); % 随机生成cicle
    pos = [Length * rand(1,1) Length * rand(1,1)];
    for k = 1:idx %将随机生成的cicle与已存在的所有cicle进行距离判断
        Distance = sqrt((pos(1)-Pos(k,1))^2+(pos(2)-Pos(k,2))^2);
        rsum = r+R(k);
        if Distance < rsum
            flag = 1;
            break;
        end
    end
    
    if flag == 1 % 如果新生成cicle与任意一个已存在cicle重叠，则进入下一轮循环，放弃此次生成的cicle
        flag = 0;
        continue;
    end
    
    if (pos(1)-r < 0) || (pos(2)-r < 0) || (pos(1)+r > Length) || (pos(2)+r > Length)%判断是否在正方形内
        continue;
    end
    
    V = Vsum + 2 * pi * r * r;
    if V > vf * Vsq %进行体积分率条件判断
        break;
    end
    
    % 至此，随机生成的cicle参数满足不重叠条件、正方形内条件与体积分率条件判断，将其正式生成几何
    cl_name = ['cl',num2str(idx)]; % cicle序列号
    model.component('comp1').geom('geom1').create(cl_name, 'Circle');
    model.geom('geom1').feature(cl_name).set('base', 'center');
    model.component('comp1').geom('geom1').feature(cl_name).set('r', num2str(r));
    model.component('comp1').geom('geom1').feature(cl_name).set('pos', pos);
    mphgeom(model,'geom1');
    
    Pos(idx,:) = pos;
    R(idx) = r;
    idx = idx +1;
    Vsum = Vsum + 2 * pi * r * r; 
    
end
mphsave(model,'random_circle'); % 保存mph文件至当前文件夹
out = model;                

模型效果图

COMSOL Desktop

一方面代码最后一行直接生成了mph文件，可以在comsol中打开
另一方面可以通过服务器打开（文件→COMSOL Server→从服务器导入App）

今天就到这里吧。最后附上参考文献
天乐树的博客：http://blog.sina.com.cn/s/blog_4a0a8b5d0101lnwf.html
感谢天乐树！