# BEM_GUI：帮助文档
BEM_ALL用户手册
V1.3
开发环境 Windows +Qtc+ MinGW
Copyright © 2022 SJTU-NAOCE. All Rights Reserved.
说明
本手册是对BEM_ALL.exe以及BEM_GUI.exe程序的一个说明文档。
- 2023-03-30 —— V1.3：
- 添加简要地描述了程序相关输入、输出文件的必要说明
- 2023-03-20 —— V1.3：
- 添加工具函数部分
- 2023-03-18 —— V1.3：
- 添加测试案例说明
目录
目录
程序简介
主要功能
频域求解程序：FDC(Frequency Domain Code)
- 经过算例测试，程序能够正确完成单体、多体频域计算内容
- 支持添加
Internal Lid、Damping Lid模型用于限制不规则频率与抑制间隙共振 - 支持添加多个Body，Body的数目不设置上限
- Body可以设置为Floating、Fixed、Internal Lid、Damping Lid四种类型
- 支持三角形和四边形面元文件
- 总面元数量受限于计算机运行内存
- 例如：8000面元，所需运行内存约 10 G
- OpenMP多线程并行
- 例如：默认4线程，可提高3倍计算效率
- 支持在Lua脚本中自定义函数
- 添加Wall Damping Model，该模型可用于抑制间隙共振而不需要添加阻尼盖面元
- 支持自定义WallDampingRang以及OffBodyPoints功能
- 对于FloatBody与FixeBody类型，单独输出FK力
- 对于PanelModel可以支持msh格式，pan格式文件
时域求解程序：TDC(Time Domain Code)
- 基于间接时域法，多浮体时域模拟
- 支持添加多浮体之间连接结构：连接缆绳，防撞垫
- 当前版本支持添加系泊系统(集中质量法)
- 支持添加自定义外力
工具函数：可以通过命令行指定参数实现
- 可以单独对单个面元文件(*.gdf,*.msh,*.pan)进行缩放、平移、对称等处理
- 可以单独对面元文件(*.gdf,*.msh,*.pan)进行法相检查并矫正(还没写完,暂时没用)
- 可以单独对面元文件(*.gdf,*.msh,*.pan)计算静力学参数
- 以上执行顺序：
correctNormal->scale->translate->symmetry->calcHydrostaticInfo - 转换面元文件格式 [*.gdf,*.msh,*.pan] -> [*.gdf,*.pan]
其他说明：
- FDC 至多设置一个Floating Body
- FDC 当前版本不支持二阶力计算
框架结构
本程序提供的频域求解模块（FDC：Frequency Domain Code）与间接时域模块（TDC：Time Domain Code）程序，下图分别展示了 FDC 与 TDC 程序的分层框架结构。清晰的层次结构有助于进行
程序后续调试与扩展，本层的模型可以利用下层模型提供的数据与函数接口，实现自身功能，并为上层模型提供函数据与函数调用接口。
BEM_ALL.exe 使用帮助
如何运行BEM_ALL.exe软件
目前提供两种版本的程序，分别是
BEM_XXX.exe与BEM_XXX_OpenBLAS.exe，两者主要的区别在于后者利用OpenBLAS库作为矩阵运算后端程序，后者计算效率更高，但并未对其进行全面测试。
运行环境
BEM_ALL.exe软件可在Windows系统上独立运行。
注意：如果需要在Linux系统上运行，则需要添加
#define LINUX宏，重新编译，所使用的第三方库也需要重新编译。
查看版本信息
- 双击程序(相关的版本权信息和必要的说明内容如下图所示)
获得使用帮助
-
在当前路径打开cmd窗口，或者其他方法自行百度
-
cmd窗口输入
BEM_ALL_V2.0.exe -h，显示以下使用帮助 -
其中显示了全部的命令行参数，常用的参数如下
参数 含义 -jor--json指定输入的配置文件(*.json) -nor--nthread指定并行线程数(如果不指定，默认为4) -cor--check添加此参数，则仅执行检查，不执行计算 -lor--log指定输出的log文件名称(默认为default.log) -eor--example当前路径下，输出模板文件 全部参数
USAGE:
BEM_ALL_V2.0.exe [-?|-h|--help] [-n|--nthread <num of threads>] [-j|--json <json input config name>] [-c|--check] [-e|--example] [-l|--log <generate log File>] [-i|--inputfile <input "*.gdf">] [-o|--outputfile <output "*.gdf *.pan">] [--calcHydrostatic] [--correctNormal] [--scale <set the scale factor,default is 1.0>] [--translate <set the translate vector,default is 0,0,0>] [--symmetry <set the symmetry type, 0 = NULL;1 = AlongX;2 = AlongY;4 = AlongZ;3 = AlongXY;default = 0;>] [--hos] [--moor] [--valid_moor]
Display usage information.
OPTIONS, ARGUMENTS:
-?, -h, --help
-n, --nthread <num of threads>
Specifies the number of parallel threads, default n=4
-j, --json <json input config name>
Specifies the name of json config file, defalue j="default.json"
-c, --check Only check the calc info
-e, --example Generate FD/TD example [ default_FD/TD.json ] file
-l, --log <generate log File>
Specific log file name, default="default.log"
-i, --inputfile <input "*.gdf">
input panel file(input "*.gdf")
-o, --outputfile <output "*.gdf *.pan">
output panel file(output "*.gdf,*.pan")
--calcHydrostatic calculate hydrostatic info for panel
--correctNormal correct Normal for the panel
--scale <set the scale factor,default is 1.0>
--translate <set the translate vector,default is 0,0,0>
--symmetry <set the symmetry type, 0 = NULL;1 = AlongX;2 = AlongY;4 = AlongZ;3 = AlongXY;default = 0;>
--hos HOS execute( only test )
--moor Moor Dyn Lib( only test )
--valid_moor Valid for MoorDynLib(single line). Input File(top_motion.txt,line.cfg)
生成模板文件
cmd窗口输入 BEM_ALL_V2.0.exe -e，则在当前路径下生成以下文件
default_FD.json
default_TD.json
UDF.lua
lines.cfg
相关文件具体解释可参考：
如何运行程序
使用CMD窗口
在cmd窗口输入以下命令：
BEM_ALL_V2.0.exe -n 8 -j test1.json -l test1_folder/test.log
以8个线程执行test1.json算例文件，并将显示信息输出至test.log，
注意：test1.json文件中会指定输出文件夹路径，如果使用命令行方式运行，则此文件需要自行创建，如果使用BEM_GUI界面，则不需要自行创建。
使用可视化界面
输入文件有哪些?
目前，软件输入文件主要有4类：
配置文件以json格式文件给出，模板文件可以通过指定--example参数输出模板文件，参考生成模板文件。这里我们是用BEM_GUI中项目树窗口查看json文件，如下图所示。
注意：若关键词存在
(*)标记，则表示该关键字下存在必填信息，若不存在，则表示该关键字对应的数值由软件自动生成。
频域计算配置模板文件(default_FD.json)
频域计算的配置文件包含五个部分，同上图。
设定算例基本信息
"1.Basic Info(*)": {
"1.Name(*)": "object",
"2.Solver Type(*)": {
"info": "0=FD_Single;1=FD_MultiBody;2=TD_Single;3=TD_MultiBody",
"value": 0
},
"3.Project Path(*)": "Project/Path"
}
项目名称
| 关键字 | 类型 | 数值 | 说明 |
|---|---|---|---|
"1.Name(*)" |
字符串(string) | "object" |
为算例指定一个名称(建议以Solver_为前缀) |
注意：部分结果文件将以此为文件名前缀，例如：offBodyPoints的输出文件便是以此名称为前缀，并且在路径下会额外生成一个以此为名称命名的json文件，其中为当前算例输出的频域配置文件
求解器类型
在 "2.Solver Type(*)"->"value"关键字下指定不同的求解器类型。
| 关键字 | 类型 | 数值 | 说明 |
|---|---|---|---|
"2.Solver Type(*)" |
整数(int) | 0 |
单浮体频域计算 |
"2.Solver Type(*)" |
整数(int) | 1 |
多浮体频域计算 |
"2.Solver Type(*)" |
整数(int) | 2 |
单浮体时域计算 |
"2.Solver Type(*)" |
整数(int) | 3 |
多浮体时域计算 |
项目路径
| 关键字 | 类型 | 数值 | 说明 |
|---|---|---|---|
"3.Project Path(*)" |
字符串(string) | "Project/Path" |
指定文件路径 |
注意：对应文件夹需要自行创建，可以使用相对路径或者绝对路径，路径必须为全英文字符，注意
\\与/的区别，例如：./case_folder .\\case_folder E:\\FDC\\case_folder E:/FDC/case_folder
设定计算工况
"2.Condition(*)": {
"1.Rho(*)": 1000.0,
"2.Water Depth(*)": {
"info": "-1 for infinite water depth",
"value": -1.0
},
"3.Frequency(*)": "1.000000:1.000000:2.000000",
"4.Direction(*)": "0.000000:45.000000:90.000000"
}
| 关键字 | 类型 | 数值 | 说明 |
|---|---|---|---|
"1.Rho(*)" |
浮点数(double) | 1000.0 |
流体密度，默认为1000\(kg/m^2\) |
"2.Water Depth(*)" |
浮点数(double) | -1.0 |
水深条件，若为无限水深，则设定为-1.0(实际计算时对应为10000m) |
"3.Frequency(*)" |
字符串(string) | "1.00:1.00:2.00" |
最小值:步长:最大值，单位为rad/s |
"4.Direction(*)" |
字符串(string) | "0.00:45.00:90.00" |
最小值:步长:最大值，单位为deg |
可以在log文件中查看相关参数，例如：
|| INFO (2023.03.14 17:33:06) || *********************************************************
|| INFO (2023.03.14 17:33:06) || ************* FDC::Condition::showInfo() ****************
|| INFO (2023.03.14 17:33:06) || *********************************************************
|| INFO (2023.03.14 17:33:06) || G =9.81
|| INFO (2023.03.14 17:33:06) || Rho =1000
|| INFO (2023.03.14 17:33:06) || Water Depth = 10000
|| INFO (2023.03.14 17:33:06) ||
|| INFO (2023.03.14 17:33:06) || Frequency Set(rad/s). count = 2
|| INFO (2023.03.14 17:33:06) || 1.000 2.000
|| INFO (2023.03.14 17:33:06) ||
|| INFO (2023.03.14 17:33:06) || Direction Set(deg). count = 3
|| INFO (2023.03.14 17:33:06) || 0.000 45.000 90.000
|| INFO (2023.03.14 17:33:06) ||
设定计算控制参数
"3.Solver Ctrl Info(*)": {
"1.Lua UDF(*)": {
"info": "the path to lua script",
"value": "UDF.lua"
},
"2.Wall Damping Model(*)": {
"1.UDF Func Name": "",
"2.Range": {
"info": "[[x_min,y_min,z_min],[x_max,y_max,z_max]]",
"value": "[[0.000000,0.000000,0.000000],[0.000000,0.000000,0.000000]]"
},
"3.Range From UDF": {
"UDF Func Name": "",
"info": "define the range in UDF;in range \"return true\",else \"return false\";"
}
}
}
Lua UDF文件定义
| 关键字 | 类型 | 数值 | 说明 |
|---|---|---|---|
"1.Lua UDF(*)" |
字符串(string) | json |
设定Lua UDF文件路径，一般放在当前目录下，Lua中写入自定义函数，若是value为空字字符串，则不使用UDF功能，在[log文件: Brief Summary](#Brief Summary Check)中有体现。 |
"2.Wall Damping Model(*)" |
字符串(string) | json |
定义壁面阻尼模型[Wall Damping Model](#使用Wall Damping Model)相关参数 |
使用Wall Damping Model
相关算例：[test09: Barge_Barge Wall_Damping(使用壁面阻尼模型)](#test09: Barge_Barge Wall_Damping(使用壁面阻尼模型))
Wall Damping Model相关参数如下：
"2.Wall Damping Model(*)": {
"1.UDF Func Name": "",
"2.Range": {
"info": "[[x_min,y_min,z_min],[x_max,y_max,z_max]]",
"value": "[[0.000000,0.000000,0.000000],[0.000000,0.000000,0.000000]]"
},
"3.Range From UDF": {
"UDF Func Name": "",
"info": "define the range in UDF;in range \"return true\",else \"return false\";"
}
}
注意：
Wall Damping Model当且仅当定义"1.Lua UDF(*)"时才会有效，在[查看log文件信息](#Brief Summary Check)中存在提示信息，查看计算中是否使用Wall Damping Model。- 若计算中使用
Wall Damping Model，则在项目文件路径下生成Wall_Damping_Panel.gdf与Wall_Damping_Panel.txt文件，用于检查设定壁面阻尼区域。
| 关键字 | 类型 | 说明 |
|---|---|---|
"1.UDF Func Name": |
字符串(string) | 指定壁面阻尼系数分布函数(UDF_Wall_Damping_mu(x, y, z, omega))，该函数存在于上一节中指定的lua文件中)，注意：函数参数 (x, y, z, omega)形式固定，函数名UDF_Wall_Damping_mu可自定义，如例1所示。 |
"2.Range" |
字符串(string) | 定义一个三维空间区域："[[x_min,y_min,z_min],[x_max,y_max,z_max]]"，在此区域内的面元将会按照UDF_Wall_Damping_mu计算阻尼系数。 |
"3.Range From UDF" |
字符串(string) | 在UDF中定义函数，此方法作为"2.Range"的扩充，可以在UDF中编写更复杂的自定义函数，例2与例3给出了两个函数实例。 |
设定流场监测点
设定流场监测点，可以监测点位置的速度势以及波面升高，对应的配置文件如下所示。
"4.Off Body Points(*)": {
"1.Off Body Points Num": 0,
"2.Off Body Points filepath(*)": "",
"3.Get Points From UDF": {
"1.UDF Name": "",
"info": "The program will preferentially load points from Lua. If you don't want to do this, leave [UDF Name] blank "
}
}
| 关键字 | 类型 | 说明 | |
|---|---|---|---|
"1.Off Body Points Num" |
整数(int) | 不需要填写，程序自动生成 | |
"2.Off Body Points filepath(*)" |
字符串(string) | 从txt文件中读取监测点坐标 | |
"3.Get Points From UDF" |
字符串(string) | 从UDF中读取监测点坐标 |
可以在log文件中查看相关参数，如下
|| INFO (2023.03.10 12:31:57) || ***********************************************
|| INFO (2023.03.10 12:31:57) || ******** FDC::OffBodyPoints::showInfo() *******
|| INFO (2023.03.10 12:31:57) || ***********************************************
|| INFO (2023.03.10 12:31:57) || OffBody Points Num = 21
|| INFO (2023.03.10 12:31:57) || 1 ( 2.000, 0.000, 0.000)
|| INFO (2023.03.10 12:31:57) || 2 ( 1.800, 0.000, 0.000)
|| INFO (2023.03.10 12:31:57) || 3 ( 1.600, 0.000, 0.000)
|| INFO (2023.03.10 12:31:57) || 4 ( 1.400, 0.000, 0.000)
|| INFO (2023.03.10 12:31:57) || 5 ( 1.200, 0.000, 0.000)
|| INFO (2023.03.10 12:31:57) || 6 ( 1.000, 0.000, 0.000)
|| INFO (2023.03.10 12:31:57) || 7 ( 0.800, 0.000, 0.000)
|| INFO (2023.03.10 12:31:57) || 8 ( 0.600, 0.000, 0.000)
|| INFO (2023.03.10 12:31:57) || 9 ( 0.400, 0.000, 0.000)
|| INFO (2023.03.10 12:31:57) || 10 ( 0.200, 0.000, 0.000)
|| INFO (2023.03.10 12:31:57) || 11 ( 0.000, 0.000, 0.000)
|| INFO (2023.03.10 12:31:57) || 12 (-0.200, 0.000, 0.000)
|| INFO (2023.03.10 12:31:57) || 13 (-0.400, 0.000, 0.000)
|| INFO (2023.03.10 12:31:57) || 14 (-0.600, 0.000, 0.000)
|| INFO (2023.03.10 12:31:57) || 15 (-0.800, 0.000, 0.000)
|| INFO (2023.03.10 12:31:57) || 16 (-1.000, 0.000, 0.000)
|| INFO (2023.03.10 12:31:57) || 17 (-1.200, 0.000, 0.000)
|| INFO (2023.03.10 12:31:57) || 18 (-1.400, 0.000, 0.000)
|| INFO (2023.03.10 12:31:57) || 19 (-1.600, 0.000, 0.000)
|| INFO (2023.03.10 12:31:57) || 20 (-1.800, 0.000, 0.000)
|| INFO (2023.03.10 12:31:57) ||
|| INFO (2023.03.10 12:31:57) || only show the first 20 points
可以通过两种方式定义检测点空间坐标：
从文件中读取
编写如下格式的txt文件
第一行为点数目，后续为空间点坐标(x y z)，以空格作为分隔符。
可以采用工具函数脚本，生成跟为复杂的空间监测点。
从UDF中读取检测点
填写"3.Get Points From UDF"->"1.UDF Name"对应的函数名称，例如：监测点定义(OffBodyPoints)
注意：
- 这种方式的优先级高于"从文件中读取"，
"1.UDF Name"对应的名称为空字符串或者未UDF文件中找到相关函数，则不能由UDF中读取监测点。- 在[log文件: Brief Summary](#Brief Summary Check)中有关键字
Using OffBodyPoints UDF or not，可通过其确定是成功从UDF中读取监测点。
设定浮体信息(Bodies)
本程序支持N浮体频域计算，当前程序最大浮体数目设定为99，以下为模板(default_FD.json)中的"5.Bodies(*)"。
计算时在CMD窗口，或者在log文件中会显示当前算例计算的浮体相关信息，如下：
|| INFO || ***********************************************************
|| INFO || ************ FDC::Solver_FD::showSolverInfo() *************
|| INFO || ***********************************************************
|| INFO || Bodies Info. count = 5
|| INFO || # Body Name Panel Num Body Type
|| INFO || 1 Barge1 836 FixedBody
|| INFO || 2 Barge2 836 FixedBody
|| INFO || 3 Internal_Lid1 200 Internal Lid
|| INFO || 4 Internal_Lid2 200 Internal Lid
|| INFO || 5 Damping Lid 400 Damping Lid(mu=0.000000|Uniform)
浮体类型
首先指定"Body Number(*)"指定浮体数目，目前程序中提供四种浮体模型，分别是
- FloatBody：浮体
- FixedBody：固定浮体
- InternalLid：内部盖(用于消除不规则频率)
- DampingLid：阻尼盖(用于抑制过大的间隙共振响应)
各个类型的浮体会输出不同的结果文件。
| 体类型 | 波浪力 | 附加质量 | 阻尼系数 | 面元文件 | 水线面文件 | 静力学参数文件 |
|---|---|---|---|---|---|---|
| FloatBody | ☑️ | ☑️ | ☑️ | ☑️ | ☑️ | ☑️ |
| FixedBody | ☑️ | ☑️ | ☑️ | ☑️ | ||
| InternalLid | ☑️ | |||||
| DampingLid | ☑️ |
浮体构成
参考程序框架结构，浮体由基础的模型构成，基础模型包括：
- "1.Basic Info(*)"：浮体的基本信息，名称以及类型等
- "2.Calculate Info(*)":
- "3.Panel Model(*)":
- "4.HydroStatic Model":
- "5.Mass Model(*)":
- "5.Damping Lid Model":
以下表格显示了不同类型浮体的不同构成。
| 体类型 | "1.Basic Info(*)" | "2.Calculate Info(*)" | "3.Panel Model(*)" | "4.HydroStatic Model" | "5.Mass Model(*)" | "5.Damping Lid Model" |
|---|---|---|---|---|---|---|
| FloatBody | ☑️ | ☑️ | ☑️ | ☑️ | ☑️ | |
| FixedBody | ☑️ | ☑️ | ☑️ | ☑️ | ||
| InternalLid | ☑️ | ☑️ | ☑️ | ☑️ | ||
| DampingLid | ☑️ | ☑️ | ☑️ | ☑️ | ☑️ |
浮体基础信息
指定浮体名称，与浮体类型
"1.Basic Info(*)": {
"1.Name(*)": "FloatBody",
"2.Type(*)": {
"info": "0=FloatBody;1=FixedBody;2=InternalLid;3=DampingLid",
"value": 0
}
}
注意：
浮体相关结果文件均是以浮体名称为前缀，例如：
BodyA受到的波浪干扰力，则结果文件为：BodyA_EXForce.txt对于多个浮体而言，存在浮体之间的相关干扰作用，例如：
BodyA运动对自身水动力系数的影响，则结果文件为：BodyA-BodyA_AddMass.txt
BodyA运动对BodyB水动力系数的影响，则结果文件为：BodyA-BodyB_AddMass.txt
浮体计算控制参数
指定计算曲面积分时的计算精度，对于四边形面元，使用高斯积分进行计算，Gauss N可以指定为2或者4，详细细节可参考博客；对于三角形面元，使用高斯积分进行计算，Tri Rule建议指定为7、10或者11，详细细节可参考博客；
"2.Calculate Info(*)": {
"0.info": "Gauss N = 2 or 4, Tri rule = 7,10,11",
"1.Gauss N(*)": 4,
"2.Tri Rule(*)": 7
}
浮体面元模型
定义浮体的面元模型，其中，存在(*)号标记项的为必填项目。
"3.Panel Model(*)": {
"1.Panel File Path(*)": "Panel/File/Path/FloatBody.GDF",
"2.Patch Num": 0,
"3.Panel Type": "PanelType::Quadrangle",
"4.Ns": 4,
"5.Mesh(*)": {
"1.Scale(*)": 1.0,
"2.Translate(*)": "[0.000000,0.000000,0.000000]",
"3.Rotate(*)": "[0.000000,0.000000,0.000000]",
"4.Symmetry(*)": {
"info": "0=Null;1=AlongX;2=AlongY;3=AlongXY;4=AlongZ",
"value": 0
}
},
"6.Coord": "[0.000000,0.000000,0.000000]"
}
"1.Panel File Path(*)": 指定面元文件路径，建议使用相对路径，目前支持生成的（*.GDF *.pan *.msh(GAMBIT) ）文件，但是面元文件需要提前进行法向矫正。可以利用BEMRosetta.exe软件进行面元前处理操作。所以推荐使用GDF格式面元文件进行计算，对于GAMBIT划分导出的msh格式面元，可以使用本程序进行格式转换，参考这一节
"5.Mesh(*)"：为节点下可以对面元模型进行各项操作，包括缩放(Scale:以坐标原点为基准)，平移(Translate)，旋转(Rotate：不建议使用，可能有bug)，对称(Symmetry:以坐标原点为基准)，按照顺序执行。
"6.Coord"：此项参数为局部坐标系原点在大地坐标系下的位置，不需要填写。
浮体静力学模型
静力学模型是根据面元模型计算得到的静力学参数（大地坐标系下），由上到下分别为
- 模型空间范围，即三个坐标值的最大值，最小值\(m\)
- 湿表面积\(m^2\)
- 水线面积\(m^2\)
- 体积\(m^3\)
- 浮心坐标\(m\)
- 稳性半径\(m\)
- 水线面形心(几何中心，又称漂心)\(m\)
以下取自[test01](#test01: Wigley I (四边形面元))算例中，程序计算Wigley I模型的静力学参数。
"4.HydroStatic Model": {
"1.Dimensions": {
"info": "From (x,y,z) to (x,y,z)",
"value": "[-1.000000,-0.150000,-0.125122] to [1.000000,0.150000,0.000000]"
},
"2.Swet(m^2)": 0.7502960635822891,
"3.Aw(m^2)": 0.41480239999999985,
"4.Volume(m^3)": 0.04168037133857894,
"5.COB(xb,yb,zb)": "[0.000000,0.000000,-0.053496]",
"6.BM(rxx,ryy,rzz)": "[0.052090,2.073846,0.000000]",
"7.Xf(xf,yf)": "[0.000000,-0.000000]"
}
注意：如果结果存在明显错误，请检查面元法向。
浮体质量模型
只有FloatBody才需要添加质量模型，质量模型包含计算浮体运动所需的相关参数，，由上到下分别为：
- 质量\(kg\)
- 惯性半径\(m\)
- 重心\(m\)
- 粘性阻尼系数\(\mu*\dot X\)
- 人工弹簧固有周期
- 开放的自由度，\(0=fixed,1=free\)
以下取自[test01](#test01: Wigley I (四边形面元))算例中，程序计算Wigley I模型的质量模型。
"5.Mass Model(*)": {
"1.Mass(kg)(*)": 40.9,
"2.Inertia Radius(m)(*)": "[0.100000,0.500000,0.500000]",
"3.COG(xg,yg,zg)(*)": "[0.000000,0.000000,0.000000]",
"4.Viscous_Damping(*)": {
"info": "viscous damping for 6 dof,mu*X_dot",
"value": "[0.000000,0.000000,0.000000,0.000000,0.000000,0.000000]"
},
"5.T of Spring(s)(*)": {
"info": "Period of Spring for (surge,sway,yaw)",
"value": "[100.000000,50.000000,100.000000]"
},
"6.Degree of Freedom(*)": {
"info": "[surge,sway,heave,roll,pitch,yaw];0=fixed;1 for free;",
"value": "[1,1,1,1,1,1]"
}
}
注意：目前程序只支持设定一个浮体。
阻尼盖模型
相关算例：
- [test07: Barge_Barge Damping Lid(使用阻尼盖)](#test07: Barge_Barge Damping Lid(使用阻尼盖))
- [test08: Barge_Barge Damping Lid(使用阻尼盖+Newtonian)](#test08: Barge_Barge Damping Lid(使用阻尼盖+Newtonian))
阻尼盖模型是一种常用的自由面阻尼修正方法，目前被广泛应用于WADAM与Hydrostar中，理论基础详见参考文献。
"5.Damping Lid Model": {
"1.Damping Lid Type(*)": {
"info": "0=Uniform;1=NonUniform;2=Newtonian;3=XB_Chen;4=UDF",
"value": 0
},
"2.Damoing Lid Coeff(*)": 0,
"3.UDF": {
"Func Name In Lua(*)": "",
"info": "this is the UDF Func Name which should be found in Lua Script"
}
}
对于阻尼盖模型而言，最重要的是定义阻尼系数，不同阻尼系数\(\mu_0\)形式如下表所示
| 阻尼盖类型 | 对应数值 | 阻尼形式 | 说明 |
|---|---|---|---|
Uniform |
0 |
\(\frac{\partial \phi}{\partial z}=\frac{(\omega+\text i\mu)^2}{g}\phi\),\(\mu=\mu_0\) | |
NonUniform |
1 |
\(\frac{\partial \phi}{\partial z}=\frac{(\omega+\text i\mu)^2}{g}\phi\),\(\mu=\mu_0f(\omega)g(x)\) | 仅测试 |
Newtonian |
2 |
\(\frac{\partial \phi}{\partial z}=k{(1+\text i\varepsilon)^2}\phi\),\(\varepsilon=\mu_0/\omega\) | WADAM |
XB_Chen |
3 |
\(\frac{\partial \phi}{\partial z}=\frac{\left(\omega^2+\text i\omega\mu\right)}{g}\phi\),\(\mu=\mu_0\) | Hydrostar |
UDF |
4 |
\(\frac{\partial \phi}{\partial z}=\frac{(\omega+\text i\mu)^2}{g}\phi\),\(\mu=UDF(x,y,z,\omega)\) | 自定义 |
注意：
阻尼系数为公式中对应的\(\mu_0\)，UDF的优先级最高，
"Func Name In Lua(*)"指定为空字符串，则表示不使用UDF，以下给出了一个UDF函数实例。在log文件中显示相关的阻尼参数
5 Damping Lid 400 Damping Lid(mu=0.000000|Uniform)
-- UDF Function For Damping Lid Model
function UDF_Damping_Lid_mu(x, y, omega)
-- write code
a = 2.45e-1
b = 2.8e-1
mu = 0.0175 * a * (math.exp(b * omega) - 1)
mu = mu * 2
return mu
end
时域计算配置模板文件(default_TD.json)
时域计算采用间接时域法，在上述频域计算结果上进行展开，计算单浮体、多浮体在不规则波浪下的运动时历，同时考虑系泊系统、连接系统等结构。
时域计算的配置文件包含六个部分。
设定时域算例基本信息
"1.Basic Info(*)"同频域算例基本信息
设定时域计算工况
设定时域计算工况信息，工况包好四个部分，如下图所示，分别是：
- "1.Incident Wave"：入射波浪相关参数
- "2.Wind"：风载荷
- "3.Current"：流载荷
- "4.Other Force"：其他外力
入射波浪参数
入射波浪包含五个部分，如下图所示：
-
入射波浪类型
波浪类型 数值 谱密度函数 说明 RegularWave0输入波高\(H(m)\)，周期\(T(s)\) Jonswap1见参考文献 ITTC2\(S(\omega)=\frac{173H_{1/3}^2}{T_1^4\omega^5}\text e^{-691/(T_1^4\cdot\omega^4)}\) P-M3\(S(\omega)=\frac{8.1\times10^{-3}g^2}{\omega^5}\text e^{-3.11/(H_{1/3}^2\cdot\omega^4)}\) WhiteNoise4见参考文献 Gauss5\(S(f)=\frac{H_s^2}{16}\frac{1}{\sqrt{2\pi\delta^2}}\text{e}^{-\frac{(f-f_p)^2}{2\delta^2}}\) UserWaveEta6自定义波浪序列 目前没用 -
指定浪向
浪向与频域结果相关，为频域结果对应浪向的索引
注意：浪向索引从
0开始，例如：频域结果浪向为0°，45°，90°，则对应索引为0，1，2"2.Wave Direction": { "info": "index of wave direction for FreDomain results (interger)", "value": 0 } -
指定随机数种子
Seed为整数，不同的随机数种子对应着不同的随机相位（对于不规则波浪而言） -
设定聚焦波
指定
"Is Focused Wave"关键字对应数值为true，则设定不规则波浪成分相位，使得其在"[x0,t0](m,s)"时间和位置达到最大值。"4.Focused Wave": { "Is Focused Wave": false, "Position and Time": { "info": "[x0,t0](m,s)", "value": "[0.000000,1000.000000]" } } -
"5.Spectrum Detials"其中包含不同波浪类型对应的参数，自行填写
以ITTC谱为例：以下展示了ITTC谱所需要指定的基本参数。
"3.For ITTC": { "1.Range": { "info": "frequence range", "value": "[0.000000,0.000000]" }, "2.Num": { "info": "Number of frequence", "value": 100 }, "3.Hs": { "info": "Significat Wave Height(m)", "value": 0.0 }, "4.T1": { "info": "Character Period(s)", "value": 0.0 }, "5.Tp": { "info": "Peak Period(s)", "value": 0.0 } }注意：这里的波谱分割按照是按照等间距圆频率分割波谱，如果想要按照等间距波数分割，则需要定义
#define BaseOnK宏，再编译。
风载荷
以下为风载荷作用流程：
"2.Wind": {
"Wind Force Num": 2,
"wind01": {
"1.ForceType(*)": {
"info": "0=NO_Force;1=UDF",
"value": 0
},
"2.OnWhichAnchor": {
"info": "specific the anchor name which the wind force on",
"value": "anchor_name"
},
"3.Force Detials": {
"1.For UDF": {
"info": "this is the UDF Fun ... be found in Lua Script",
"value": "UDF_Function_Name"
}
}
},
"wind02": {
...
-
-
在Lua脚本中定义外力函数
例如：定义了一个名为
UDF_Other_Force(time, dir)，time为时间，单位为s，返回\((F_x,F_y,F_z)\)，大地坐标系下，三个方向外力(可随时间、浪向变化)。-- UDF Function For Other Forces function UDF_Other_Force(time, dir) -- write code -- time(s) dir(int index) Fx = 100.*math.sin(0.1*time) Fy = 0 Fz = 0 return {Fx, Fy, Fz} end -
在
"3.Force Detials"->"1.For UDF"->"value"中填写函数名称。 -
在
"1.ForceType(*)"节点下指定1，使用UDF。 -
在
"2.OnWhichAnchor"节点中指定作用锚点
-
流载荷
其他外力
用于定义额外自己输入的外力
设定时域计算控制参数
控制参数中定义相关的控制参数
"3.Solver Ctrl Info(*)": {
"1.Create Test(*)": true,
"2.Output Interval(*)": {
"info": "Write File Every [ N ] Step",
"value": 500
},
"3.RK4 Solver(*)": {
"1.Time Step(s)": 0.01,
"2.Time Max(s)": 50.0,
"3.Ramping Time(s)": 10.0
},
"4.Lua UDF(*)": {
"info": "the path to lua script",
"value": "UDF.lua"
}
}
"1.Create Test(*)"
指定"1.Create Test(*)" 关键字为true或者false，如果指定为true，则在项目路径下生成测试文件，详见：时域计算测试文件说明
注意：生成测试文件会有点耗时，与频率间隔以及时间步长相关
"2.Output Interval(*)"
指定"2.Output Interval(*)"表示写入数据文件的时间步间隔，例如指定步数为500，则表示每隔500个时间步输出数据至文件。
注意：频繁地打开、写入数据文件至文件比较耗时，本程序内部会设置一个暂存容器，容器大小即为
"2.Output Interval(*)"指定的数值，当容器存满数据时再写入至文件。
"3.RK4 Solver(*)"
时域计算采用4阶数龙格库塔法求解二阶微分方程组，相关理论基础参考博客：数值计算：四阶龙格-库塔法 for 二阶微分方程
"3.RK4 Solver(*)": {
"1.Time Step(s)": 0.01,
"2.Time Max(s)": 50.0,
"3.Ramping Time(s)": 10.0
}
本程序中需要设定3个参数，分别是：
-
"1.Time Step(s)"：时间步大小 -
"2.Time Max(s)"：最大求解时间 -
"3.Ramping Time(s)"：缓冲时间，为了避免产生较大的初始运动幅值，设定缓冲函数为以下形式，设定缓冲时间\(t_r\)，使得在\([0,t_r]\)时间内，作用力缓慢增加，使得作用在物体上的外力缓慢增加，\(\text e^{-5}=0.0067\)\[\begin{align} f_R(t)=1-\text e^{-5\left(\frac{t}{t_r}\right)^3} \end{align} \] 函数图形如下：\(t_r=50s\)
"4.Lua UDF(*)"
与类似[Lua UDF文件定义](#Lua UDF文件定义)。
设定浮体信息
浮体信息参数定义如下：
"4.Bodies(*)": {
"Body Number(*)": 3,
"body01": {
"1.Name": "barge1",
"2.Body Folder Name": "body1"
},
"body02": {
"1.Name": "barge2",
"2.Body Folder Name": "body2"
},
"body03": {
"1.Name": "barge3",
"2.Body Folder Name": "body3"
}
}
时域计算在频域结果基础上展开，因此，需要在"4.Bodies(*)"节点下定义浮体的详细信息：包括
浮体名称
为浮体指定一个名称，后续的锚点定义通过浮体名称进行关联。
参考算例：test01_TD_Wigley_III(白噪声，时域计算)
浮体频域数据文件夹
在"2.Body Folder Name"节点下指定频域结果的文件路径：
注意：控制浮体的自由度，可以在频域结果文件夹下的
*_HydroStatic.json*中修改"6.Degree of Freedom(*)": { "info": "[surge,sway,heave,roll,pitch,yaw];0=fixed;1 for free;", "value": "[1,1,1,1,1,1]" }
设定锚点信息(Anchors)
定义"Anchors Number(*)"数量。
"5.Anchors(*)": {
"Anchors Number(*)": 3,
"anchor01": {
"1.Anchor Name(*)": "an1",
"2.AnchorType": {
"info": "0 = OnEarth(Fixed);1 = OnBody(Floating Body)",
"value": 1
},
"3.Parent(*)": "barge1",
"4.Position(*)": {
"info": "position on body local coord",
"value": "[0.000000,0.000000,0.000000]"
}
},
"anchor02": {
...
}
锚点Anchor属性
"1.Anchor Name(*)"
为当前的锚点指定一个名称，后续均需要通过名称进行关联。
"2.AnchorType"
| AnchorType指定 | 数值 | 坐标系 | 说明 |
|---|---|---|---|
OnEarth(Fixed) |
0 |
大地坐标系 | |
OnBody(Floating Body) |
1 |
局部坐标系 | 表示与浮体坐标系进行关联 |
"3.Parent(*)"
当且仅当指定"2.AnchorType"指定为OnBody时才会有效，指定为浮体名称，浮体名称必须为浮体中对应的名称，参考：设定浮体信息
当且仅当指定"2.AnchorType"指定为OnEarth时，"3.Parent(*)"建议指定为"Earth"虽然没啥用。
"4.Position(*)"
-
如果定义父亲节点为Earth，position为大地坐标系坐标
-
如果定义父亲节点浮体名称，position为局部坐标
局部坐标系\(O'-x'y'z'\)和大地坐标系\(O-xyz\)存在以下关系
物体在\(O'-x'y'z'\)存在相对位移\(\vec X(x,y,z)\) ，转动\(\vec\theta(\alpha,\beta,\gamma)\) 相对位置关系，如下：
\[\vec{OP}=\vec{OO'}+\vec{O'P}\\
\]
设定附加连接结构
"6.Additional Bodies(*)"节点下指定附加的连接结构，包括以下三类：
- 连接缆绳(Cables)
- 防撞垫(Fenders)
- 系泊系统(Moorings)
- 系泊系统第三方库MoorDynLib
面元文件
GDF面元文件
GDF面元格式，GDF格式是WAMIT软件沿用的模型文件格式，只支持四边面元，本程序在其基础上进行修改，使得其支持三角形面元：以下展示了四边形GDF文件开头格式，来源于[test01: Wigley I (四边形面元)](#test01: Wigley I (四边形面元))，在读取.gdf面元文件给出的读取面元的MATLAB程序
4 Type=PanelType::Quadrangle
1.000000 9.806650
0 0
640
-1.000000 0.000000 -0.125122
-1.000000 0.000000 -0.109482
-0.950000 0.003668 -0.109863
-0.950000 0.000000 -0.125122
-1.000000 0.000000 -0.109482
-1.000000 0.000000 -0.093842
-0.950000 0.006951 -0.094604
-0.950000 0.003668 -0.109863
- 第1行：
Ns+空格+字符串，Ns指定为4(四边形面元)或者3(三角形面元)，若是不指定Ns，则默认为4(四边形面元)； - 第2~3行：为WAMIT使用GDF的固定格式，本程序中无实际意义，但不可省略。
- 第4行：为面元数目
- 后续为面元的各个顶点，以四边形为例，则4个点一组组成一个面元，四个点的顺序决定法向方向
以下展示了四边形GDF文件开头格式，来源于[test02: Wigley I (三角形面元)](#test02: Wigley I (三角形面元))
3 Type=PanelType::Triangle
1.000000 9.806650
0 0
1154
-0.496737 -0.118563 0.000000
-0.514314 -0.106908 -0.044663
-0.536197 -0.113018 0.000000
-0.418618 -0.111470 -0.066330
-0.385360 -0.118403 -0.062261
-0.382171 -0.099653 -0.092400
0.536197 -0.113018 0.000000
0.513190 -0.107393 -0.043903
0.496737 -0.118563 0.000000
msh格式文件
不建议使用，可能存在法向不统一的问题，建议使用GDF格式面元文件，转换方式参考：面元文件格式转换。
以下为GAMBIT导出的msh文件格式。
其中，
- 行6：
2B8为十六进制数，表示三维空间点的数目。 - 行21：
2a0为十六进制数，表示面元的数目。
(0 "GAMBIT to Fluent File")
(0 "Dimension:")
(2 3)
(10 (0 1 2B8 1 3))
(10 (1 1 2B8 1 3)(
4.1662500000e-001 3.0050000000e-001 -1.8500000000e-001
-4.1662500000e-001 3.0050000000e-001 -1.8500000000e-001
3.1246875000e-001 3.0050000000e-001 -1.8500000000e-001
2.0831250000e-001 3.0050000000e-001 -1.8500000000e-001
...
...
-4.1662500000e-001 3.1823551231e-001 -1.5140452750e-001
-4.1662500000e-001 3.0677016746e-001 -1.6548572283e-001
-4.1662500000e-001 2.9279395019e-001 -1.7612677537e-001
-4.1662500000e-001 2.7706420302e-001 -1.8275125783e-001
))
(0 "Faces:")
(13(0 1 2a0 0))
(13(2 1 2a0 3 0)(
4 1e 89 7a 1c 0 0
4 89 90 79 7a 0 0
4 90 97 78 79 0 0
4 97 9e 77 78 0 0
...
...
4 49 48 176 16f 0 0
4 48 47 17d 176 0 0
4 47 46 184 17d 0 0
4 46 1d 24 184 0 0
))
(0 "Cells:")
(12 (0 1 0 0))
(0 "Zones:")
(45 (2 wall wall)())
(45 (4 interior default-interior)())
pan格式文件
不建议使用，可能存在法向不统一的问题，建议使用GDF格式面元文件，转换方式参考：面元文件格式转换。
Total numbers of points and panels
900 953
Points and the conectivities
-105.748808 0.000000 14.000000
105.100000 0.000000 14.000000
-89.533834 0.000000 0.000000
96.657691 0.000000 0.000000
...
...
136 147 100 101
147 158 99 100
158 169 98 99
169 180 97 98
UDF文件
UDF(User Defined Function)借助lua脚本实现自定义函数功能，可以在UDF脚本中自定相关函数，可以通过程序生成模板文件
UDF文件中可以自定义的函数类型包括以下5类：
- [阻尼系数分布函数(Damping Lid)](#阻尼系数分布函数(Damping Lid))
- [壁面阻尼系数分布函数(Wall Damping)](#壁面阻尼系数分布函数(Wall Damping))
- [壁面阻尼区域定义(Wall Damping Range)](#壁面阻尼区域定义(Wall Damping Range))
- 监测点定义(OffBodyPoints)
- 时域外力定义(OtherForce)
阻尼系数分布函数(Damping Lid)
function UDF_Damping_Lid_mu(x, y, omega)
-- write code
a = 2.45e-1
b = 2.8e-1
mu = 0.0175 * a * (math.exp(b * omega) - 1)
mu = mu * 2
return mu
end
壁面阻尼系数分布函数(Wall Damping)
用于[壁面阻尼模型](#使用Wall Damping Model)(Wall Damping Model)定义阻尼系数的空间(x,y,z)频域\(\omega\)分布函数。
-- UDF_Function_For Wall Damping Model
A = {0.0014, 0.0020, 0.0032, 0.0042, 0.0054} -- 幅值
omega_res = {5.74, 6.54, 7.42, 8.34, 9.24} -- 共振频率
function UDF_Wall_Damping_mu(x, y, z, omega)
-- write code
-- x(m),y(m),z(m),omega(rad/s)
mu = 0
k = 15
for i = 1, 5 do
mu = mu + A[i] * math.exp(-k * math.pow(omega - omega_res[i], 2))
end
return mu * 2
end
壁面阻尼区域定义(Wall Damping Range)
用于[壁面阻尼模型](#使用Wall Damping Model)(Wall Damping Model)定义阻尼区域范围。
例1：壁面阻尼系数分布函数
UDF_Wall_Damping_mu(x, y, z, omega)
-- UDF_Function_For Wall Damping Model
function UDF_Wall_Damping_mu(x, y, z, omega)
-- write code
-- x(m),y(m),z(m),omega(rad/s)
mu = 0.1
return mu
end
例2：壁面阻尼区域范围定义函数
isInWallDampingRange(x, y, z)
-- 以下给出了判断点是否在区域内的一个示例(三维矩形空间)
-- x,y,z为空间点坐标
-- 在区域内返回true，反之返回false
function isInWallDampingRange(x, y, z)
min = {-100, -0.5, -100} -- 三维坐标最小值
max = {100, 0.5, 100} -- 三维坐标最大值
if (x >= min[1] and x <= max[1] and y >= min[2] and y <= max[2] and z >=
min[3] and z <= max[3]) then
return true
else
return false
end
end
例3：壁面阻尼区域范围定义函数(多个区域取并集)
-- 三维空间1
function isInWallDampingRange_1(x, y, z)
min = {-100, -0.08+0.45, -100} -- 三维坐标最小值
max = {100, 0.08+0.45, 100} -- 三维坐标最大值
if (x >= min[1] and x <= max[1] and y >= min[2] and y <= max[2] and z >=
min[3] and z <= max[3]) then
return true
else
return false
end
end
-- 三维空间2
function isInWallDampingRange_2(x, y, z)
min = {-100, -0.08-0.45, -100} -- 三维坐标最小值
max = {100, 0.08-0.45, 100} -- 三维坐标最大值
if (x >= min[1] and x <= max[1] and y >= min[2] and y <= max[2] and z >=
min[3] and z <= max[3]) then
return true
else
return false
end
end
-- 将三维空间1与三维空间2 并集
function isInWallDampingRange(x, y, z)
return isInWallDampingRange_1(x, y, z) or isInWallDampingRange_2(x, y, z)
end
监测点定义(OffBodyPoints)
用于从UDF中读取OffBodyPoints定义监测点
-- 以下给出了在一条线段上均分得到Num个点的例子
-- 其中start_point为起始点，end_point为末端点，Num>=2为点数目
function getOffBodyPoints()
-- write code
start_point = {-20, 0, 0}
end_point = {20, 0, 0}
Num = 10
delta = {}
for k = 1, 3 do delta[k] = (end_point[k] - start_point[k]) / (Num - 1) end
points = {}
for i = 1, Num do
points[i] = {}
for k = 1, 3 do
points[i][k] = start_point[k] + delta[k] * (i - 1)
end
end
return points
end
时域外力定义(OtherForce)
用于定义外力，外力可以为时间(time(s))与浪向(dir)的函数，
注意：浪向为索引值，参考：入射波浪参数说明
-- UDF Function For Other Forces
function UDF_Other_Force(time, dir)
-- write code
-- time(s) dir(int index)
Fx = 0.1
Fy = 0.2
Fz = 0.3
return {Fx, Fy, Fz}
end
Lua基础的数学运算符
-- ************************************** Lua Usage ***************************************
-- ************************* 以下为编写函数可能需要的运算符与数学函数 *************************
-- 1. lua 运算符
-- + 加法 A + B
-- - 减法 A - B
-- * 乘法 A * B
-- / 除法 B / A
-- % 取余 B % A
-- ^ 乘幂 A^2
-- - 负号 -A
-- // 整除运算符(>=lua5.3) 5//2 输出结果 2
-- 2. lua 关系运算符
-- == 等于，检测两个值是否相等，相等返回 true，否则返回 false
-- ~= 不等于，检测两个值是否相等，不相等返回 true，否则返回 false
-- > 大于，如果左边的值大于右边的值，返回 true，否则返回 false
-- < 小于，如果左边的值大于右边的值，返回 false，否则返回 true
-- >= 大于等于，如果左边的值大于等于右边的值，返回 true，否则返回 false
-- <= 小于等于， 如果左边的值小于等于右边的值，返回 true，否则返回 false
-- 逻辑运算符： 与and 或or 非not
-- 3. lua数学函数
-- pi 圆周率 math.pi 3.1415926535898
-- abs 取绝对值 math.abs(-2012) 2012
-- ceil 向上取整 math.ceil(9.1) 10
-- floor 向下取整 math.floor(9.9) 9
-- max 取参数最大值 math.max(2,4,6,8) 8
-- min 取参数最小值 math.min(2,4,6,8) 2
-- pow 计算x的y次幂 math.pow(2,16) 65536
-- sqrt 开平方 math.sqrt(65536) 256
-- mod 取模 math.mod(65535,2) 1
-- modf 取整数和小数部分 math.modf(20.12) 20 0.12
-- randomseed 设随机数种子 math.randomseed(os.time())
-- random 取随机数 math.random(5,90) 5~90
-- rad 角度转弧度 math.rad(180) 3.1415926535898
-- deg 弧度转角度 math.deg(math.pi) 180
-- exp e的x次方 math.exp(4) 54.598150033144
-- log 计算x的自然对数 math.log(54.598150) 4
-- log10 计算10为底，x的对数 math.log10(1000) 3
-- frexp 将参数拆成x * (2 ^ y)的形式 math.frexp(160) 0.625 8
-- ldexp 计算x * (2 ^ y) math.ldexp(0.625,8) 160
-- sin 正弦 math.sin(math.rad(30)) 0.5
-- cos 余弦 math.cos(math.rad(60)) 0.5
-- tan 反正切 math.tan(math.rad(45)) 1
-- asin 反正弦 math.deg(math.asin(0.5)) 30
-- acos 反余弦 math.deg(math.acos(0.5)) 60
-- atan 反正切 math.deg(math.atan(1)) 45
系泊系统配置文件
本程序采用集中质量法计算系泊系统的动态响应，基于开源的计算程序库Github链接，同时对该程序库进行验证。目前该程序库作为系泊分析模块已经集成至开源风机模拟软件 OpenFAST。
配置文件，配置文件说明参考官方说明文档：
----------------------- LINE TYPES ------------------------------------------
1 NTypes - number of LineTypes
Name Diam MassDen EA BA/-zeta Can Cat Cdn Cdt
(-) (m) (kg/m) (N) (N-s/-) (-) (-) (-) (-)
main 0.09 77.7066 384.243E6 -0.8 1.0 0.0 1.6 0.1
---------------------- CONNECTION PROPERTIES --------------------------------
6 NConnects - number of connections including anchors and fairleads
Node Type X Y Z M V FX FY FZ CdA Ca
1 fixed 853.87 0.0 -320.0 0 0 0 0 0 0 0
2 fixed -426.94 739.47 -320.0 0 0 0 0 0 0 0
3 fixed -426.94 -739.47 -320.0 0 0 0 0 0 0 0
4 vessel 5.2 0.0 -70.0 0 0 0 0 0 0 0
5 vessel -2.6 4.5 -70.0 0 0 0 0 0 0 0
6 vessel -2.6 -4.5 -70.0 0 0 0 0 0 0 0
---------------------- LINE PROPERTIES --------------------------------------
3 NLines - number of line objects
Line LineType UnstrLen NumSegs NodeAnch NodeFair Flags/Outputs
1 main 950 40 1 4 p
2 main 950 40 2 5 p
3 main 950 40 3 6 p
---------------------- SOLVER OPTIONS ---------------------------------------
0.001 dtM - time step to use in mooring integration (s)
3.0e6 kBot - bottom stiffness (Pa/m)
3.0e5 cBot - bottom damping (Pa-s/m)
320 WtrDpth - water depth (m)
1.0 dtIC - time interval for analyzing convergence during IC gen (s)
60.0 TmaxIC - max time for IC gen (s)
4.0 CdScaleIC - factor by which to scale drag coefficients during dynamic relaxation (-)
0.001 threshIC - threshold for IC convergence (-)
------------------------ OUTPUTS --------------------------------------------
FairTen1
FairTen2
FairTen3
AnchTen3
L2N4pX
------------------------- need this line --------------------------------------
输出文件有哪些？
目前，软件输出文件主要包含4类：
log文件
log文件中包含程序计算时输出的必要信息，是检查计算参数、调试以及查看错误信息的必要依据
Brief Summary Check
以下为控制参数的简要总结：
|| INFO (2023.03.14 17:33:06) || ***********This is a Brief Summary***********
|| INFO (2023.03.14 17:33:06) || Using Lua Script or not:0
|| INFO (2023.03.14 17:33:06) || Using Wall Damping Model or not:0
|| INFO (2023.03.14 17:33:06) || Using UDF Range or not:0
|| INFO (2023.03.14 17:33:06) || Using OffBodyPoints UDF or not:0
|| INFO (2023.03.14 17:33:06) || ***********End of the Brief Summary***********
| 关键字 | 描述 | 说明 |
|---|---|---|
Using Lua Script or not |
是否使用Lua脚本 | 1=Yes; 0=No |
Using Wall Damping Model or not |
是否使用壁面阻尼模型 | 1=Yes; 0=No |
Using UDF Range or not |
是否使用壁面阻尼模型中自定义空间范围 | 1=Yes; 0=No |
Using OffBodyPoints UDF or not |
是否使用自定义监测点 | 1=Yes; 0=No |
其他提示信息
详见log文件
频域结果文件
以下结果文件，可以采用相关的MATLAB后处理函数进行数据后处理
BodyName为指定的浮体名称
浮体几何信息
BodyName_HydroStatic.json：浮体的静力学参数
BodyName_m.GDF：浮体面元模型
BodyName_m.txt：面元模型详细信息
数据抬头如下：
640 4 (Num:Patches Number || Type:PanelType::Quadrangle)
#Centroid(x,y,z),#Normal(x,y,z),#ds,#P1(x,y,z),#P2(x,y,z),#P3(x,y,z),#P4(x,y,z)
BodyName_WaterLine.txt：水线面
数据抬头如下：
80 (Water Lines Number)
# Node1(x,y,z) Node2(x,y,z)
波浪力结果
BodyName_EXForce1st.txt：入射力+绕射力
BodyName_FKForce1st.txt：入射力
数据抬头如下：
100 3 (Num_Fre Num_Dir)
Fre Dir real(F1) imag(F1) real(F2) imag(F2) real(F3) imag(F3) real(F4) imag(F4) real(F5) imag(F5) real(F6) imag(F6)
运动结果
BodyName_Motion.txt：六自由度运动
数据抬头如下：
100 3 (Num_Fre Num_Dir)
Fre Dir real(X1) imag(X1) real(X2) imag(X2) real(X3) imag(X3) real(X4) imag(X4) real(X5) imag(X5) real(X6) imag(X6)
水动力系数结果
BodyName1-BodyName2_AddMass.txt：附加质量，无因次化方式见数据抬头
100 4.168037e-02(Num_Fre | V | Added Mass(A/(rho*V)))
Fre A11-A66
BodyName1-BodyName2_Damping.txt：阻尼系数，无因次化方式见数据抬头
100 4.168037e-02(Num_Fre | V | Damping Coef(D/(omega*rho*V)))
Fre D11-D66
注意：BodyName1-BodyName2表示Body1运动对于Body2的水动力系数的影响，此项主要用与计算时延函数
空间监测点
ProjectName_OffBodyPoints.txt：空间点坐标以及对应的位置的速度势
500 3 21 (Num_Fre Num_dir Num_OffBodyPoints) |(x,y,z)|(omega,dir,real(phi),imag(phi),abs(i*phi*omega/g))
2.000 0.000 0.000
1.800 0.000 0.000
1.600 0.000 0.000
1.400 0.000 0.000
1.200 0.000 0.000
1.000 0.000 0.000
0.800 0.000 0.000
0.600 0.000 0.000
0.400 0.000 0.000
0.200 0.000 0.000
0.000 0.000 0.000
-0.200 0.000 0.000
-0.400 0.000 0.000
-0.600 0.000 0.000
-0.800 0.000 0.000
-1.000 0.000 0.000
-1.200 0.000 0.000
-1.400 0.000 0.000
-1.600 0.000 0.000
-1.800 0.000 0.000
-2.000 0.000 0.000
2.000000e-02 0.000000e+00 4.204404e-02 -4.904905e+02 9.999806e-01
2.000000e-02 0.000000e+00 3.893519e-02 -4.904891e+02 9.999779e-01
2.000000e-02 0.000000e+00 3.483940e-02 -4.904867e+02 9.999730e-01
2.000000e-02 0.000000e+00 2.964404e-02 -4.904839e+02 9.999672e-01
2.000000e-02 0.000000e+00 2.508884e-02 -4.904821e+02 9.999636e-01
...
...
时域结果文件
时域计算测试文件
如何使用工具函数
处理面元文件
-i, --inputfile <input "*.gdf">
input panel file(input "*.gdf")
-o, --outputfile <output "*.gdf *.pan">
output panel file(output "*.gdf,*.pan")
--calcHydrostatic calculate hydrostatic info for panel
--correctNormal correct Normal for the panel
--scale <set the scale factor,default is 1.0>
--translate <set the translate vector,default is 0,0,0>
--symmetry <set the symmetry type, 0 = NULL;1 = AlongX;2 = AlongY;4 = AlongZ;3 = AlongXY;default = 0;>
| 命令行参数 | 指定参数 | 说明 |
|---|---|---|
-i or --inputfile |
☑️ | 指定输入文件名称(可以为*.GDF,*.msh,*.pan) |
-o or --outputfile |
☑️ | 指定输出文件名称(可以为*.GDF,*.msh,*.pan)，默认在输入文件名称添加 _Out后缀 |
--calcHydrostatic |
添加此参数，则计算静力学参数并输出同名*.json文件 |
|
--correctNormal |
添加此参数，则进行法向矫正(有Bug，暂时没用) |
|
--scale |
☑️ | 对面元文件进行缩放 |
--translate |
☑️ | 对面元文件进行平移 |
--symmetry |
☑️ | 对面元文件进行对称0 = NULL;1 = AlongX;2 = AlongY;4 = AlongZ;3 = AlongXY; |
面元文件格式转换
将GAMBIT导出的.msh面元文件转换为*.pan或者*.gdf
BEM_ALL.exe -i panel.msh -o panel.gdf
BEM_ALL.exe -i panel.msh -o panel.pan
计算面元文件静力学参数
计算gdf面元文件静力学参数
BEM_ALL.exe -i panel.gdf --calcHydrostatic
cmd窗口显示以下输出
C:\Users\chentianwen\Desktop\example_with_results\ToolFunction>..\BEM_ALL_V2.0.exe -i panel.GDF --calcHydrostatic
|| INFO (2023.03.17 00:16:52) || create Log [ default.log ]
|| INFO (2023.03.17 00:16:52) || Current EXE Version: BEM_ALL_V2.0
|| INFO (2023.03.17 00:16:52) || setPoints as PanelType::Quadrangle; Num = 640
|| INFO (2023.03.17 00:16:52) || PanelModel::update()
|| INFO (2023.03.17 00:16:52) || PanelModel initialize()
|| INFO (2023.03.17 00:16:52) || Output config file [panel.GDF.json] successfully !!!
|| INFO (2023.03.17 00:16:52) || output panel file name is not specific, so we decide to add the [Out] subfix
|| INFO (2023.03.17 00:16:52) || Output [panel.GDF_Out] successfully
|| INFO (2023.03.17 00:16:52) || Finished Operation for [ panel.GDF ]
处理面元文件
%% 缩放面元
BEM_ALL.exe -i panel.gdf -o panel_sacled.gdf --scale 100.0
%% 平移面元
BEM_ALL.exe -i panel.gdf -o panel_trans.gdf --translate 100.0,200,300
%% 对称面元(关于以Y轴为法向的面对称)
BEM_ALL.exe -i panel.gdf -o panel_sym.gdf --symmetry 2
使用第三方库
C++编写dll
参考博客：Lua调用C++生成的DLL库
Fortran编写dll
建议使用Fortran编写静态库，再使用C++做接口，参考上一节C++编写dll
BEM_GUI.exe使用帮助
软件GUI界面
软件基本操作流程：
软件基本操作流程参见：B站视频
特殊操作
利用bat执行多个算例
BEM_GUI一次只能运行一个算例，可以通过此以下方法，批量生成bat文件(windows平台)，运行bat脚本按顺序执行多个算例。
选择指定其他JSON文件，多选需要运行的json配置文件，即可生成如下的run.bat文件
"C:/Users/chentianwen/Desktop/Bem_GUI_exe_Simple/solver/BEM_ALL_V2.0.exe" -n 4 -j C:/Users/chentianwen/Desktop/Bem_GUI_exe_Simple/example_with_results/FDC cases/test01_Wigley.json --log C:/Users/chentianwen/Desktop/Bem_GUI_exe_Simple/example_with_results/FDC cases/test01_wigley/default.log
"C:/Users/chentianwen/Desktop/Bem_GUI_exe_Simple/solver/BEM_ALL_V2.0.exe" -n 4 -j C:/Users/chentianwen/Desktop/Bem_GUI_exe_Simple/example_with_results/FDC cases/test01_Wigley_Lid.json --log C:/Users/chentianwen/Desktop/Bem_GUI_exe_Simple/example_with_results/FDC cases/test01_wigley_lid/default.log
"C:/Users/chentianwen/Desktop/Bem_GUI_exe_Simple/solver/BEM_ALL_V2.0.exe" -n 4 -j C:/Users/chentianwen/Desktop/Bem_GUI_exe_Simple/example_with_results/FDC cases/test02_Wigley_tri.json --log C:/Users/chentianwen/Desktop/Bem_GUI_exe_Simple/example_with_results/FDC cases/test02_wigley_tri/default.log
"C:/Users/chentianwen/Desktop/Bem_GUI_exe_Simple/solver/BEM_ALL_V2.0.exe" -n 4 -j C:/Users/chentianwen/Desktop/Bem_GUI_exe_Simple/example_with_results/FDC cases/test03_Barge.json --log C:/Users/chentianwen/Desktop/Bem_GUI_exe_Simple/example_with_results/FDC cases/test03_barge/default.log
"C:/Users/chentianwen/Desktop/Bem_GUI_exe_Simple/solver/BEM_ALL_V2.0.exe" -n 4 -j C:/Users/chentianwen/Desktop/Bem_GUI_exe_Simple/example_with_results/FDC cases/test04_Barge_Lid.json --log C:/Users/chentianwen/Desktop/Bem_GUI_exe_Simple/example_with_results/FDC cases/test04_barge_lid/default.log
"C:/Users/chentianwen/Desktop/Bem_GUI_exe_Simple/solver/BEM_ALL_V2.0.exe" -n 4 -j C:/Users/chentianwen/Desktop/Bem_GUI_exe_Simple/example_with_results/FDC cases/test05_Barge_Barge.json --log C:/Users/chentianwen/Desktop/Bem_GUI_exe_Simple/example_with_results/FDC cases/test05_barge_barge/default.log
从文件夹中算例数据
菜单栏->后处理->打开数据文件夹
选择文件夹中的带有项目名称的json文件，即可进入绘图选项面板，如下图所示
测试算例
收敛性及并行测试
硬件条件
计算模型与工况
| 频率 | 浪向 | 检测点 |
|---|---|---|
0.1:0.1:10(100) |
0:45:90(3) |
0 |
收敛性与并行测试结果
收敛性测试算例文件夹：ConvergentTest
并行测试算例文件夹：Parallel 、Parallel_OpenBLAS
| 面元数目 | 积分方法 | 并行线程数 | Eigen(s) | OpenBLAS(s) | Eigen耗时比 |
|---|---|---|---|---|---|
| 160 | \(2\times 2\ Gauss\) | 4 | 0.967 | - | 0.003 |
| 320 | \(2\times 2\ Gauss\) | 4 | 3.985 | - | 0.013 |
| 400 | \(2\times 2\ Gauss\) | 4 | 6.204 | - | 0.020 |
| 640 | \(2\times 2\ Gauss\) | 4 | 19.017 | - | 0.061 |
| 800 | \(2\times 2\ Gauss\) | 4 | 29.283 | - | 0.094 |
| 1200 | \(2\times 2\ Gauss\) | 4 | 67.932 | - | 0.219 |
| 2400 | \(2\times 2\ Gauss\) | 1 | 834.589 | 601.239 | 2.969 |
| 2400 | \(2\times 2\ Gauss\) | 2 | 479.91 | 330.901 | 1.547 |
| 2400 | \(2\times 2\ Gauss\) | 4 | 310.262 | 201.323 | 1 |
| 2400 | \(2\times 2\ Gauss\) | 8 | 242.376 | 148.594 | 0.781 |
| 2400 | \(2\times 2\ Gauss\) | 16 | 219.537 | 124.019 | 0.708 |
| 4000 | \(2\times 2\ Gauss\) | 4 | 1035.53 | - | 3.338 |
频域算例
单浮体频域计算(无限水深)
以下为求解在无限水深条件下，单浮体
test01: Wigley I (四边形面元)
test02: Wigley I (三角形面元)
test03: Barge (驳船模型)
使用内部盖模型
添加内部盖模型需要指定"body*"模块中"1.Basic Info(*)"中"2.Type(*)"节点为2(InternalLid)
"body02": {
"1.Basic Info(*)": {
"1.Name(*)": "Internal_Lid",
"2.Type(*)": {
"info": "0=FloatBody;1=FixedBody;2=InternalLid;3=DampingLid",
"value": 2
}
},
...
test04: Barge Lid(单浮体+内部盖)
-
模型 类型 Barge FloatBody Lid Internal Lid -
配置文件：test04_Barge_Lid.json
-
计算内容：针对Barge船型(四角形面元)进行频域计算
-
计算结果：添加内部盖(Internal Lid)能够有效消除不规则频率
多浮体频域计算
test05: Barge-Barge(双浮体+监测点)
test06: Wigley_Box(双浮体+验证)
使用阻尼盖模型
添加阻尼盖模型需要指定"body*"模块中"1.Basic Info(*)"中"2.Type(*)"节点为3(DampingLid)，并且在"5.Damping Lid Model"中选择相应的阻尼盖类型以及定义阻尼系数，可参考使用阻尼盖模型。
test07: Barge_Barge Damping Lid(使用阻尼盖)
-
-
************************************************ ***** FDC::Solver_FD::showSolverInfo() ********* ************************************************ Bodies Info. count = 5 # Body Name Panel Num Body Type 1 Barge1 836 FixedBody 2 Barge2 836 FixedBody 3 Internal_Lid1 200 Internal Lid 4 Internal_Lid2 200 Internal Lid 5 Damping Lid 400 Damping Lid(mu=0.000000|Uniform)
test08: Barge_Barge Damping Lid(使用阻尼盖+Newtonian)
算例设置同[test07](#test07: Barge_Barge Damping Lid(使用阻尼盖))，选择Newtonian作为阻尼系数频率分布，阻尼形式见如何使用阻尼盖模型
-
"5.Damping Lid Model": { "1.Damping Lid Type(*)": { "info": "0=Uniform;1=NonUniform;2=Newtonian;3=XB_Chen;4=UDF", "value": 2 }, "2.Damoing Lid Coeff(*)": 0.003, "3.UDF": { "Func Name In Lua(*)": "", "info": "this is the UDF Func Name which should be found in Lua Script" } }
使用壁面阻尼模型
test09: Barge_Barge Wall_Damping(使用壁面阻尼模型)
-
壁面阻尼模型参数设置如下
"2.Wall Damping Model(*)": { "1.UDF Func Name": "UDF_Wall_Damping_mu", "2.Range": { "info": "[[x_min,y_min,z_min],[x_max,y_max,z_max]]", "value": "[[0.000000,0.000000,0.000000],[0.000000,0.000000,0.000000]]" }, "3.Range From UDF": { "UDF Func Name": "isInWallDampingRange", "info": "define the range in UDF;in range \"return true\",else \"return false\";" } }相关的UDF函数如下：
-- UDF_Function_For Wall Damping Model A = {0.0014, 0.0020, 0.0032, 0.0042, 0.0054} -- 幅值 omega_res = {5.74, 6.54, 7.42, 8.34, 9.24} -- 共振频率 function UDF_Wall_Damping_mu(x, y, z, omega) -- write code -- x(m),y(m),z(m),omega(rad/s) mu = 0 k = 15 for i = 1, 5 do mu = mu + A[i] * math.exp(-k * math.pow(omega - omega_res[i], 2)) end return mu * 2 end -- 以下给出了判断点是否在区域内的一个示例(三维矩形空间) -- x,y,z为空间点坐标 -- 在区域内返回true，反之返回false function isInWallDampingRange(x, y, z) min = {-100, -0.08, -100} -- 三维坐标最小值 max = {100, 0.08, 100} -- 三维坐标最大值 if (x >= min[1] and x <= max[1] and y >= min[2] and y <= max[2] and z >= min[3] and z <= max[3]) then return true else return false end end
特殊算例
以下为其他类型的算例
test10: 方形box (检测浮体附近空间波形)
test11: 半潜平台(面元数目10000+)
test12: N浮体频域水动力计算
test13: 半球(有限水深)
test14: 风机下半平台
test15: 不知道怎么描述
test16: 月池共振(添加阻尼盖)
时域算例
test01_TD_Wigley_III(白噪声，时域计算)
test02_TD_Barge_Barge(双浮体，Cable，Fender，Mooring)
工具函数
BEMRosetta进行面元处理
BEMRosetta是一个很有用的GDF面元可视化处理工具，可以对面元文件进行：
- 平移、旋转、缩放等几何操作
- 显示并统一面元法向
- 划分网格、内部盖网格
- 以及其他NB的功能
该软件在Github上已开源，BEMRosetta仓库链接
MATLAB后处理数据
读取波浪力(ReadEXForce.m)
function out = ReadEXForce(filename)
% 读取波浪力
fp=fopen(filename,'r');
line=fgetl(fp);
str=strsplit(line,{' '});
N_fre=str2double(str{1});
N_dir=str2double(str{2});
line=fgetl(fp);
C=textscan(fp,"%f %f %f %f %f %f %f %f %f %f %f %f %f %f");
fclose(fp);
out.fre=C{1,1}(1:N_dir:end,:);
out.dir=C{1,2}(1:1:N_dir,:);
out.EXForce=cell(N_dir,6);
for k=1:1:6
EXF(:,k)=complex(C{1,2*k+1},C{1,2*k+2});
end
for i=1:1:N_dir
for k=1:1:6
out.EXForce{i,k}=EXF(i:N_dir:end,k);
end
end
end
读取水动力系数(ReadHydroCoef.m)
读取附加质量、阻尼系数结果
function out = ReadHydroCoef(filename)
% 读取水动力系数
fp=fopen(filename,'r');
C=textscan(fp,"%f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f",'headerlines',2);
fclose(fp);
out.fre=C{1,1};
out.AD=cell(6,6);
for i=1:1:6
for j=1:1:6
out.AD{i,j}=C{1,(i-1)*6+j+1};
end
end
end
读取监测点结果(ReadOffBodyPoints.m)
function out= ReadOffBodyPoints(filename)
%读取监测点数据
fp=fopen(filename,'r');
line=fgetl(fp);
str=strsplit(line,{' '});
N_fre=str2double(str{1});
N_dir=str2double(str{2});
N_op =str2double(str{3});
C=textscan(fp,"%f %f %f",N_op);
C1=textscan(fp,"%f %f %f %f %f");
fclose(fp);
out.Points=[C{1,1},C{1,2},C{1,3}];
out.Eta=cell(1,N_dir);
out.Phi=cell(1,N_dir);
Fre=C1{1,1};
Dir=C1{1,2};
Phi=complex(C1{1,3},C1{1,4});
Eta=C1{1,5};
out.Fre=Fre(1:N_dir*N_op:end,:);
out.Dir=sort(Dir(1:N_fre*N_op:end,:));
Phi=reshape(Phi,N_op,[]);
Eta=reshape(Eta,N_op,[]);
for i=1:1:N_dir
out.Phi{1,i}=Phi(:,i:N_dir:end).';
out.Eta{1,i}=Eta(:,i:N_dir:end).';
end
for i=1:1:N_dir
out.eta{1,i}=complex(0,1)/9.81*out.Phi{1,i}.*repmat(out.Fre,1,N_op);
end
end
读取面元文件
读取.txt面元文件
function TxT = ReadPanelTxT(filename)
% 读取面元txt文件
tic
fp=fopen(filename,'r');
line=fgetl(fp);
str=strsplit(line,{' '});
Num=str2double(str{1});
Type=str2double(str{2});
line=fgetl(fp);
if Type==4
C=textscan(fp,'%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f');
else
C=textscan(fp,'%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f');
end
fclose(fp);
TxT.Centroid=[C{1,1},C{1,2},C{1,3}];
TxT.Normal=[C{1,4},C{1,5},C{1,6}];
TxT.ds=[C{1,7}];
TxT.P1=[C{1,8},C{1,9},C{1,10}];
TxT.P2=[C{1,11},C{1,12},C{1,13}];
TxT.P3=[C{1,14},C{1,15},C{1,16}];
if Type==4
TxT.P4=[C{1,17},C{1,18},C{1,19}];
for i=1:1:Num
TxT.Patches(4*i-3:4*i,:)=[TxT.P1(i,:);TxT.P2(i,:);TxT.P3(i,:);TxT.P4(i,:)];
end
else
for i=1:1:Num
TxT.Patches(3*i-2:3*i,:)=[TxT.P1(i,:);TxT.P2(i,:);TxT.P3(i,:)];
end
end
TxT.Type=Type;
TxT.Num=Num;
toc
end
读取.gdf面元文件
function [out] = ReadGDF(filename)
%READGDF 读取GDF文件,并计算相关参数
fp=fopen(filename,'rt');
C=textscan(fp,'%f %f %f ','headerlines',4);
out.x=C{1,1};out.y=C{1,2};out.z=C{1,3};
fclose(fp);
out.Num=length(out.x)/4;
for i=1:1:out.Num*4
index=ceil(i/4)*4;
Point1=[out.x(index-3),out.y(index-3),out.z(index-3)];
Point2=[out.x(index-2),out.y(index-2),out.z(index-2)];
Point3=[out.x(index-1),out.y(index-1),out.z(index-1)];
Point4=[out.x(index),out.y(index),out.z(index)];
out.Patches(index-3:index,:)=[Point1;Point2;Point3;Point4];
end
out.ds=zeros(out.Num,1);
out.Centroid=zeros(out.Num,3);
out.Normal=zeros(out.Num,3);
for i=1:1:out.Num
P1=out.Patches(4*i-3,:);
P2=out.Patches(4*i-2,:);
P3=out.Patches(4*i-1,:);
P4=out.Patches(4*i-0,:);
n1=P1-P3;n2=P2-P4;
sintheta=sqrt(1-(n1*n2'/(norm(n1)*norm(n2)))^2);
out.ds(i)=0.5*(norm(n1)*norm(n2)*sintheta);
out.Centroid(i,:)=(P1+P2+P3+P4)/4;
out.Normal(i,:)=cross(n1,n2)/(norm(cross(n1,n2)));
end
end
读取实数，复数矩阵文件
读取实数(.txt)文件(ReadMat.m)
function mat = ReadMat(filename,M)
tic
fp=fopen(filename,'r');
line=fgetl(fp);
str=strsplit(line,{' '});
C=textscan(fp,'%f,%f');
fclose(fp);
x=reshape(C{1,1}.',M,[]).';
fx=reshape(C{1,2}.',M,[]).';
toc
mat.x=x;
mat.fx=fx;
end
读取复数(.txt)文件(ReadMatC.m)
function mat = ReadMatC(filename)
tic
fp=fopen(filename,'r');
line=fgetl(fp);
str=strsplit(line,{' '});
M=str2double(str{1});
N=str2double(str{2});
C=textscan(fp,'%f,%f');
fclose(fp);
mat=complex(reshape(C{1,1},N,M),reshape(C{1,2},N,M));
mat=mat.';
toc
end
读取系泊系统文件
读取系泊缆绳位置信息(readLineOut.m)
function out = readLineOut(filename)
% 读取系泊out文件
fp=fopen(filename);
count=1;
line=fgetl(fp); % 从文件读行
label=strsplit(line,'\t');
label(end)=[];
Nnode=(length(label)-1)/3;
dim=[];
for i=1:1:length(label)
dim=[dim,'%f '];
end
C=textscan(fp,dim,'headerlines',2);
fclose(fp);
out.t=C{1,1};
Nt=length(out.t);
out.Node=cell(Nt,1);
tmp=zeros(Nnode,3);
for tt=1:1:Nt
for k=1:1:Nnode
tmp(k,:)=[C{1,3*k-1}(tt),C{1,3*k}(tt),C{1,3*k+1}(tt)];
end
out.Node{tt,1}=tmp;
end
fprintf("readlineout [ %s ]\n",filename);
end
读取系泊缆绳受力信息(readLineOutForce.m)
function out = readLineOutForce(filename)
% 读取out文件
fp=fopen(filename);
line=fgetl(fp); % 从文件读行
label=strsplit(line,'\t');
label(end)=[];
NN=length(label);
dim=[];
for i=1:1:length(label)
dim=[dim,'%f '];
end
C=textscan(fp,dim,'headerlines',1);
fclose(fp);
out.t=C{1,1};
out.fx=[];
for i=2:1:NN
out.fx(:,i-1)=C{1,i };
end
fprintf("readlineoutforce [ %s ]\n",filename);
end
时历数据统计分析
计算有义值(Cal_Hs.m)
function Hs = Cal_Hs(Eta)
%CAL_HS 计算有义波高
N=length(Eta);
Pset=UpCrossZeroPoints(Eta);
N_P=length(Pset);
H=[];
for k=1:1:N_P-1
eta_tmp=Eta(Pset(k):Pset(k+1));
H=[H;max(eta_tmp)-min(eta_tmp)];
end
%降序排序
H=sort(H,'descend');
%最大三分之一平均值
Num=floor(length(H)/3);
Hs=mean(H(1:Num));
function Pset=UpCrossZeroPoints(Eta)
%计算上跨零点
Pset=[];
for i=1:1:N-1
if(Eta(i)<=0&&Eta(i+1)>=0)
Pset=[Pset;i];
end
end
end
end
谱分析
function [omega,res] = getWelch(data,dt,N)
% 功能：FFT加窗处理
% 输入：data 数据段
% dt 时间间隔
% N 每段处理数据长度N
fs=1/dt;
% 每段数据之间的重叠区域宽度为：N*1/2
% [p,f] = pwelch(data,ones(N,1),N/2,N,fs);
[p,f] = pwelch(data,hamming(N),round(N*0.8),N,fs);
omega = f*2*pi;
res = p/2/pi;
end
function out = myFFT(time,data)
% 功能：FFT
% 输入：time：时间
% data：数据
% 输出：out.Fre:频率 rad
% out.Ampti:幅值
% out.Phase:相位
N=length(data);
fft_Data=fft(data);
Ampti = abs(fft_Data) / (N)*2;
Phase = atan2(imag(fft_Data),real(fft_Data));
dt=time(2)-time(1);
Fre=(0:1:(length(time)-1))/(N*dt)*2*pi;
out.Fre=Fre;
out.Ampti=Ampti;
out.Phase=Phase;
end
生成自由面监测点
clc;clear;
folderpath='FolderWaterLines';
WL1=getWaterLine([folderpath,'\','Barge1_WaterLine.txt']);
WL2=getWaterLine([folderpath,'\','Barge2_WaterLine.txt']);
WL3=getWaterLine([folderpath,'\','Barge3_WaterLine.txt']);
%%
WaterLines={WL1,WL2,WL3};
Xrange=[-3:0.2:-2,-2:0.1:2,2:0.2:3];
Yrange=[-5:0.2:-2,-2:0.1:2,2:0.2:5];
Xrange=unique(Xrange);
Yrange=unique(Yrange);
outfilename='Three_barges';
ofPoints = genFsPoints(Xrange,Yrange,WaterLines,outfilename);
%% 相关函数
function [out] = getWaterLine(filename)
fp=fopen(filename,'rt');
line=fgetl(fp);
str=strsplit(line,{' '});
N=str2double(str{1});
line=fgetl(fp);
C=textscan(fp,"%d,%f,%f,%f,%f,%f,%f");
out.Node1=[C{1,2},C{1,3},C{1,4}];
out.Node2=[C{1,5},C{1,6},C{1,7}];
out.Vec=out.Node2-out.Node1;
fclose(fp);
end
function ofPoints = genFsPoints(Xrange,Yrange,WaterLines,outfilename)
%GENFSPOINTS 输入WaterLine文件，划分自由面网格空间点
% Xrange：X轴范围
% Yrange：Y轴范围
% WaterLines：水线面文件：注：单一的WL文件为封闭的单连通区域
% outfilename：输出至文件
if(isempty(WaterLines) )
warnning('WaterLines is empty!!!');
return;
end
% 划分网格
[X,Y]=meshgrid(Xrange,Yrange);
Grid.X=X';
Grid.Y=Y';
for i=1:1:length(Xrange)
for j=1:1:length(Yrange)
pos=[Grid.X(i,j),Grid.Y(i,j),0];
for k=1:1:length(WaterLines)
if(IsInside(WaterLines{k},pos))
Grid.X(i,j)=nan;
Grid.Y(i,j)=nan;
end
end
end
end
fprintf('划分网格成功！\n');
fprintf('显示绘图！\n');
figure(996)
plot(Grid.X,Grid.Y,'ro');hold on
xlabel("X");ylabel("Y");axis equal;
% grid
ofPoints=[];
for i=1:1:length(Xrange)
for j=1:1:length(Yrange)
if(~isnan(Grid.X(i,j)))
ofPoints=[ofPoints;[Grid.X(i,j),Grid.Y(i,j),0]];
end
end
end
eval(['save ',outfilename,num2str(size(ofPoints,1),'_%d'),' Grid']);
fprintf('保存mat网格数据！\n');
genPointFile([outfilename,num2str(size(ofPoints,1),'_%d.txt')],ofPoints);
fprintf('保存txt网格数据！\n');
end
function re=IsInside(WL,Pos)
% 判断点Pos是否在封闭曲线WL内
N=size(WL.Node1,1);
re_sum=zeros(1,3);
for i=1:1:N
N1=WL.Node1(i,:);
N2=WL.Node2(i,:);
center=(N1+N2)/2;
r=center-Pos;
vec=N2-N1;
tmp=cross(vec,r);
tmp=tmp/norm(tmp);
N1p=N1-Pos;
N2p=N2-Pos;
cosTheta=dot(N1p,N2p)/(norm(N1p)*norm(N2p));
theta=acos(cosTheta);
re_sum=re_sum+tmp*theta;
end
if norm(re_sum)<1.e-3
re=false;
else
re=true;
end
% 位于边界
if isnan(norm(re_sum))
re=true;
end
end
function genPointFile(filename,points)
% 导出至文件
N=size(points,1);
fp=fopen(filename,'w');
fprintf(fp,'%d\n',N);
fprintf(fp,'%10f %10f %10f\n',points');
fclose(fp);
end
参考文献
[1] 陈天文. 旁靠外输系统水动力性能数值模拟研究[D].上海交通大学,2023.
