MATLAB求解浸润角

1、概述

　　通过后处理软件获得流场数据的密度云图，做出附着于壁面上液滴或气泡在接触点上的切线，测量可以获得浸润角度，手工做切线存在较大的自由度，所得的角度往往并不精确，可以考虑用圆拟合液滴或者气泡，通过拟合出的圆半径及圆心位置，则可以精确的计算出相应的角度，相比做切向的方法，准确度有较大的提升，但是如果还是通过手工来找准这个拟合圆，则数据准确性也受影响，可以考虑采用程序，通过图像处理的方式获取相应的信息，下面是大致的实现思路：

• 从数据文件中获取流场数据
• 从密度云图图像提取相界面处密度等值线
• 通过等值线散点拟合圆，获取圆数据
• 求解浸润角

2、求解效果

下面是一个算例的求解效果图

3、代码实现

clc;clear;
close all;

% load color data
load mycolor.mat;

% extract data
file_info = dir('*.plt');
tempB = file_info.name( isstrprop( file_info.name, 'digit' ) );
num = str2double( tempB );

% file name
filename = strcat( 'result', num2str(num), '.plt' );
[var,var_name,var_num] = tecplot2mat_2D(filename);

% load variables
for j = 1:var_num
    eval([var_name{1,j},'=var(:,:,j)'';']);
end

% figure 1
figure('units','centimeters','position',[32 18 18 6])
hold on
axis equal
box on
axis off
ax = gca;

pcolor( X, Y, Density );
[ C, h ] = contour( X, Y, Density, [ 1.0 ] , 'k-', 'linewidth', 2 );
shading interp

caxis( [ 0.5, 7 ] )
colormap(mycolor);
colorbar

% setup the axis
ax.LineWidth = 1.5;
ax.FontSize = 12;
ax.FontName = 'Times New Roman';

% delete white space
set(gca,'looseinset',get(gca,'tightinset'))
set(gca,'looseinset',[0 0 0 0])

% points for polyfit circle
index = C(1,:) > min(min(X)) & C(1,:) < max(max(X));
C = C(:,index);
index = C(2,:) > min(min(Y)) & C(2,:) < max(max(Y));
C = C(:,index);
% plot( C(1,:), C(2,:), 'b.', 'markersize', 7 )

% polyfit circle
[ xc, yc, R, a ] = circfit( C(1,:), C(2,:) );
theta = 0:0.1:2*pi;
Circle1 = xc+R*cos(theta);
Circle2 = yc+R*sin(theta);

plot( xc, yc, 'k.', 'markersize', 15 );
plot( Circle1, Circle2, 'm:', 'linewidth', 2 );

% save picture
picname = strcat( './Diameter_', num2str( num ), '.tiff');
print( gcf, '-dtiff', '-r150', picname );

% contact angle
if yc >= R
    angle = 0;
else
    angle = 90 - asin( yc / R ) * 180 / pi;
end

fprintf( 'contact angle = % 10.5f 
', angle );

 子函数 circfit 

function [xc,yc,R,a] = circfit(x,y)
%CIRCFIT Fits a circle in x,y plane
% [XC, YC, R, A] = CIRCFIT(X,Y)
% Result is center point (yc,xc) and radius R.A is an
% optional output describing the circle’s equation:
% x^2+y^2+a(1)*x+a(2)*y+a(3)=0

n = length(x);
xx = x.*x;
yy = y.*y;
xy = x.*y;
A = [ sum(x), sum(y), n; sum(xy), sum(yy), sum(y); sum(xx), sum(xy), sum(x) ];
B = [ -sum(xx+yy); -sum(xx.*y+yy.*y); -sum(xx.*x+xy.*y) ];
a = AB;
xc = -.5*a(1);
yc = -.5*a(2);
R = sqrt((a(1)^2+a(2)^2)/4-a(3));
end