基于FreeFEM++的有限元编程--4

2 第二类方程--- 弹性方程

2.1 弹性方程描述

在物理学上，弹性是指物体在外力作用下发生形变，当外力撤消后能恢复原来大小和形状的性质。弹性理论是描述一个物体在外力的作用下如何运动或发生形变，弹性方程是描述变形量与外力之间的关系，通常表示为σ=Eε，其中E表示杨氏模量。下面是线性弹性方程的变分形式：

 

下面脚本是对该问题比较粗糙的写法：

real L=10,H=1,dn=5;

mesh Sh=square(L*dn,H*dn,[x*L,y*H]); // square OK of rectangles

plot(Sh, wait=1);

real mu=1,lambda=1;

fespace Vh(Sh,P1);

Vh u1,u2,v1,v2; // One for each component

func g1=0; 

func g2=-1.;

problem elasticity(u1,u2,v1,v2)=

int2d(Sh)(2.*mu*(dx(u1)*dx(v1)+dy(u2)*dy(v2)+(dy(u1)+dx(u2))*(dy(v1)+dx(v2))/2.)

+ lambda*(dx(u1)+dy(u2))*(dx(v1)+dy(v2)))

-int1d(Sh,2)(g1*v1+g2*v2) // ARGH ....

+on(4,u1=0,u2=0);

elasticity;

plot(u1,cmm="u1 distribute", wait=1);

plot(u2,cmm="u2 distribute",  wait=2);

plot([u1,u2],wait=1); // you can plot vectors :)

 

 

 

 

比较高级的方法是采用宏macro，采用宏之后写的脚本更清晰易懂。

需要注意的是，每个宏都需要一个//结束符

real L=10,H=1,dn=5;

mesh Sh=square(L*dn,H*dn,[x*L,y*H]);

real mu=1,lambda=1;

fespace Vh(Sh,[P1,P1]); // vectorial FE space

macro u [u1,u2]     // end of macro

macro v [v1,v2]     // end of macro

macro e(u) [dx(u[0]),(dx(u[1])+dy(u[0]))/ sqrt(2),dy(u[1])] //

macro A [[2.*mu+lambda,0,lambda],[0,2.*mu,0],[lambda,0,2.*mu+lambda]]

//

macro g [0,-1.] //

Vh u,v; // One for each component

problem elasticity(u,v)=

    int2d(Sh)((A*e(u))'*e(v))-int1d(Sh,2)(g'*v)

+on(4,u1=0,u2=0);

elasticity;

plot(u,wait=1);