LongLong's Blog
分享IT技术，分享生活感悟，热爱摄影，热爱航天。
Android中使用蓝牙通信
蓝牙(Bluetooth)是一种无线数据和语音通信开放的全球规范，基于低成本的近距离无线连接，为固定和移动设备建立通信环境的一种特殊的近距离无线技术连接。相比WIFI，蓝牙的优势在于连接不需要借助路由器或AP设备，可通过简单配对过程实现设备之间的快速连接。
以下将通过将Android系统作为服务端实现一个基于蓝牙的Socket服务器，客户端通过蓝牙连接后与其实现Socket通信。
1. Android的Service
为了让蓝牙服务在启动后一直能够保持运行状态，需要使用Android中提供的Service，其能够在后台也保持运行状态。如实现的BluetoothService需要继承自Service类，并覆盖onBind和onStartCommand两个方法来实现服务的逻辑。
public class BluetoothService extends Service { public IBinder onBind(Intent intent) { return null; } public int onStartCommand(Intent intent, int flags, int startId) { return START_STICKY; } }在Activity中可以启动Service，启动的代码如下
startService(new Intent(getBaseContext(), BluetoothService.class));
2. 启动蓝牙
首先需要在App中设置访问设备蓝牙的权限，在AndroidManifest.xml中增加以下权限申请
<uses-permission android:name="android.permission.BLUETOOTH"/> <uses-permission android:name="android.permission.BLUETOOTH_ADMIN"/>使用BluetoothAdapter类中的getDefaultAdapter方法可以获得系统的蓝牙适配器，如果系统没有蓝牙设备将会返回null，需要对此进行判断防止发生空指针异常。正常获得系统的蓝牙适配器后调用enable方法开启系统的蓝牙功能
BluetoothAdapter adapter = BluetoothAdapter.getDefaultAdapter(); if (adapter == null) { return -1; } //这里可以判断是否已经开启，不过重复开启也不会导致程序出现错误 adapter.enable();为了让其他设备能够发现此蓝牙服务，需要请求开启蓝牙的可发现性，使用申请提示待用户需要点击确认后，可启动蓝牙的可发现性一段时间（可设置，最大300秒）
Intent intent = new Intent(BluetoothAdapter.ACTION_REQUEST_DISCOVERABLE); //早期版本的Android，支持设置0时长期保持可发现状态 intent.putExtra(BluetoothAdapter.EXTRA_DISCOVERABLE_DURATION, 300); startActivity(intent);
3. 发现其他蓝牙设备
调用蓝牙适配器的startDiscovery方法可以搜索附近的蓝牙设备，但从Android 6开始启动蓝牙搜索需要获取定位的权限，同时定位的权限是不支持通过AndroidManifest中进行预先申请的，需要使用申请提示待用户需要点击确认，在Activity中增加以下代码
ActivityCompat.requestPermissions(this, new String[]{Manifest.permission.ACCESS_COARSE_LOCATION, Manifest.permission.ACCESS_FINE_LOCATION, Manifest.permission.ACCESS_BACKGROUND_LOCATION}, 1);startDiscovery方法为异步运行，若想要在发现设备后执行相应的操作，需要设置回调类
//过滤器，仅关注发现事件 IntentFilter filter = new IntentFilter(); filter.addAction(BluetoothDevice.ACTION_FOUND); //注册回调类 registerReceiver(new BroadcastReceiver() { //实现回调函数，在发现设备后输出设备的名称或地址 @Override public void onReceive(Context context, Intent intent) { BluetoothDevice device = intent.getParcelableExtra(BluetoothDevice.EXTRA_DEVICE); if (device != null) { //当没有设置设备名称是getName方法将会返回null String name = device.getName(); Log.i("bluetooth", name != null ? name : device.toString()); } } }, filter);
4. 实现蓝牙的SocketServer
实现与普通的TCP SocketServer非常类似，首先通过listen方法得到一个Socket，然后循环调用accept方法等待连接，当连接成功后获得输入输出流实现接收和发送数据，并通过一个线程来启动Server
new Thread(new Runnable() { @Override public void run() { try { UUID uuid = UUID.fromString("b383ee98-fd44-4972-8e3e-48bf542b0a7b"); //设置名称和uid来供客户端连接时进行识别 BluetoothServerSocket socket = adapter.listenUsingRfcommWithServiceRecord("Rfc Server", uuid); Log.i("bluetooth", uuid.toString()); Log.i("bluetooth", adapter.getName()); Log.i("bluetooth", adapter.getAddress()); while (true) { BluetoothSocket client = socket.accept(); Log.i("bluetooth", "Accept"); if (client != null) { OutputStream os = client.getOutputStream(); PrintStream out = new PrintStream(os); out.println("Hello!"); } } } catch (IOException e) { e.printStackTrace(); } } }).start();
5. 通过RPC与蓝牙Service通信
可以在蓝牙Service中再启动一个HTTP Server，然后通过HTTP进行RPC，同时在接到RPC调用后需要使用建立的蓝牙连接则需要对蓝牙Server和HTTP Server之间对连接操作进行同步。大致方法就是将建立的蓝牙连接保存在一个Set容器中，将此容器做为两个线程的共享对象，并以此进行操作连接的同步。具体实现如下
//实例化共享的连接集合 final Set<BluetoothSocket> clients = new HashSet<BluetoothSocket>(); //在蓝牙服务线程中 while (true) { BluetoothSocket client = socket.accept(); Log.i("bluetooth", "Accept"); if (client != null) { //进行同步 synchronized (clients) { clients.add(client); OutputStream os = client.getOutputStream(); PrintStream out = new PrintStream(os); out.println("Hello from Bluetooth!"); } } } //在HTTP服务中 AsyncHttpServer server = new AsyncHttpServer(); server.get("/", new HttpServerRequestCallback() { @Override public void onRequest(AsyncHttpServerRequest request, AsyncHttpServerResponse response) { response.send("Hello!!!"); //进行同步 synchronized (clients) { //对每个连接发送消息 for (BluetoothSocket c : clients) { try { OutputStream os = c.getOutputStream(); PrintStream out = new PrintStream(os); out.println("Hello from Http!"); } catch (IOException e) { e.printStackTrace(); } } } } });
6. 实现蓝牙的Client
这里使用Python实现了一个蓝牙的Client，使用pybluez库作为操作蓝牙的基础库，连接完成后接收数据，并关闭连接
import bluetooth; sock = bluetooth.BluetoothSocket(bluetooth.RFCOMM); host = 'D4:12:43:66:C4:A3'; port = 3; sock.connect((host, port)); print(sock.recv(1024)); print(sock.recv(1024)); sock.close();Client将会从Server收到一个Hello from Bluetooth!的字符串，并输出，同时通过HTTP调用HTTP Server，会再输出一个Hello from Http!
7. 总结
总体感觉Android的坑还是挺多的，各个版本之间的差异还是比较明显的，同时在相关的文档中又没有很清楚的说明，如果需要适配多个版本的Android系统的额外工作量还是比较大的，同时需要的各种访问权限需要在UI中手动确认，作为服务程序每次启动和重启都要手动点击确认确实还是不太方便。
有限元方法在计算流体力学中的应用
在传统的CFD中广泛使用的方法一般为有限容积法(FVM)，基于的是控制容积内积分形式的守恒方程，对整个求解区域进行离散得到整个求解域的方程组。而有限元法则是从变分原理（或加权残值）出发，求使泛函在整个求解域内取极值的节点物理量。有限元法的优点为对于复杂边界有很好的适应性，缺点为对于流体力学的方程中无相应的变分原理，而是直接使用加权残值法处理，求解过程没有明显的物理含义。
1. 流体力学的基本方程组
对于一般不可压缩流体，主要的控制方程为流体的连续性方程和NS（Navier-Stokes ）方程
连续性方程
\[ \frac{\partial \rho}{\partial t} + \nabla \rho \boldsymbol{v} = 0\] 动量方程（Navier-Stokes 方程）
\[ \frac{\partial \rho\boldsymbol{v}}{\partial t} + \nabla \cdot \rho\boldsymbol{v} \boldsymbol{v} = -\nabla p + \mu \nabla^2 \boldsymbol{v} + \boldsymbol{f}\] 方程中一共有速度的3个分量和压力4个未知量，并且有4个方程，可以完成求解。
如果流体中涉及能量传递，则需要补充一个能量方程进行温度场的求解；如果流体是可压缩的，则需要补充一个状态方程完成密度场的求解。
2. 流体力学方程组的有限元格式
考虑二维不可压缩流体，使用标准的形函数对待求解的节点速度和节点压力进行插值 \[ \{\boldsymbol{u}\} = [\boldsymbol{N}]\{\boldsymbol{u^e}\} \] \[ \{\boldsymbol{v}\} = [\boldsymbol{N}]\{\boldsymbol{v^e}\} \] \[ \{p\} = [\boldsymbol{N}]\{p^e\} \] 带入到连续性方程和动量方程 \[ \frac{\partial[\boldsymbol{N}]}{\partial x}\{\boldsymbol{u^e}\} + \frac{\partial[\boldsymbol{N}]}{\partial y}\{\boldsymbol{v^e}\} = \{0\}\] \[ [\boldsymbol{N}]\frac{\partial \{\boldsymbol{u^e}\}}{\partial t} + [\boldsymbol{N}]\{\boldsymbol{u^e}\}\frac{\partial[\boldsymbol{N}]}{\partial x}\{\boldsymbol{u^e}\} + [\boldsymbol{N}]\{\boldsymbol{v^e}\}\frac{\partial[\boldsymbol{N}]}{\partial y}\{\boldsymbol{u^e}\} + \frac{1}{\rho}\frac{\partial[\boldsymbol{N}]}{\partial x}\{p^e\} \] \[- \frac{\mu}{\rho}(\frac{\partial^2[\boldsymbol{N}]}{\partial x^2} + \frac{\partial^2[\boldsymbol{N}]}{\partial y^2})\{\boldsymbol{u^e}\} = \{0\}\] \[ [\boldsymbol{N}]\frac{\partial \{\boldsymbol{v^e}\}}{\partial t} + [\boldsymbol{N}]\{\boldsymbol{u^e}\}\frac{\partial[\boldsymbol{N}]}{\partial x}\{\boldsymbol{v^e}\} + [\boldsymbol{N}]\{\boldsymbol{v^e}\}\frac{\partial[\boldsymbol{N}]}{\partial y}\{\boldsymbol{v^e}\} + \frac{1}{\rho}\frac{\partial[\boldsymbol{N}]}{\partial y}\{p^e\} \] \[- \frac{\mu}{\rho}(\frac{\partial^2[\boldsymbol{N}]}{\partial x^2} + \frac{\partial^2[\boldsymbol{N}]}{\partial y^2}) \{\boldsymbol{v^e}\} = \{0\}\] 对以上方程进行加权积分，并进行分布积分，最终得到 \[\iint[\boldsymbol{N}]^T(\frac{\partial[\boldsymbol{N}]}{\partial x}\{\boldsymbol{u^e}\} + \frac{\partial[\boldsymbol{N}]}{\partial y}\{\boldsymbol{v^e}\})dxdy=\{0\}\] \[\iint[\boldsymbol{N}]^T[\boldsymbol{N}]dxdy\frac{\{\boldsymbol{u^e}\}}{\partial t} + \iint[\boldsymbol{N}]^T[\boldsymbol{N}]\{\boldsymbol{u^e}\}\frac{\partial[\boldsymbol{N}]}{\partial x}dxdy\{\boldsymbol{u^e}\}\] \[+\iint[\boldsymbol{N}]^T[\boldsymbol{N}]\{\boldsymbol{v^e}\}\frac{\partial[\boldsymbol{N}]}{\partial y}dxdy\{\boldsymbol{u^e}\} + \frac{1}{\rho}\iint[\boldsymbol{N}]^T\frac{\partial[\boldsymbol{N}]}{\partial x}dxdy\{p\}\] \[+\frac{\mu}{\rho}(\iint\frac{\partial[\boldsymbol{N}]^T}{\partial x}\frac{\partial[\boldsymbol{N}]}{\partial x}dxdy + \iint\frac{\partial[\boldsymbol{N}]^T}{\partial y}\frac{\partial[\boldsymbol{N}]}{\partial y}dxdy)\{\boldsymbol{u^e}\}=\{0\}\] \[\iint[\boldsymbol{N}]^T[\boldsymbol{N}]dxdy\frac{\{\boldsymbol{v^e}\}}{\partial t} + \iint[\boldsymbol{N}]^T[\boldsymbol{N}]\{\boldsymbol{u^e}\}\frac{\partial[\boldsymbol{N}]}{\partial x}dxdy\{\boldsymbol{v^e}\}\] \[+\iint[\boldsymbol{N}]^T[\boldsymbol{N}]\{\boldsymbol{v^e}\}\frac{\partial[\boldsymbol{N}]}{\partial y}dxdy\{\boldsymbol{v^e}\} + \frac{1}{\rho}\iint[\boldsymbol{N}]^T\frac{\partial[\boldsymbol{N}]}{\partial x}dxdy\{p\}\] \[+\frac{\mu}{\rho}(\iint\frac{\partial[\boldsymbol{N}]^T}{\partial x}\frac{\partial[\boldsymbol{N}]}{\partial x}dxdy + \iint\frac{\partial[\boldsymbol{N}]^T}{\partial y}\frac{\partial[\boldsymbol{N}]}{\partial y}dxdy)\{\boldsymbol{v^e}\}=\{0\}\] 从而得到最终的方程组 \[\begin{bmatrix} \boldsymbol{m_{11}} & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & \boldsymbol{m_{33}} \end{bmatrix} \frac{\partial}{\partial t} \begin{bmatrix} \boldsymbol{u^e} \\ p \\ \boldsymbol{v^e}\end{bmatrix} + \begin{bmatrix} \boldsymbol{k_{11}} & \boldsymbol{k_{12}} & 0 \\ \boldsymbol{k_{21}} & 0 & \boldsymbol{k_{23}} \\ 0 & \boldsymbol{k_{32}} & \boldsymbol{k_{33}} \end{bmatrix} \begin{bmatrix} \boldsymbol{u^e} \\ p \\ \boldsymbol{v^e}\end{bmatrix} = \{0\}\] 将其记为 \[[\boldsymbol{M}]\{\dot{\boldsymbol{\Phi}}\}+[\boldsymbol{K}]\{\boldsymbol{\Phi}\}=\{0\}\] 对时间采用中心差分格式 \[[\boldsymbol{M}]\frac{\boldsymbol{\Phi_{2}}-\boldsymbol{\Phi_{1}}}{\Delta t} + [\boldsymbol{K}]\frac{\boldsymbol{\Phi_{1} + \Phi_{2}}}{2}=\{0\}\] 得到最终的方程组 \[([\boldsymbol{M}]+\frac{[\boldsymbol{K}]}{2}\Delta t)\{\boldsymbol{\Phi_{2}}\}=([\boldsymbol{M}]-\frac{[\boldsymbol{K}]}{2}\Delta t)\{\boldsymbol{\Phi_{1}}\}\] 即每进行一个时间步长求解需要求解一次以上的方程组。
3. 有限元程序
有限元程序分别为单元部分和求解器，其中单元部分与弹性力学类似只是更换了形函数的计算过程，求解器为标准的有限元求解器只是更改了自由度编号的规则（弹性力学为按照节点排列自由度，这里把相同的自由度排在一起），编程语言使用Julia。
单元部分的核心代码为生成单元刚度矩阵，利用OO特性相关的逻辑并不依赖于单元的具体形式，可使用于任何二维单元，其代码如下
function getStiffMatrix(self::Element2D) ndof = getDofNum(self); nnode = getNodeNum(self); Ke = zeros(Float64, ndof, ndof); nu = self.property.nu; rho = self.property.rho; #对于每一个高斯点进行数值积分 for p in self.points #形函数导数矩阵(局部坐标系) Der = getDerMatrix(self, p); N = getShapeMatrix(self, p); Nx = zeros(1, nnode); Ny = zeros(1, nnode); u0 = 0.0; v0 = 0.0; #计算高斯点上的速度 for i = 1 : nnode node = self.nodes[i]; u0 += node.vals[U::Dof] * N[1, i]; v0 += node.vals[V::Dof] * N[1, i]; i += 1; end #计算形函数矩阵的导数 for i = 1 : nnode Nx[1, i] = Der[1, i]; Ny[1, i] = Der[2, i]; end J = getJacobi(self, p); k1 = 1; k2 = nnode + 1; k3 = nnode * 2 + 1; k4 = nnode * 3 + 1; #积分变换系数和权值 w = abs(det(J)) * last(p); Ke[k1 : k2 - 1, k1 : k2 - 1] += (N' * u0 * Nx + N' * v0 * Ny + nu * Nx' * Nx + nu * Ny' * Ny) * w; Ke[k1 : k2 - 1, k2 : k3 - 1] += (N' * Nx / rho) * w; Ke[k2 : k3 - 1, k1 : k2 - 1] += (N' * Nx) * w; Ke[k2 : k3 - 1, k3 : k4 - 1] += (N' * Ny) * w; Ke[k3 : k4 - 1, k2 : k3 - 1] += (N' * Ny / rho) * w; Ke[k3 : k4 - 1, k3 : k4 - 1] += (N' * u0 * Nx + N' * v0 * Ny + nu * Nx' * Nx + nu * Ny' * Ny) * w; end return Ke; end求解器部分的核心代码为进行一个时间步长的求解，与标准的有限元求解过程基本一致
function solve(self::Model, dt::Float64) ndof = getDofNum(self); nnode = getNodeNum(self); K = zeros(Float64, ndof, ndof); M = zeros(Float64, ndof, ndof); R = zeros(Float64, ndof, 1); inode::Int32 = 1; mnode::Dict{Node, Int32} = Dict{Node, Int32}(); #建立节点到节点编号的映射表 for node in self.nodes mnode[node] = inode; R[inode] = node.vals[U::Dof]; R[inode + nnode] = node.vals[P::Dof]; R[inode + nnode * 2] = node.vals[V::Dof]; inode += 1; end for element in self.elements minode::Dict{Int32, Int32} = Dict{Int32, Int32}(); einode::Int32 = 1; ennode = length(element.nodes); #建立单元的节点自由度编号与总体节点自由度编号映射表 for node in element.nodes inode = mnode[node]; for dofn = 1 : getDofNum(node) minode[einode + (dofn - 1) * ennode] = inode + (dofn - 1) * nnode; end einode += 1; end #总体刚度和质量矩阵的集成(使用得到的自由度编号映射) Ke = getStiffMatrix(element); Me = getMassMatrix(element); (m, n) = size(Ke); for i = 1 : m mi = minode[i]; for j = 1 : n mj = minode[j]; K[mi, mj] += Ke[i, j]; M[mi, mj] += Me[i, j]; end end end #计算最终的刚度矩阵和载荷向量 R = (M - K * dt / 2) * R; K = (M + K * dt / 2); #施加约束 for constraint in self.constraints for node in constraint.nodes inode = mnode[node]; dofn::Int32 = 1; for d in instances(Dof) if haskey(constraint.values, d) idof = inode + (dofn - 1) * nnode; #对对角元增加罚函数 K[idof, idof] += 1e20; R[idof] = constraint.values[d] * K[idof, idof]; end dofn += 1; end end end #进行方程组的求解 self.result = K \ R; #更新节点上的数值 inode = 1; for node in self.nodes node.vals[U::Dof] = self.result[inode]; node.vals[P::Dof] = self.result[inode + nnode]; node.vals[V::Dof] = self.result[inode + nnode * 2]; inode += 1; end end
4. 圆柱绕流问题求解
圆柱绕流问题为流体力学中的经典问题，用以上完成的程序进行求解，可以验证程序的正确性
(1) 前处理
虽然Patran作为前处理软件来说并不优秀，但由于对BDF文件格式较为熟悉，故还是选择了Patran。通过几何建模得到一个带有圆孔的长方形区域，并生成三角形网格。对长方形区域的4个边界和圆孔的边界分别设置边界条件（具体边界条件的设置并不重要，主要是为了可以在生成的BDF文件中找到对应边界上的节点）。
虽然不需要用到Patran中的单元属性，但如果不对单元属性进行设置会无法得到完整的BDF文件，为此需要对单元设置一个随意的单元属性。
在分析中选择只输入分析模型文件，同时由于对BDF文件格式解析上的方便，选择small格式的卡片格式。
对BDF文件进行解析的代码如下
local nodes::Dict{Int32, Node} = Dict{Int32, Node}(); local elements::Array{Array{Int32}} = Array{Array{Int32}}([]); for line in eachline("A9.bdf") if line[1] != '$' local len = length(line); local n = len ÷ 8 + 1; local card = strip(line[1 : 8]); #解析单元和节点 if card == "CTRIA3" local nids::Array{Int32} = Array{Int32}([]); local id = parse(Int32, String(strip(line[9 : 16]))); for nid in split(strip(line[25 : end]), r"\s+") push!(nids, parse(Int32, String(nid))); end push!(elements, nids); elseif card == "GRID" #对数字的科学计数法进行处理，BDF文件中指数为负数时省略了e local x::Float64 = parse(Float64, replace(String(strip(line[25 : 32])), r"\-(\d+)$" => s"e-\1")); local y::Float64 = parse(Float64, replace(String(strip(line[33 : 40])), r"\-(\d+)$" => s"e-\1")); local id = parse(Int32, String(strip(line[9 : 16]))); node = Node(x, y); nodes[id] = node; end end end(2) 边界条件的设定
进口出处的节点设定为恒定速度u = 5，v = 0，上下边界处速度均设定为0，出口处设定为恒定速度u = 5，v = 0保证区域内的质量守恒。相关的代码如下
local spc1 = SPC(); addConstraint(spc1, U::Dof, 5.0); addConstraint(spc1, V::Dof, 0.0); addConstraint(spc1, P::Dof, 0.0); for i = 83 : 101 addNode(spc1, nodes[i]); nodes[i].vals[U::Dof] = 5.0; end for i = 108 : 200 addNode(spc1, nodes[i]); nodes[i].vals[U::Dof] = 5.0; end local spc2 = SPC(); addConstraint(spc2, U::Dof, 0.0); addConstraint(spc2, V::Dof, 0.0); addConstraint(spc2, P::Dof, 0.0); for i = 1 : 82 addNode(spc2, nodes[i]); end for i = 102 : 181 addNode(spc2, nodes[i]); end for i = 201 : 250 addNode(spc2, nodes[i]); end addConstraint(m, spc1); addConstraint(m, spc2);(3) 求解
求解的过程较为简单，只需要安装设定的时间补偿循环一定的次数即可，每步完成后记录相应的结果，为了方便进行后处理，结果保存为3个文件——节点信息文件，单元信息文件，计算结果文件
#通过命令行给定计算的步数 local n = parse(Int32, ARGS[1]); for i = 1 : n solve(m, 0.05); save(m, "data"); end(4) 后处理
由于Julia没有提供比较合适的绘图库（有的一种本质上也为调用的matplotlib），这里使用Python提供的matplotlib进行分析结果的可视化。
计算结果可以绘制流场图和速度云图，绘制流场图的代码如下
import matplotlib.pyplot as plt; x = []; y = []; u = []; v = []; p = []; m = []; for line in open('data.1'): row = line.strip().split(" "); x.append(float(row[0])); y.append(float(row[1])); for line in open('data.2'): row = line.strip().split(" "); u0 = float(row[0]); v0 = float(row[2]); u.append(u0); p.append(float(row[1])); v.append(v0); for line in open('data.3'): row = line.strip().split(" "); m.append([int(row[0]) - 1, int(row[1]) - 1, int(row[2]) - 1]); fig, ax = plt.subplots(); ax.set_aspect('equal'); q = ax.quiver(x, y, u, v); plt.savefig(sys.argv[1] + '.png');从以0.05的时间步长计算400步，生成的结果图中可以观察到，圆柱附近涡的形成、发展和脱落的过程，可参见下面的动图。
