Дисциплина:
Раздел:
Выполнил:

Методы конечных элементов

Методы конечных элементов

ЧГУ, матфак, 4 курс Осипов.П.Г

www.domath.ru

>    # Подключим пакеты: пакет графики, пакет графического инструментария и пакет линейной алгебры

>    restart;

>    with(plots):
with(plottools):
with(LinearAlgebra):

  #входные данные (длина радиуса, угол, средние значения производных на соотвествующих отрезках)
>    R:=2:Fi:=Pi/18:

   unAF:=0: unFB:=0:
   unBH:=-1: unHC:=-1:

         unAK:=1/R: unKL:=1/R:

  #Согласно начальным данным зададим координаты узлов конечных элементов A,F,B,K и т.д., как на рисунке 1. Данные координаты задаются чисто из геометрических соображений.

>    A:=[-R*sin(Fi),R*cos(Fi)]:

F:=[-((R+1)/2)*sin(Fi),((R+1)/2)*cos(Fi)]:

B:=[-sin(Fi),cos(Fi)]:

K:=[0,R]:
E:=[0,(1+R)/2]:
H:=[0,1]:

L:=[R*sin(Fi),R*cos(Fi)]:

G:=[((R+1)/2)*sin(Fi),((R+1)/2)*cos(Fi)]:
C:=[sin(Fi),cos(Fi)]:

>    #Зададим локальную нумерацию узлов.
#Например, локальный узел x[1,4] совпадает с точкой H[1]
#(1-означает лок.нум., а 4-означает глоб.нум.
#Такая нумерация предложена на рисунке.

>    x[1,1]:=F[1]:
x[1,2]:=F[1]:
x[1,3]:=B[1]:
x[1,4]:=H[1]:
x[1,5]:=L[1]:
x[1,6]:=G[1]:
x[1,7]:=G[1]:
x[1,8]:=C[1]:


x[2,1]:=K[1]:
x[2,2]:=E[1]:
x[2,3]:=E[1]:
x[2,4]:=E[1]:
x[2,5]:=K[1]:
x[2,6]:=L[1]:
x[2,7]:=E[1]:
x[2,8]:=G[1]:

x[3,1]:=A[1]:
x[3,2]:=K[1]:
x[3,3]:=F[1]:
x[3,4]:=B[1]:
x[3,5]:=E[1]:
x[3,6]:=E[1]:
x[3,7]:=H[1]:
x[3,8]:=H[1]:

>    y[1,1]:=F[2]:
y[1,2]:=F[2]:
y[1,3]:=B[2]:
y[1,4]:=H[2]:
y[1,5]:=L[2]:
y[1,6]:=G[2]:
y[1,7]:=G[2]:
y[1,8]:=C[2]:


y[3,1]:=A[2]:
y[3,2]:=K[2]:
y[3,3]:=F[2]:
y[3,4]:=B[2]:
y[3,5]:=E[2]:
y[3,6]:=E[2]:
y[3,7]:=H[2]:
y[3,8]:=H[2]:


y[2,1]:=K[2]:
y[2,2]:=E[2]:
y[2,3]:=E[2]:
y[2,4]:=E[2]:
y[2,5]:=K[2]:
y[2,6]:=L[2]:
y[2,7]:=E[2]:
y[2,8]:=G[2]:

  #Введем некоторые элементы, значения, которых изменяются по правилу круговой перестановки.

>    for j from 1 to 8 do

    a[1,j]:=y[2,j]-y[3,j]:
    a[2,j]:=y[3,j]-y[1,j]:
    a[3,j]:=y[1,j]-y[2,j]:

     b[1,j]:=x[3,j]-x[2,j]:
     b[2,j]:=x[1,j]-x[3,j]:
     b[3,j]:=x[2,j]-x[1,j]:

#Вычислим площади всех треугольных элементов по координатам их вершин. Площадь #вычисляется, исходя из свойства определителя матрицы. Такая матрица составлена из координат #вершин соотвествующего треугольника.   

M[j]:=Matrix(3,3,[ [x[1,j],y[1,j],1], [x[2,j],y[2,j],1], [x[3,j],y[3,j],1] ]):
S[j]:=(1/2)*Determinant(M[j]):

    for k from 1 to 3 do
      for i from 1 to 3 do
 
#Введем значения интегралов.
           g[k,i,j]:=(a[i,j]*a[k,j]+b[i,j]*b[k,j])/(4*S[j]):
         od:
      od:
  od:

>    #Построение элементов матрицы глобальной системы уравнений
f[1,1]:=g[3,3,1]:
f[2,1]:=g[1,3,1]: f[2,2]:=g[1,1,1]+g[1,1,2]+g[3,3,3]:
f[3,1]:=0: f[3,2]:=g[1,3,3]:f[3,3]:=g[1,1,3]+g[3,3,4]:
f[4,1]:=g[2,3,1]:f[4,2]:=g[1,2,1]+g[1,3,2]:f[4,3]:=0:f[4,4]:=g[2,2,1]+g[3,3,2]+g[2,2,5]:
f[5,1]:=0: f[5,2]:=g[1,2,2]+g[2,3,3]:f[5,3]:=g[1,2,3]+g[2,3,4]:f[5,4]:=g[2,3,2]+g[2,3,5]:
f[5,5]:=g[2,2,2]+g[2,2,3]+g[2,2,4]+g[3,3,5]+g[3,3,6]+g[2,2,7]:
f[6,1]:=0: f[6,2]:=0: f[6,4]:=0: f[6,3]:=g[1,3,4]: f[6,5]:=g[1,2,4]+g[2,3,7]: f[6,6]:=g[1,1,4]+g[3,3,7]+g[3,3,8]:
f[7,1]:=0: f[7,2]:=0: f[7,3]:=0: f[7,6]:=0: f[7,4]:=g[1,2,5]: f[7,5]:=g[1,3,5]+g[2,3,6]: f[7,7]:=g[1,1,5]+g[2,2,6]:
f[8,1]:=0: f[8,2]:=0: f[8,3]:=0: f[8,4]:=0: f[8,5]:=g[1,3,6]+g[1,2,7]: f[8,6]:=g[1,3,7]+g[2,3,8]: f[8,7]:=g[1,2,6]: f[8,8]:=g[1,1,6]+g[1,1,7]+g[2,2,8]:
f[9,1]:=0: f[9,2]:=0: f[9,3]:=0: f[9,4]:=0: f[9,5]:=0: f[9,7]:=0: f[9,6]:=g[1,3,8]: f[9,8]:=g[1,2,8]: f[9,9]:=g[1,1,8]:

>    #Свойство системы глобальных уравнений - симметричность. Добиться симметричности можно легко прибегнув к циклу.

>    for t from 1 to 9 do
     for m from 1 to 9 do
                      f[t,m]:=f[m,t]:
      od:
   od:

>    #Процедура, которая вычисляет расстояние между двумя точками P,Q (входные параметры проц.)
l:=proc(P,Q)
               sqrt((P[1]-Q[1])^2+(P[2]-Q[2])^2);
 end proc:

  #Формирование правой части глобальной системы уравнений
>    z[1]:=(unAK*l(K,A)+unAF*l(A,F))/2:
z[2]:=(unAF*l(A,F)+unFB*l(F,B))/2:
z[3]:=(unFB*l(F,B)+unBH*l(B,H))/2:
z[4]:=(unAK*l(A,K)+unKL*l(K,L))/2:
z[5]:=0:
z[6]:=(unBH*l(B,H)+unHC*l(H,C))/2:
z[7]:=(unKL*l(K,L)+unLG*l(L,G))/2:
z[8]:=(unLG*l(L,G)+unGC*l(G,C))/2:
z[9]:=(unGC*l(G,C)+unHC*l(C,H))/2:

  # Процедура, которая вычисляет радиус точки T(T-входной параметр процедуры, T[1] - первая координата, T[2] - вторая).
>    r:=proc(T)
           sqrt((T[1])^2+(T[2])^2);
end proc:

  #Остальные данные для u[7], u[8], u[9 ] дополним граничными условиями.

#Функция evalf(f,k) вычисляет выражение f с точностью k. Для нас будет достаточно k=2.  

>    u[7]:=evalf(ln(r(L)),2):
u[8]:=evalf(ln(r(G)),2):
u[9]:=evalf(ln(r(C)),2):

  # Глобальная система уравнений FU=Z .  
>    for v from 1 to 9 do
                    e[v]:=sum(f[v,w]*u[w],w=1..9)-z[v]=0:
                    print(evalf(e[v])):
 od:

.9823122995*u[1]-.3062103226*u[2]-.6761019780*u[4]-.8715574278e-1 = 0.

-.3062103226*u[1]+2.444522142*u[2]-.2187216579*u[3]-.16e-9*u[4]-1.919590161*u[5] = 0.

-.2187216579*u[2]+1.702158610*u[3]-.96e-9*u[5]-1.483436952*u[6]+.8715574275e-1 = 0.

-.6761019780*u[1]-.16e-9*u[2]+1.964624602*u[4]-.6124206440*u[5]-.6429496658 = 0.

-1.919590161*u[2]-.96e-9*u[3]-.6124206440*u[4]+4.889044284*u[5]-.4374433186*u[6]-.7783268314 = 0.

-1.483436952*u[3]-.4374433186*u[5]+3.404317223*u[6]+.1743114856 = 0.

-.6761019785*u[4]+.96e-9*u[5]+.4695736571-.2500000000*unLG = 0.

-1.919590161*u[5]+.16e-9*u[6]+.7789196124-.2500000000*unLG-.2500000000*unGC = 0.

-1.483436953*u[6]-.152825828e-2-.2500000000*unGC = 0.

  #Решение глобальной системы уравнений относительно U. Решение может быть получено с помощью solve.   

>    for v from 1 to 6 do
                       u[v]:=evalf(solve(e[v],u[v])):
end do:


#Вывод на экран полученныхвычислений.

for v from 1 to 9 do
                          print(v , u[v]):
end do;

1, .6904734468

2, .4045558678

3, .2395173451e-2

4, .6910579289

5, .4047696324

6, .1852169983e-2

7, .6931471806

8, .4054651081

9, 0.

  #Вычисление точных значений Uточное=ln(r), r -радиус вектор точки.   

>    u1[1]:=evalf(ln(r(A))): u1[2]:=evalf(ln(r(F))): u1[3]:=evalf(ln(r(B))): u1[4]:=evalf(ln(r(K))): u1[5]:=evalf(ln(r(E))): u1[6]:=evalf(ln(r(H))): u1[7]:=evalf(ln(r(L))): u1[8]:=evalf(ln(r(G))): u1[9]:=evalf(ln(r(C))):

for v from 1 to 9 do
                     print(u1[v]):
end do;

.6931471806

.4054651081

0.

.6931471806

.4054651081

0.

.6931471806

.4054651081

0.

  #Вычисление погрешностей

>    delta:=sqrt(add((u1[v]-u[v])^2,v=1..9));

delta := .4689522256e-2

  #Зададим число узлов.

>    lam:=6:

>    v[1,2]:=u[2]: v[2,2]:=u[5]: v[3,2]:=u[4]: v[1,3]:=u[3]: v[2,3]:=u[5]: v[3,3]:=u[2]:
v[1,4]:=u[6]: v[2,4]:=u[5]: v[3,4]:=u[3]: v[1,5]:=u[7]: v[2,5]:=u[4]: v[3,5]:=u[5]:
v[1,6]:=u[8]: v[2,6]:=u[7]: v[3,6]:=u[5]: v[1,7]:=u[8]: v[2,7]:=u[5]: v[3,7]:=u[6]:

>    for j from 2 to 7 do
  
         Vx[j]:=evalf(add(a[i,j]*v[i,j],i=1..3)/(2*S[j]));
         Vy[j]:=evalf(add(b[i,j]*v[i,j],i=1..3)/(2*S[j]));

 od:

>    VX:=evalf((1/lam)*sum(Vx[r],r=2..7));

>    VY:=evalf((1/lam)*sum(Vy[r],r=2..7));

VX := -.1702332775e-2

VY := .6893645715

>    V:=evalf(sqrt(VX^2+VY^2));

>    VT:=evalf(2/(1+R));

V := .6893666733

VT := .6666666667

>    deltaV:=abs(V-VT);

deltaV := .227000066e-1

>    #Далее следует построение объектов для графического отображения полученного решения.

>    l0 :=point([-R*sin(Fi),R*cos(Fi),0], color=red,symbolsize=150),
     point([-((R+1)/2)*sin(Fi),((R+1)/2)*cos(Fi),0], color=red,symbolsize=150),
     point([-sin(Fi),cos(Fi),0], color=red, symbolsize=150),
     point([0,R,0], color=red, symbolsize=150),
     point([0,(1+R)/2,0], color=red, symbolsize=150),
     point([0,1,0], color=red,  symbolsize=150),
     point([R*sin(Fi),R*cos(Fi),0], color=red,  symbolsize=150),
     point([((R+1)/2)*sin(Fi),((R+1)/2)*cos(Fi),0], color=red,  symbolsize=150),
     point([sin(Fi),cos(Fi),0], color=red, symbolsize=150),
     point([0,0,0],  color=red, symbolsize=150):

>    l :=point([-R*sin(Fi),R*cos(Fi),u1[1]], color=green,symbol=CROSS,
                                               symbolsize=15 ),
     point([-((R+1)/2)*sin(Fi),((R+1)/2)*cos(Fi),u1[2]], color=green,
                                            symbol=CROSS, symbolsize=15 ),
     point([-sin(Fi),cos(Fi),u1[3]],
                                            symbol=CROSS, symbolsize=15 ),
     point([0,R,u1[4]], color=green,
                                            symbol=CROSS, symbolsize=15 ),
     point([0,(1+R)/2,u1[5]], color=green,
                                             symbol=CROSS, symbolsize=15 ),
     point([0,1,u1[6]], color=green,
                                             symbol=CROSS, symbolsize=15 ),
     point([R*sin(Fi),R*cos(Fi),u1[7]], color=green,
                                             symbol=CROSS, symbolsize=15 ),
     point([((R+1)/2)*sin(Fi),((R+1)/2)*cos(Fi),u1[8]], color=green,
                                             symbol=CROSS, symbolsize=15 ),
     point([sin(Fi),cos(Fi),u1[9]], color=green,
                                              symbol=CROSS, symbolsize=15):

>    ll :=point([-R*sin(Fi),R*cos(Fi),u[1]],color=yellow,symbol=CROSS, symbolsize=25),
     point([-((R+1)/2)*sin(Fi),((R+1)/2)*cos(Fi),u[2]],color=yellow,symbol=CROSS, symbolsize=25),
     point([-sin(Fi),cos(Fi),u[3]],color=yellow,symbol=CROSS, symbolsize=25),
     point([0,R,u[4]],color=yellow,symbol=CROSS, symbolsize=25),
     point([0,(1+R)/2,u[5]], color=yellow,symbol=CROSS, symbolsize=25),
     point([0,1,u[6]], color=yellow,symbol=CROSS, symbolsize=25),
     point([R*sin(Fi),R*cos(Fi),u[7]], color=yellow,symbol=CROSS, symbolsize=25),
     point([((R+1)/2)*sin(Fi),((R+1)/2)*cos(Fi),u[8]], color=yellow,symbol=CROSS, symbolsize=25),
     point([sin(Fi),cos(Fi),u[9]], color=yellow,symbol=CROSS, symbolsize=25):

>    p:=textplot3d([[-R*sin(Fi),R*cos(Fi),0,`a`],
[-((R+1)/2)*sin(Fi),((R+1)/2)*cos(Fi),0,`b`],
[-sin(Fi),cos(Fi),0,`c`],
[0,R,0,`d`],
[0,(1+R)/2,0,`e`],
[0,1,0,`j`],
[R*sin(Fi),R*cos(Fi),0,`f`],
[((R+1)/2)*sin(Fi),((R+1)/2)*cos(Fi),0,`g`],
[sin(Fi),cos(Fi),0,`h`]
],color=black):

>    display(l0,l,ll,p, axes=boxed, view=[-0.5..0.5, -0.2..2.1, -0.2..0.8], orientation=[165,65]);

[Maple Plot]