空气动力学数值方法：有限元法(FEM)：有限元法概论

空气动力学数值方法：有限元法(FEM)：有限元法概论1空气动力学基础1.1流体力学基本概念流体力学是研究流体（液体和气体）的运动和静止状态的学科。在空气动力学中，我们主要关注气体的行为，尤其是空气。流体的基本概念包括：密度（ρ）:单位体积的流体质量。压力（P）:流体垂直作用于单位面积上的力。速度（v）:流体在某一点的运动速度。温度（T）:流体的热状态，影响其物理性质如密度和粘度。粘度（μ）:流体内部摩擦力的度量，影响流体流动的阻力。1.2连续介质假设在空气动力学中，连续介质假设是一个基本的理论前提，它认为流体可以被视为连续的、无间隙的介质，而不是由离散的分子组成。这一假设允许我们使用连续函数来描述流体的物理量，如速度、压力和温度，从而简化了流体动力学方程的求解。例如，我们可以使用偏微分方程来描述流体的运动，而不是处理复杂的分子动力学。1.3纳维-斯托克斯方程纳维-斯托克斯方程是描述粘性流体运动的偏微分方程组。在空气动力学中，这些方程是理解和预测空气流动的关键。方程的一般形式为：ρ其中，v是流体的速度矢量，ρ是流体的密度，μ是流体的粘度，f是作用在流体上的外力，如重力。这个方程描述了流体的动量守恒。1.3.1示例代码下面是一个使用Python和SciPy库来求解二维纳维-斯托克斯方程的简化示例。这个例子使用了有限差分方法，而不是有限元法，但可以帮助理解数值求解流体动力学方程的基本步骤。importnumpyasnp

fromscipy.sparseimportdiags

fromscipy.sparse.linalgimportspsolve

#定义网格大小和时间步长

nx,ny=100,100

dx,dy=1.0/(nx-1),1.0/(ny-1)

nt=100

dt=0.01

#定义流体的物理属性

rho=1.0

mu=0.1

#初始化速度和压力场

u=np.zeros((ny,nx))

v=np.zeros((ny,nx))

p=np.zeros((ny,nx))

#定义边界条件

u[0,:]=1.0

#构建拉普拉斯算子

A=diags([-2,1,1],[0,-1,1],shape=(ny-2,ny-2)).toarray()/(dx*dy)

A[0,:]=0

A[-1,:]=0

A[0,0]=1/dy

A[-1,-1]=1/dy

#时间迭代

forninrange(nt):

un=u.copy()

vn=v.copy()

#更新速度场

u[1:-1,1:-1]=un[1:-1,1:-1]-un[1:-1,1:-1]*dt/dx*(un[1:-1,1:-1]-un[1:-1,0:-2])\

-vn[1:-1,1:-1]*dt/dy*(un[1:-1,1:-1]-un[0:-2,1:-1])\

-dt/(rho*dx*dy)*(p[1:-1,2:]-p[1:-1,0:-2])\

+mu*(dt/dx**2+dt/dy**2)*(un[1:-1,2:]-2*un[1:-1,1:-1]+un[1:-1,0:-2]\

+un[2:,1:-1]-2*un[1:-1,1:-1]+un[0:-2,1:-1])

v[1:-1,1:-1]=vn[1:-1,1:-1]-un[1:-1,1:-1]*dt/dx*(vn[1:-1,1:-1]-vn[1:-1,0:-2])\

-vn[1:-1,1:-1]*dt/dy*(vn[1:-1,1:-1]-vn[0:-2,1:-1])\

-dt/(rho*dx*dy)*(p[2:,1:-1]-p[0:-2,1:-1])\

+mu*(dt/dx**2+dt/dy**2)*(vn[1:-1,2:]-2*vn[1:-1,1:-1]+vn[1:-1,0:-2]\

+vn[2:,1:-1]-2*vn[1:-1,1:-1]+vn[0:-2,1:-1])

#应用边界条件

u[0,:]=1.0

u[-1,:]=0.0

u[:,0]=0.0

u[:,-1]=0.0

v[0,:]=0.0

v[-1,:]=0.0

v[:,0]=0.0

v[:,-1]=0.0

#更新压力场

b=np.zeros((ny-2,nx-2))

b[1:-1,1:-1]=-rho*(1/dt*(un[1:-1,2:]-un[1:-1,0:-2])/(2*dx)\

+1/dt*(vn[2:,1:-1]-vn[0:-2,1:-1])/(2*dy))

p[1:-1,1:-1]=spsolve(diags([-1,1,1,-1],[0,-1,1,-nx+1],shape=(ny-2,ny-2)),b.flatten()).reshape(ny-2,nx-2)1.3.2代码解释这段代码首先定义了网格的大小、时间步长和流体的物理属性。然后，它初始化了速度和压力场，并设置了边界条件。在时间迭代中，它使用了显式有限差分方法来更新速度场，然后通过求解泊松方程来更新压力场，以满足不可压缩流体的连续性方程。1.4边界条件与初始条件在求解流体动力学方程时，边界条件和初始条件是至关重要的。边界条件描述了流体在边界上的行为，如速度、压力或温度的值。初始条件则指定了在计算开始时流体的状态。1.4.1边界条件示例在上述代码示例中，我们设置了以下边界条件：上边界（y=0）:下边界（y=1）:左边界（x=0）:u=右边界（x=1）:u=1.4.2初始条件示例在代码示例中，初始条件是所有速度和压力场的值都设为零，除了上边界的速度设为1.0。这表示流体开始时静止，除了在上边界有一个恒定的入流速度。通过这些边界条件和初始条件，我们可以设置一个简单的流体流动场景，并使用数值方法来求解流体动力学方程，从而预测流体的行为。2有限元法原理2.1有限元法的历史与应用有限元法（FiniteElementMethod,FEM）起源于20世纪40年代末，最初由工程师们在解决结构工程问题时提出。1943年，R.Courant在解决弹性问题时首次使用了类似于有限元法的离散化技术。然而，直到1956年，O.C.Zienkiewicz和Y.K.Cheung在《工程计算中的有限元法》一文中详细阐述了有限元法的理论基础，这一方法才开始被广泛接受和应用。有限元法的应用领域广泛，包括但不限于：-结构分析：评估桥梁、建筑物、飞机等结构的强度和稳定性。-流体动力学：模拟空气动力学、水动力学等流体流动问题。-热传导分析：预测材料在不同温度下的热传导行为。-电磁学：分析电磁场的分布和电磁设备的性能。2.2离散化过程2.2.1原理有限元法的核心在于将连续的物理域离散化为有限数量的子域，即“有限元”。每个子域（或单元）的物理行为通过简单的数学模型来近似，这些模型通常基于多项式函数。通过组合所有单元的模型，可以构建整个物理域的近似解。2.2.2步骤网格划分：将分析区域划分为多个单元，每个单元的形状可以是三角形、四边形、六面体等。选择基函数：在每个单元内，选择适当的基函数来表示未知函数的近似解。建立单元方程：利用基函数和单元的物理定律，建立每个单元的方程。组装全局方程：将所有单元方程组合成一个全局方程系统。求解方程：使用数值方法求解全局方程，得到整个物理域的解。2.3加权残值法2.3.1原理加权残值法是有限元法中用于建立弱形式方程的一种方法。它基于残差（即方程的不满足程度）的加权积分，通过选择适当的加权函数，使残差在整体上最小化，从而得到方程的近似解。2.3.2应用考虑一个简单的线性微分方程：d在有限元法中，我们首先将方程的强形式转化为弱形式，通过乘以一个测试函数（或加权函数）并进行积分，得到：Ω其中，v是测试函数，Ω是分析区域。2.3.3代码示例假设我们使用Python的SciPy库来求解上述方程的一个具体实例，其中fx=x，边界条件为uimportnumpyasnp

fromscipy.sparseimportdiags

fromscipy.sparse.linalgimportspsolve

#定义网格参数

n=100#网格点数

h=1/(n+1)#网格步长

#创建网格

x=np.linspace(0,1,n+2)[1:-1]

#定义微分方程的右侧函数

deff(x):

returnx

#创建矩阵A和向量b

data=np.array([[-1,2,-1]]).repeat(n,axis=0)

A=diags(data,[-1,0,1],shape=(n,n)).toarray()

b=h**2*f(x)

#应用边界条件

A[0,0]=1

A[-1,-1]=1

b[0]=0

b[-1]=1

#求解方程

u=spsolve(A,b)

#输出结果

print(u)2.4有限元方程的建立2.4.1原理有限元方程的建立基于变分原理或能量原理，通过最小化能量泛函来求解物理问题。对于线性弹性问题，能量泛函通常表示为应变能和外力做功的差。2.4.2应用考虑一个弹性体的平衡方程，其能量泛函为：Π其中，σ是应力张量，ε是应变张量，b是体力，u是位移，t是表面力。2.4.3代码示例使用Python和FEniCS库来建立并求解一个简单的弹性问题的有限元方程。fromfenicsimport*

#创建网格和函数空间

mesh=UnitSquareMesh(8,8)

V=VectorFunctionSpace(mesh,'P',1)

#定义边界条件

defboundary(x,on_boundary):

returnon_boundary

bc=DirichletBC(V,Constant((0,0)),boundary)

#定义变分问题

u=TrialFunction(V)

v=TestFunction(V)

f=Constant((0,-10))

T=Constant((1,0))

a=dot(grad(u),grad(v))*dx

L=dot(f,v)*dx+dot(T,v)*ds

#求解方程

u=Function(V)

solve(a==L,u,bc)

#输出结果

plot(u)

interactive()以上代码示例展示了如何使用FEniCS库来建立和求解一个二维弹性问题的有限元方程。通过定义网格、函数空间、边界条件、变分形式和求解器，可以得到位移场的数值解，并通过绘图功能直观展示结果。3有限元法在空气动力学中的应用3.1网格生成技术网格生成是有限元分析中的关键步骤，它将连续的物理域离散化为一系列有限的、互不重叠的子域，即单元。在空气动力学中，网格的质量直接影响到计算的准确性和效率。网格可以是结构化的，如矩形网格，或非结构化的，如三角形和四面体网格。3.1.1生成结构化网格结构化网格通常用于形状规则的区域，如圆柱或矩形。生成结构化网格的一个简单方法是使用均匀分割。例如，考虑一个二维矩形区域，我们可以将其分割为一系列矩形单元。#生成结构化网格的示例代码

importnumpyasnp

#定义矩形区域的边界

x_min,x_max=0,1

y_min,y_max=0,1

#定义网格的分辨率

nx,ny=10,10

#生成网格

x=np.linspace(x_min,x_max,nx)

y=np.linspace(y_min,y_max,ny)

X,Y=np.meshgrid(x,y)

#打印网格节点坐标

print("网格节点坐标:")

print(X,Y)3.1.2生成非结构化网格非结构化网格适用于复杂几何形状，如飞机或汽车的表面。生成非结构化网格通常需要使用专门的网格生成软件，如Gmsh或TetGen。这里以Gmsh为例，展示如何生成非结构化网格。//Gmsh脚本示例

Point(1)={0,0,0,1.0};

Point(2)={1,0,0,1.0};

Point(3)={1,1,0,1.0};

Point(4)={0,1,0,1.0};

Line(1)={1,2};

Line(2)={2,3};

Line(3)={3,4};

Line(4)={4,1};

LineLoop(5)={1,2,3,4};

PlaneSurface(6)={5};

//设置网格细化程度

Mesh.CharacteristicLengthMin=0.1;

Mesh.CharacteristicLengthMax=0.1;

//生成网格

Mesh.Algorithm=6;

Mesh.Algorithm3D=1;

Mesh.Optimize=1;

Mesh.HighOrderOptimize="Standard";3.2空气动力学问题的建模空气动力学问题的建模涉及将流体动力学方程（如Navier-Stokes方程）转化为有限元形式。这通常包括选择适当的基函数、定义边界条件和应用加权残值法。3.2.1基函数的选择基函数用于表示解的近似形式。在空气动力学中，常用的基函数包括线性、二次和高阶多项式。选择基函数时，需要考虑问题的复杂性和计算资源。3.2.2边界条件的定义边界条件描述了流体在边界上的行为，如速度、压力或热流。在有限元法中，边界条件的正确设置对于获得准确的解至关重要。3.2.3加权残值法的应用加权残值法是将微分方程转化为积分形式，然后在每个单元上应用。这涉及到选择权重函数，通常是基函数的变体，以确保方程在每个单元上都得到满足。3.3求解器与算法有限元法求解空气动力学问题时，需要使用高效的求解器和算法。这包括直接求解器、迭代求解器和多网格方法。3.3.1直接求解器直接求解器通过求解线性方程组来获得解。对于小规模问题，直接求解器可以提供准确的解，但计算成本较高。3.3.2迭代求解器迭代求解器通过一系列近似解逐步逼近真实解。对于大规模问题，迭代求解器通常更有效，但需要更精细的参数调整。3.3.3多网格方法多网格方法结合了不同分辨率的网格，以加速迭代求解器的收敛。这种方法在处理复杂流体动力学问题时特别有效。3.4后处理与结果分析后处理是有限元分析的最后一步，它涉及可视化和分析计算结果。这包括流线、压力分布和速度矢量的可视化，以及计算升力、阻力等空气动力学参数。3.4.1可视化流线流线是流体流动的轨迹，可以使用流体速度场来计算。在后处理中，流线的可视化可以帮助理解流体的流动模式。#使用matplotlib可视化流线的示例代码

importmatplotlib.pyplotasplt

importnumpyasnp

#生成速度场数据

x=np.linspace(0,1,10)

y=np.linspace(0,1,10)

X,Y=np.meshgrid(x,y)

U=np.cos(2*np.pi*X)*np.sin(2*np.pi*Y)

V=-np.sin(2*np.pi*X)*np.cos(2*np.pi*Y)

#绘制流线

plt.streamplot(X,Y,U,V)

plt.title('流线图')

plt.show()3.4.2分析空气动力学参数计算升力和阻力等空气动力学参数是有限元分析的重要部分。这通常涉及到积分压力和剪切力分布。#计算升力和阻力的示例代码

importnumpyasnp

#假设已知压力和剪切力分布

pressure=np.array([1,2,3,4,5])

shear_force=np.array([0.1,0.2,0.3,0.4,0.5])

#计算升力和阻力

lift=np.sum(pressure*np.sin(np.deg2rad(90)))+np.sum(shear_force*np.cos(np.deg2rad(90)))

drag=np.sum(pressure*np.cos(np.deg2rad(90)))+np.sum(shear_force*np.sin(np.deg2rad(90)))

print("升力:",lift)

print("阻力:",drag)以上内容详细介绍了有限元法在空气动力学中的应用，包括网格生成技术、空气动力学问题的建模、求解器与算法，以及后处理与结果分析。通过这些步骤，可以有效地分析和预测流体在复杂几何形状上的行为。4高级有限元法技术4.1自适应网格细化(subdir4.1)自适应网格细化是有限元分析中一种提高求解精度和效率的技术。它基于对解的局部误差估计，动态调整网格的密度，以在需要更高精度的区域增加网格密度，而在误差较小的区域保持较低的网格密度。这种方法可以显著减少计算资源的消耗，同时保持或提高求解的准确性。4.1.1原理自适应网格细化通常包括以下步骤：1.初始网格生成：创建一个粗网格作为求解的起点。2.求解：在当前网格上进行有限元求解。3.误差估计：评估解的局部误差，这可以通过后验误差估计器完成。4.网格调整：根据误差估计结果，细化误差较大的区域的网格，同时可能在误差较小的区域进行网格粗化。5.重复：重复步骤2至4，直到满足预定的误差标准或达到计算资源的限制。4.1.2示例假设我们正在使用Python的FEniCS库进行自适应网格细化。以下是一个简化示例，展示如何在求解一个简单的偏微分方程时应用自适应网格细化：fromfenicsimport*

#创建一个初始网格

mesh=UnitSquareMesh(8,8)

#定义函数空间

V=FunctionSpace(mesh,'P',1)

#定义边界条件

defboundary(x,on_boundary):

returnon_boundary

bc=DirichletBC(V,Constant(0),boundary)

#定义变分问题

u=TrialFunction(V)

v=TestFunction(V)

f=Constant(1)

a=dot(grad(u),grad(v))*dx

L=f*v*dx

#求解

u=Function(V)

solve(a==L,u,bc)

#自适应网格细化

error_control=AdaptiveMeshRefinementControl()

error_control.error_control_parameter=0.5

error_control.marking_fraction=0.5

#重复求解和网格细化

foriinrange(5):

error_estimate=compute_error_estimate(u)

mesh=error_control.adapt(mesh,error_estimate)

V=FunctionSpace(mesh,'P',1)

bc=DirichletBC(V,Constant(0),boundary)

u=Function(V)

solve(a==L,u,bc)在这个示例中，我们首先创建了一个8x8的初始网格，然后定义了函数空间、边界条件和变分问题。通过AdaptiveMeshRefinementControl类，我们控制了网格的自适应细化过程，重复求解和细化网格5次，以达到更高的精度。4.2多尺度分析(subdir4.2)多尺度分析是处理具有多个尺度特征的复杂问题的一种方法。在空气动力学中，这可能涉及到从微观尺度的流体动力学到宏观尺度的气动性能分析。多尺度分析通过在不同尺度上使用不同的模型和方法，能够更准确地模拟和预测空气动力学现象。4.2.1原理多尺度分析的核心在于将问题分解为不同尺度的子问题，然后通过适当的耦合策略将这些子问题的解组合起来。这通常涉及到：1.尺度分离：识别问题中的不同尺度，并为每个尺度建立独立的模型。2.尺度间耦合：定义尺度间的信息传递机制，确保不同尺度的模型能够相互作用。3.尺度间迭代：在不同尺度的模型之间进行迭代求解，直到达到收敛。4.2.2示例在多尺度分析中，一个常见的例子是使用微观尺度的LatticeBoltzmann方法（LBM）来模拟流体动力学，同时在宏观尺度上使用有限元法来分析结构的响应。以下是一个简化示例，展示如何在Python中使用LBM和FEniCS库进行多尺度分析：#LBM模拟流体动力学

#假设我们使用一个LBM库，如lbmpy

#lbmpy_simulation=lbmpy.LBM(...)

#lbmpy_simulation.run()

#使用FEniCS进行结构响应分析

fromfenicsimport*

#创建一个结构网格

mesh=UnitSquareMesh(32,32)

#定义函数空间

V=VectorFunctionSpace(mesh,'P',2)

#定义边界条件

defboundary(x,on_boundary):

returnon_boundary

bc=DirichletBC(V,Constant((0,0)),boundary)

#定义变分问题

u=TrialFunction(V)

v=TestFunction(V)

f=Constant((0,-1))

a=inner(grad(u),grad(v))*dx

L=dot(f,v)*dx

#求解

u=Function(V)

solve(a==L,u,bc)

#将LBM的流体动力学结果作为边界条件输入到FEniCS中

#假设lbmpy_simulation.get_velocity()返回流体速度

velocity=lbmpy_simulation.get_velocity()

bc=DirichletBC(V,velocity,boundary)

solve(a==L,u,bc)在这个示例中，我们首先使用LBM库模拟流体动力学，然后使用FEniCS库分析结构响应。通过将LBM的流体速度作为边界条件输入到FEniCS中，我们实现了微观尺度和宏观尺度之间的耦合。4.3并行计算(subdir4.3)并行计算是利用多处理器或分布式计算资源来加速有限元分析的一种方法。在处理大型或复杂问题时，通过并行化可以显著减少计算时间。4.3.1原理并行计算通常基于以下策略：1.数据并行：将问题的数据集分割到多个处理器上，每个处理器处理数据集的一部分。2.任务并行：将问题分解为多个独立的任务，每个任务可以在不同的处理器上并行执行。3.混合并行：结合数据并行和任务并行，以实现更高效的并行计算。4.3.2示例使用Python的Dolfinx库，我们可以轻松地实现有限元分析的并行计算。以下是一个示例，展示如何在多个处理器上并行求解一个偏微分方程：fromdolfinximportfem,mesh,io,log

importufl

importmpi4py.MPIasmpi

#创建并行网格

domain=mesh.create_unit_square(mpi.COMM_WORLD,32,32)

#定义函数空间

V=fem.FunctionSpace(domain,("Lagrange",1))

#定义边界条件

defboundary(x):

returnnp.logical_or(np.isclose(x[0],0),np.isclose(x[0],1))

boundary_dofs=fem.locate_dofs_geometrical(V,boundary)

bc=fem.dirichletbc(mpi.COMM_WORLD,(0.0,0.0),boundary_dofs)

#定义变分问题

u=ufl.TrialFunction(V)

v=ufl.TestFunction(V)

f=fem.Constant(domain,(0.0,-1.0))

a=ufl.inner(ufl.grad(u),ufl.grad(v))*ufl.dx

L=ufl.dot(f,v)*ufl.dx

#求解

problem=fem.petsc.LinearProblem(a,L,bcs=[bc],petsc_options={"ksp_type":"preonly","pc_type":"lu"})

uh=problem.solve()在这个示例中，我们使用Dolfinx库创建了一个并行网格，并定义了函数空间、边界条件和变分问题。通过mpi4py库，我们能够利用多个处理器并行求解问题，从而加速计算过程。4.4高阶元与谱元方法(subdir4.4)高阶元和谱元方法是有限元法的扩展，它们通过使用更高阶的多项式基函数来提高解的精度。谱元方法特别适用于求解具有光滑解的问题，可以提供极高的精度。4.4.1原理高阶元和谱元方法的关键在于：1.基函数选择：使用更高阶的多项式基函数，如Lagrange多项式或Chebyshev多项式。2.积分规则：采用高精度的数值积分规则，如Gauss-Legendre积分，以确保变分形式的准确计算。3.求解策略：由于高阶元和谱元方法的矩阵通常更密集，可能需要使用更高效的求解器，如直接求解器或迭代求解器。4.4.2示例使用Python的FEniCS库，我们可以选择高阶的多项式基函数来提高求解的精度。以下是一个示例，展示如何使用高阶元求解一个偏微分方程：fromfenicsimport*

#创建一个网格

mesh=UnitSquareMesh(32,32)

#定义高阶函数空间

V=FunctionSpace(mesh,'CG',5)#使用5阶连续Galerkin元

#定义边界条件

defboundary(x,on_boundary):

returnon_boundary

bc=DirichletBC(V,Constant(0),boundary)

#定义变分问题

u=TrialFunction(V)

v=TestFunction(V)

f=Constant(1)

a=dot(grad(u),grad(v))*dx

L=f*v*dx

#求解

u=Function(V)

solve(a==L,u,bc)在这个示例中，我们使用了5阶的连续Galerkin元（CG）来定义函数空间，这将提供比标准1阶元更高的精度。通过选择适当的基函数和积分规则，我们可以有效地提高有限元分析的准确性，尤其是在处理具有复杂几何或高精度要求的问题时。5案例研究与实践5.1飞机翼型分析在飞机设计中，翼型的空气动力学性能至关重要。有限元法(FEM)可以用于分析翼型在不同飞行条件下的应力、应变和流体动力学特性。例如，考虑一个NACA0012翼型，在特定的飞行速度和攻角下进行分析。5.1.1数据样例NACA0012翼型的几何参数可以通过以下公式计算：y其中，t是翼型厚度，x是沿翼弦的位置，c是翼弦长度，t/5.1.2代码示例#导入必要的库

importnumpyasnp

importmatplotlib.pyplotasplt

#定义NACA0012翼型的几何参数

defnaca0012(x,t=0.12):

c=1.0#翼弦长度

yt=t/0.2*(0.2969*np.sqrt(x/c)-0.1260*(x/c)-0.3516*(x/c)**2+0.2843*(x/c)**3-0.1015*(x/c)**4)

returnyt

#生成翼型的坐标点

x=np.linspace(0,1,100)

y_upper=naca0012(x)

y_lower=-naca0012(1-x)

#绘制翼型

plt.figure(figsize=(10,5))

plt.plot(x,y_upper,'b',label='上表面')

plt.plot(1-x,y_lower,'r',label='下表面')

plt.legend()

plt.title('NACA0012翼型')

plt.xlabel('x位置')

plt.ylabel('y厚度')

plt.grid(True)

plt.show()5.2涡轮机叶片设计涡轮机叶片的设计需要精确计算叶片在高速旋转下的应力和变形，以确保其结构安全和效率。FEM可以模拟叶片在不同工况下的行为，帮助优化设计。5.2.1数据样例涡轮机叶片的几何数据通常包括叶片的轮廓、厚度分布和材料属性。例如，叶片的材料可能为钛合金，其弹性模量为110GPa，泊松比为0.34。5.2.2代码示例#导入有限元分析库

importfenics

#定义几何和网格

mesh=fenics.UnitSquareMesh(10,10)

#定义材料属性

E=110e9#弹性模量

nu=0.34#泊松比

mu=E/(2*(1+nu))

lmbda=E*nu/((1+nu)*(1-2*nu))

#定义边界条件

defboundary(x,on_boundary):

returnon_boundary

bc=fenics.DirichletBC(