matlab网格划分程序与matlab有限元的结合

最新推荐文章于 2024-08-05 13:52:13 发布
最新推荐文章于 2024-08-05 13:52:13 发布
阅读量5.8k 收藏 48
点赞数 6

原文来链接：https://blog.csdn.net/conjimmy/article/details/19070495
保存自用；
distmesh是一个较好的网格划分程序，具体可以参考：http://persson.berkeley.edu/distmesh/
2.matlab有限元可以参考徐荣桥的书
3.这里本人打算画一个园内包含一个椭圆的模型：
具体程序如下：

a.

网格划分：

fh=@§ 0.05+0.3*dellipse(p,[0.5,0.2]);
fd=@§ ddiff(dcircle(p,0,0,1),dellipse(p,[0.5,0.2]));
[p,t]=distmesh2d(fd,fh,0.05,[-1,-1;1,1],[-1,-1;-1,1;1,-1;1,1]);

b.在distmesh2d.m最后插入语句，导出数据格式（节点坐标，单元材料属性，边界约束条件）

fid = fopen(‘exam4_2.dat3’, ‘w’) ;
[ilength,jlength]=size( p );
fprintf(fid,’%d\n’,ilength);
for i=1:1:ilength
fprintf(fid,’%d\t%f\t%f\n’,i,p(i,1),p(i,2));
end
[ilength,jlength]=size( t );
fprintf(fid,’%d\n’,ilength);
for i=1:1:ilength
fprintf(fid,’%d\t%d\t%d\t%d\t%d\n’,i,t(i,1),t(i,2),t(i,3),1);
end
fclose(fid);
文件’exam4_2.dat3内容如下：

5105（节点个数）
1 -0.707106 -0.707106 （每个节点坐标）
2 -0.707106 0.707107
…

9869(三角单元个数)
1 5006 4951 4934 1 （标号，i，j，k，材料属性的编号）
2 197 147 100 1
3 413 2 366 1
4 413 2 323 1

…

5(5种材料属性)
1 2.00e9 0.25 100.0 22.0e3(编号，杨氏模量，泊松比，厚度[平面应力填1]，密度)
2 2.60e9 0.20 1.0 23.0e3
3 2.85e10 0.20 1.0 25.0e3
4 1.85e10 0.20 1.0 23.0e3
5 2.85e10 0.20 1.0 22.0e3
2（约束个数）
1 2 0.0 （节点，方向【2，为y方向】，0.0（0为固定））
1 1 0.0

4 2 0.0 （节点，方向【2，为y方向】，0.0（0为固定））
41 0.0

c。调用有限元程序exam_2.m(见徐荣桥书)

exam4_2 exam4_2.dat3

d.后处理

exam4_2_post

最大主应力云图如下：

最大剪应力云图如下：

最小主应力云图如下：

exam_2。m源程序如下：

function exam4_2( file_in )
% 本程序为第四章的第二个算例，采用三角形单元计算隧道在自重作用下的变形和应力
% exam4_2( filename )
% 输入参数：
% file_in ---------- 有限元模型文件

% 定义全局变量
% gNode ------------- 节点坐标
% gElement ---------- 单元定义
% gMaterial --------- 材料性质
% gBC1 -------------- 第一类约束条件
% gK ---------------- 整体刚度矩阵
% gDelta ------------ 整体节点坐标
% gNodeStress ------- 节点应力
% gElementStress ---- 单元应力
global gNode gElement gMaterial gBC1 gK gDelta gNodeStress gElementStress gNF

if nargin < 1
    file_in = 'exam4_2.dat' ;
end

% 检查文件是否存在
if exist( file_in ) == 0
    disp( sprintf( '错误：文件 %s 不存在', file_in ) )
    disp( sprintf( '程序终止' ) )
    return ;
end

% 根据输入文件名称生成输出文件名称
[path_str,name_str,ext_str] = fileparts( file_in ) ;
ext_str_out = '.mat' ;
file_out = fullfile( path_str, [name_str, ext_str_out] ) ;

% 检查输出文件是否存在
if exist( file_out ) ~= 0 
    answer = input( sprintf( '文件 %s 已经存在，是否覆盖? ( Yes / [No] ):  ', file_out ), 's' ) ;
    if length( answer ) == 0
        answer = 'no' ;
    end
    
    answer = lower( answer ) ;
    if answer ~= 'y' | answer ~= 'yes' 
        disp( sprintf( '请使用另外的文件名，或备份已有的文件' ) ) ;
        disp( sprintf( '程序终止' ) );
        return ; 
    end
end

% 建立有限元模型并求解，保存结果
FemModel( file_in ) ;          % 定义有限元模型
SolveModel ;                   % 求解有限元模型
SaveResults( file_out ) ;      % 保存计算结果

% 计算结束
disp( sprintf( '计算正常结束，结果保存在文件 %s 中', file_out ) ) ;
disp( sprintf( '可以使用后处理程序 exam4_2_post.m 显示计算结果' ) ) ;

return ;

function FemModel(filename)
% 定义有限元模型
% 输入参数：
% filename — 有限元模型文件
% 返回值：
% 无
% 说明：
% 该函数定义平面问题的有限元模型数据：
% gNode ------- 节点定义
% gElement ---- 单元定义
% gMaterial — 材料定义，包括弹性模量，梁的截面积和梁的抗弯惯性矩
% gBC1 -------- 约束条件

global gNode gElement gMaterial gBC1 gNF

% 打开文件
fid = fopen( filename, 'r' ) ;

% 读取节点坐标
node_number = fscanf( fid, '%d', 1 ) ;
gNode = zeros( node_number, 2 ) ;
for i=1:node_number
    dummy = fscanf( fid, '%d', 1 ) ;
    gNode( i, : ) = fscanf( fid, '%f', [1, 2] ) ;
end

% 读取单元定义
element_number = fscanf( fid, '%d', 1 ) ;
gElement = zeros( element_number, 4 ) ;
for i=1:element_number
    dummy = fscanf( fid, '%d', 1 ) ;
    gElement( i, : ) = fscanf( fid, '%d', [1, 4] ) ;
end

% 读取材料信息
material_number = fscanf( fid, '%d', 1 ) ;
gMaterial = zeros( material_number, 4 ) ;
for i=1:material_number
    dummy = fscanf( fid, '%d', 1 ) ;
    gMaterial( i, : ) = fscanf( fid, '%f', [1,4] ) ;
end

% 读取边界条件
bc1_number = fscanf( fid, '%d', 1 ) ;
gBC1 = zeros( bc1_number, 3 ) ;
for i=1:bc1_number
    gBC1( i, 1 ) = fscanf( fid, '%d', 1 ) ;
    gBC1( i, 2 ) = fscanf( fid, '%d', 1 ) ;
    gBC1( i, 3 ) = fscanf( fid, '%f', 1 ) ;
end

% 关闭文件
fclose( fid ) ;
 % 集中力
%     节点号   自由度号   集中力值
gNF = [  1,       1,         -800e3;
   1,       2,         -800e3;
   4,       1,         -800e3;
   4,       2,         -800e3;
       ] ;

return

function SolveModel
% 求解有限元模型
% 输入参数：
% 无
% 返回值：
% 无
% 说明：
% 该函数求解有限元模型，过程如下
% 1. 计算单元刚度矩阵，集成整体刚度矩阵
% 2. 计算单元的等效节点力，集成整体节点力向量
% 3. 处理约束条件，修改整体刚度矩阵和节点力向量
% 4. 求解方程组，得到整体节点位移向量
% 5. 计算单元应力和节点应力

global gNode gElement gMaterial gBC1 gK gDelta gNodeStress gElementStress gNF

% step1. 定义整体刚度矩阵和节点力向量
[node_number,dummy] = size( gNode ) ;
gK = sparse( node_number * 2, node_number * 2 ) ;
f = sparse( node_number * 2, 1 ) ;

% step2. 计算单元刚度矩阵，并集成到整体刚度矩阵中
[element_number,dummy] = size( gElement ) ;
for ie=1:1:element_number
    disp( sprintf(  '计算刚度矩阵，当前单元: %d', ie  ) ) ;
    k = StiffnessMatrix( ie ) ;
    AssembleStiffnessMatrix( ie, k ) ;
end

% step3. 把集中力直接集成到整体节点力向量中
[nf_number, dummy] = size( gNF ) ;
for inf=1:1:nf_number
n = gNF( inf, 1 ) ;

 d = gNF( inf, 2 ) ;
   f( (n-1)*2 + d ) = gNF( inf, 3 ) ;

end

% step3. 计算自重产生的等效节点力

% for ie=1:1:element_number
% disp( sprintf( ‘计算自重的等效节点力，当前单元: %d’, ie ) ) ;
% egf = EquivalentGravityForce( ie ) ;
% i = gElement( ie, 1 ) ;
% j = gElement( ie, 2 ) ;
% m = gElement( ie, 3 ) ;
% f( (i-1)*2+1 : (i-1)*2+2 ) = f( (i-1)*2+1 : (i-1)*2+2 ) + egf( 1:2 ) ;
% f( (j-1)*2+1 : (j-1)*2+2 ) = f( (j-1)*2+1 : (j-1)*2+2 ) + egf( 3:4 ) ;
% f( (m-1)*2+1 : (m-1)*2+2 ) = f( (m-1)*2+1 : (m-1)*2+2 ) + egf( 5:6 ) ;
% end

% step4. 处理约束条件，修改刚度矩阵和节点力向量。采用乘大数法
[bc_number,dummy] = size( gBC1 ) ;
for ibc=1:1:bc_number
    n = gBC1(ibc, 1 ) ;
    d = gBC1(ibc, 2 ) ;
    m = (n-1)*2 + d ;
    f(m) = gBC1(ibc, 3)* gK(m,m) * 1e15 ;
    gK(m,m) = gK(m,m) * 1e15 ;
end

% step 5. 求解方程组，得到节点位移向量
gDelta = gK \ f ;

% step 6. 计算单元应力
gElementStress = zeros( element_number, 6 ) ;
for ie=1:element_number 
    disp( sprintf(  '计算单元应力，当前单元: %d', ie  ) ) ;
    es = ElementStress( ie ) ;
    gElementStress( ie, : ) = es ;
end

% step 7. 计算节点应力(采用绕节点加权平均)
gNodeStress = zeros( node_number, 6 ) ;       
for i=1:node_number                         
    disp( sprintf(  '计算节点应力，当前结点: %d', i  ) ) ;
    S = zeros( 1, 3 ) ;                         
    A = 0 ;
    for ie=1:1:element_number
        for k=1:1:3
            if i == gElement( ie, k ) 
                area= ElementArea( ie ) ;
                S = S + gElementStress(ie,1:3 ) * area ;
                A = A + area ;
                break ;
            end
        end
    end
    gNodeStress(i,1:3) = S / A ;
    gNodeStress(i,6) = 0.5*sqrt( (gNodeStress(i,1)-gNodeStress(i,2))^2 + 4*gNodeStress(i,3)^2 ) ;
    gNodeStress(i,4) = 0.5*(gNodeStress(i,1)+gNodeStress(i,2)) + gNodeStress(i,6) ;
    gNodeStress(i,5) = 0.5*(gNodeStress(i,1)+gNodeStress(i,2)) - gNodeStress(i,6) ;
end

return

function k = StiffnessMatrix( ie )
% 计算单元刚度矩阵
% 输入参数:
% ie ---- 单元号
% 返回值:
% k ---- 单元刚度矩阵

global gNode gElement gMaterial 
k = zeros( 6, 6 ) ;
E  = gMaterial( gElement(ie, 4), 1 ) ;
mu = gMaterial( gElement(ie, 4), 2 ) ;
h  = gMaterial( gElement(ie, 4), 3 ) ;
xi = gNode( gElement( ie, 1 ), 1 ) ;
yi = gNode( gElement( ie, 1 ), 2 ) ;
xj = gNode( gElement( ie, 2 ), 1 ) ;
yj = gNode( gElement( ie, 2 ), 2 ) ;
xm = gNode( gElement( ie, 3 ), 1 ) ;
ym = gNode( gElement( ie, 3 ), 2 ) ;
ai = xj*ym - xm*yj ;
aj = xm*yi - xi*ym ;
am = xi*yj - xj*yi ;
bi = yj - ym ;
bj = ym - yi ;
bm = yi - yj ;
ci = -(xj-xm) ;
cj = -(xm-xi) ;
cm = -(xi-xj) ;
area = (ai+aj+am)/2 ;
B = [bi  0 bj  0 bm  0
      0 ci  0 cj  0 cm
     ci bi cj bj cm bm] ;
B = B/2/area ;
D = [ 1-mu    mu      0
       mu    1-mu     0
        0      0   (1-2*mu)/2] ;
D = D*E/(1-2*mu)/(1+mu) ;
k = transpose(B)*D*B*h*abs(area) ;

return

function B = MatrixB( ie )
% 计算单元的应变矩阵B
% 输入参数:
% ie ---- 单元号
% 返回值:
% B ---- 单元应变矩阵
global gNode gElement
xi = gNode( gElement( ie, 1 ), 1 ) ;
yi = gNode( gElement( ie, 1 ), 2 ) ;
xj = gNode( gElement( ie, 2 ), 1 ) ;
yj = gNode( gElement( ie, 2 ), 2 ) ;
xm = gNode( gElement( ie, 3 ), 1 ) ;
ym = gNode( gElement( ie, 3 ), 2 ) ;
ai = xjym - xmyj ;
aj = xmyi - xiym ;
am = xiyj - xjyi ;
bi = yj - ym ;
bj = ym - yi ;
bm = yi - yj ;
ci = -(xj-xm) ;
cj = -(xm-xi) ;
cm = -(xi-xj) ;
area = (ai+aj+am)/2 ;
B = [bi 0 bj 0 bm 0
0 ci 0 cj 0 cm
ci bi cj bj cm bm] ;
B = B/2/area ;
return

function area = ElementArea( ie )
% 计算单元面积
% 输入参数:
% ie ---- 单元号
% 返回值:
% area ---- 单元面积
global gNode gElement gMaterial

xi = gNode( gElement( ie, 1 ), 1 ) ;
yi = gNode( gElement( ie, 1 ), 2 ) ;
xj = gNode( gElement( ie, 2 ), 1 ) ;
yj = gNode( gElement( ie, 2 ), 2 ) ;
xm = gNode( gElement( ie, 3 ), 1 ) ;
ym = gNode( gElement( ie, 3 ), 2 ) ;
ai = xj*ym - xm*yj ;
aj = xm*yi - xi*ym ;
am = xi*yj - xj*yi ;
area = abs((ai+aj+am)/2) ;

return

function AssembleStiffnessMatrix( ie, k )
% 把单元刚度矩阵集成到整体刚度矩阵
% 输入参数:
% ie — 单元号
% k — 单元刚度矩阵
% 返回值:
% 无
global gElement gK
for i=1:1:3
for j=1:1:3
for p=1:1:2
for q=1:1:2
m = (i-1)*2+p ;
n = (j-1)*2+q ;
M = (gElement(ie,i)-1)*2+p ;
N = (gElement(ie,j)-1)*2+q ;
gK(M,N) = gK(M,N) + k(m,n) ;
end
end
end
end
return

function egf = EquivalentGravityForce( ie )
% 计算单元自重的等效节点力
% 输入参数
% ie ----- 单元号
% 返回值
% egf ----- 自重的等效节点力
global gElement gMaterial

area = ElementArea( ie ) ;
h    = gMaterial( gElement( ie, 4 ), 3 ) ;
ro   = gMaterial( gElement( ie, 4 ), 4 ) ;
w    = area * h * ro ;
egf  = -w/3 * [0; 1; 0; 1; 0; 1] ;

return

function es = ElementStress( ie )
% 计算单元的应力分量
% 输入参数
% ie ----- 单元号
% 返回值
% es ----- 单元应力分量列阵（1×6）： [sx, sy, txy, s1, s2, tmax]
global gElement gDelta gMaterial

es = zeros( 1, 6 ) ;   % 单元的应力分量
de = zeros( 6, 1 ) ;   % 单元节点位移列阵
E  = gMaterial( gElement(ie, 4), 1 ) ;
mu = gMaterial( gElement(ie, 4), 2 ) ;
D = [ 1-mu    mu      0
       mu    1-mu     0
        0      0   (1-2*mu)/2] ;
D = D*E/(1-2*mu)/(1+mu) ;
B = MatrixB( ie ) ;
for j=1:1:3
    de( 2*j-1 ) = gDelta( 2*gElement( ie, j )-1 ) ;
    de( 2*j   ) = gDelta( 2*gElement( ie, j )   ) ;
end
es(1:3) = D * B * de ;
es(6) = 0.5*sqrt((es(1)-es(2))^2 + 4*es(3)^2 ) ;
es(4) = 0.5*(es(1)+es(2)) + es(6) ;
es(5) = 0.5*(es(1)+es(2)) - es(6) ;

return

function SaveResults( file_out )
% 显示计算结果
% 输入参数：
% file_out — 保存结果文件
% 返回值：
% 无

global gNode gElement gMaterial gBC1 gDelta gNodeStress gElementStress
save( file_out, 'gNode', 'gElement', 'gMaterial', 'gBC1', 'gDelta', 'gNodeStress', 'gElementStress' ) ;

return

未来，我来～
关注
用matlab表示e有限元,matlab网格划分程序与matlab有限元的结合
weixin_35682992的博客
03-17 892
1.2.matlab有限元可以参考徐荣桥的书3.这里本人打算画一个园内包含一个椭圆的模型：具体程序如下：a.网格划分：>>fh=@(p) 0.05+0.3*dellipse(p,[0.5,0.2]);>> fd=@(p) ddiff(dcircle(p,0,0,1),dellipse(p,[0.5,0.2]));>> [p,t]=distmesh2d(fd,fh...
Matlab网格划分
weixin_30783913的博客
09-26 5494
之前转载了一篇博客http://blog.sina.com.cn/s/blog_6163bdeb0102dvay.html，讲Matlab网格划分程序Distmesh，看了看程序，感觉程序写得有很多值得学的地方，所以又自己重新看了一看，加了一些注释，最后再总结一下学到的东西吧。 源代码的地址已经改变，如下http://people.sc.fsu.edu/~jburkardt/m_src/dist...
【物理应用】基于有限元求解平面弹性问题matlab实现
最新发布
matlab_dingdang的博客
08-05 1017
有限元法 (FEM) 作为一种强大的数值计算方法，在工程领域得到了广泛应用，尤其在解决结构力学问题方面表现突出。本文将深入探讨基于有限元法的平面弹性问题求解，并利用 MATLAB 语言进行编程实现。文章首先阐述平面弹性问题的基本理论，包括应力应变关系、平衡方程以及边界条件；接着介绍有限元法的基本原理，包括单元划分、形函数、刚度矩阵和节点力等概念；最后，以一个具体案例为例，展示基于 MATLAB 实现平面弹性问题的求解过程，包括模型建立、单元划分、矩阵组装、求解线性方程组以及结果后处理等环节。
MATLAB网格划分
07-24
用MATLAB语言划分网格，计算流体的基础
matlabFEM.zip_MATLAB 有限元_matlab程序FEM_matlab网格划分_有限元Matlab_有限元网格
07-14
matlab有限元计算实例，包括网格划分，有限元分析，数据处理，可视化等
3_1.rar_Matlab网格划分工具_dist_有限元 matlab_有限元 分网_有限元网格
07-15
matlab网格划分工具，有限元分析前处理，分网
三角形网格matlab有限元,matlab有限元网格划分,matlab
09-10
**三、MATLAB有限元网格划分** 在MATLAB中，实现有限元网格划分通常包括以下步骤： 1. **定义几何区域**：根据COORX和COORY数据文件定义模型的几何边界。 2. **生成元素连接**：通过算法将节点连接成三角形元素，...
grid.rar_matlab网格划分_流体力学MATLAB_流体有限元
07-15
MATLAB作为一款强大的数值计算工具，常被用于进行网格划分，特别是对于流体有限元分析。本资源"grid.rar"包含两个MATLAB程序，`generate_triangle_grid.m`和`triangle_grid.m`，旨在帮助用户理解和应用网格划分技术...
DistMesh - A Simple Mesh Generator in MATLAB （有限元等网格剖分）
06-07
DistMesh, A Simple Mesh Generator in MATLAB, 有限元网格剖分不错工具，可做项目测试用，二维三维都可剖分，其主页介绍：http://persson.berkeley.edu/distmesh/
matlab网格划分
06-08
matlab网格划分，互相学习，但愿有帮助
grid.rar_grid forces_matlab划分网格_网格_网格划分_网格划分 MATLAB
07-13
有限元分析与应用中的程序，用于网格划分，和单个网格分析
matlab naca0012网格划分
03-27
matlab对naca0012翼型的网格划分文件，其中interfunction是翼型形状函数
Matlab中二维网格划分
weixin_42575505的博客
02-18 1224
Matlab可以使用meshgrid函数来实现二维网格划分，具体操作步骤如下：1.输入两个变量的值域，例如 x= linspace(-1,1,5)；y = linspace(-1,1,7)；2.调用meshgrid函数，[X,Y] = meshgrid(x,y)；3.可以使用surf函数绘制网格，surf(X,Y,Z)；4.最后使用mesh函数绘制网格，mesh(X,Y,Z)。 ...
matlab网格自动划分,有限元实现三维自动划分网格.m
weixin_36304352的博客
03-21 3306
% http://myangie.blog.sohu.com/6192144.html% 在开发有限元程序的时候，% 首先遇到的问题就是实现自动划分网格(meshgrid)，% matlab可以实现这一功能，它的基本原理就是% 根据图论和邻接矩阵的原理利用gplot函数和sparse稀疏矩阵实现的，% 但是它只能实现二维！本人重新编写了一个程序，% 可以实现三维的网格自动划分% (包括三角形单元、...
【matlab教程】20、简单网格细分
qq_38939822的博客
11-12 3901
问题 matlab处理网格时，有时需要将网格细分 思路参考：https://hanspond.github.io/2018/11/27/%E7%AE%80%E5%8D%95%E7%BD%91%E6%A0%BC%E7%BB%86%E5%88%86%201to4%20Mesh%20Subdivision/index.html 将每一个面的每一条边取中点并连接，形成新的三角形，原来一个面细分成四个面 原博主的代码只适合封闭曲面的细分，如果曲面是不封闭的，则拓扑结构会乱，如： 代码： %% created by
Matlab网格划分程序Distmesh讲解
06-01
Distmesh是MATLAB中的一种网格划分程序，可以用于生成各种形状的二维和三维网格。Distmesh的主要优点是可以方便地控制网格的密度和形状，以及能够自适应地调整网格的大小和形状以适应特定的几何形状。 Distmesh的基本原理是将待划分的区域分割成一系列小区域，并在每个小区域内生成一个网格节点。然后，使用一些算法（如Delaunay三角剖分）将这些节点连接起来，形成网格。 以下是使用Distmesh进行网格划分的基本步骤： 1. 定义几何形状：定义一个表示待划分区域边界的函数，该函数返回一个nx2或nx3的矩阵，其中n表示边界点的数量，每一行表示一个点的坐标。 2. 定义密度函数：定义一个表示网格密度的函数，该函数接受一个nx2或nx3的矩阵作为输入，返回一个n维向量，表示每个点的密度。 3. 运行Distmesh程序：使用定义的几何形状和密度函数作为输入，运行Distmesh程序生成网格。程序会自动调整网格大小和形状以适应几何形状和密度函数。 4. 可选：对生成的网格进行后处理，如去除不必要的节点、平滑网格等。 以下是使用Distmesh进行网格划分的示例代码： ```matlab % 定义几何形状 fd=@(p) drectangle(p,-1,1,-1,1); [p,e,t]=initmesh(fd,'Box','Hmax',0.2); % 定义密度函数 fh=@(p) 0.05+0.3*(sin(5*p(:,1)).*sin(5*p(:,2))); % 运行Distmesh程序 [p,e,t]=distmesh2d(fd,fh,0.2,[-1,-1;1,1],[]); ``` 该代码生成一个边长为2，中心坐标为(0,0)的正方形区域，并在正方形内部生成一个波浪形状的密度函数。运行Distmesh程序后，程序会自适应地生成一个网格，其中节点的密度与密度函数呈正比例关系，节点分布在波浪形状密度函数的高密度区域。 Distmesh还提供了许多其他的功能和选项，例如控制边界点的位置、指定网格的最大和最小尺寸等，可以根据具体需要进行调整。
评论 2
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值