Расчет прогиба пластинки под действием гидродинамических сил в протяженном канале

Санкт-Петербургский  Государственный Политехнический  Университет

Кафедра Гидроаэродинамики

 

 

 

 

 

 

 

Курсовая работа

Расчет прогиба  пластинки под действием гидродинамических  сил в протяженном канале

 

 

 

 

 

 

Работу  выполнил:

Андреев Виктор.

Группа  6054/11

 

 

 

 

 

 

 

 

 

Санкт-Петербург 2011

1. Постановка задачи

1.1. Общая формулировка

Рисунок 1. Схематическая постановка задачи

Рассматривается пластинка, разделяющая два протяженных, плоских канала с течением Пуазейля (Рисунок 1). В общем случае среды, текущие по двум каналам, могут отличаться друг от друга. Для простоты задачи будем считать, что среды одинаковые и, как следствие: ρ₂ = ρ₁ и ϻ₂= ϻ₁.

В ходе решения задачи использовались следующие параметры:

• длина канала и пластинки L=0,08 м.
• высота нижнего канала H₂=0,02 м
• высота верхнего канала H₁=0,04 м
• толщина пластинки δ=0,01 м
• среднерасходная скорость на входе в каждый канал,  u=0,008 м/с
• модуль Юнга .

Граничные условия:

• Вход – развитый профиль скорости V=V(y),
• Выход – р=0

 

Будем рассматривать задачу, в которой правый конец свободный, а левый закреплен.

Несмотря на то, что в  двух каналах среда одинаковая, при  разных параметрах Q и H возникают условия, когда давления сверху и снизу пластины разные. Ожидается, что давление в нижнем канале будет больше, чем давление в верхнем канале. Из-за этого возникает перепад давления (Рисунок 3), благодаря которому пластинка будет отклоняться от своего начального положения (Рисунок 2).

Рисунок 2. Прогиб пластины под действием разности давлений.

Рассматриваются две  постановки задачи:

1). W(x=L)<<H₁,H₂

2). W(x=L)=О(H₁,H₂)

Сначала будем  рассматривать случай, когда прогиб пластинки мал по сравнению с высотой каналов, следовательно, его влиянием на течения в каналах можно пренебречь (однонаправленная задача). Решение однонаправленной задачи состоит из следующих этапов:

• решается гидродинамическая задача о течении в каналах,
• находится перепад давлений ,
• решается задача о прогибе пластины

1.2. Течение в каналах

Из решения  задачи о течении Пуазейля [2]:

На рисунке 3 представлены зависимости давления сверху и снизу пластины от продольной координаты, рассчитанные с использованием аналитической формулы (1). 

Рисунок 3. Давление в каналах

Прогиб пластины определяется разницей давлений на верхней и нижней поверхностях пластины (Рисунок 3). Эта разница определяется соотношением .

Зная разницу  давлений на верхней и нижней поверхности  пластины можно найти прогиб пластины по формуле (3):

1.3. Задача о прогибе пластины

1.3.1. Общее уравнение

Рассмотрим модель тонкой пластинки. Уравнение моментов можно записать в виде:

где – модуль Юнга,

 – момент инерции, 

 – поперечная нагрузка,

 –момент силы.

Дифференцируя два раза первое уравнение системы (4) и подставляя вместо второй производной  от момента правую часть второго  уравнения, получим:

Рисунок 4. Действие растягивающей силы

Растягивающая сила (Рисунок 4):

При учете растягивающей  силы уравнение  обобщается:     

Численное решение  уравнения  производится с помощью пятидиагональной прогонки.

1.3.2. Граничные условия для пластинки

Вход - условие жесткой заделки

Выход - “мягкие” граничные условия

1.3.3. Аппроксимация уравнений

 

1.3.4. Аппроксимация граничных условий

1.3.5. Вид уравнения для прогонки

    , где

               

 

1.3.6. Аналитическое решение

Пусть T(x) – мала, вторым слагаемым в уравнении (8) можно пренебречь:

Интегрируя выражение  четыре раза, с учетом граничных условий и , получим

 

 

 

1.3.7. Метод прогонки для пятиточечных уравнений[1]

Наиболее часто  встречаются системы пятиточечных уравнений следующего вида:

Такого вида системы возникают при аппроксимации  краевых задач для обыкновенных дифференциальных уравнений четвертого порядка, а также при реализации разностных схем для уравнений в частных производных. Матрица A системы является пятидиагональной квадратичной матрицей размерности и имеет не более ненулевых элементов.

 

 

Для решения  системы  используем метод исключения Гаусса. Учитывая структуру системы, получим, что обратный ход метода Гаусса должен осуществляться по формулам

для .

 Для реализации  необходимо знать , а также определить коэффициенты .

Сначала найдем формулы для  . Используя , выразим через  
Получим:

Подставляя  и в , получим

Сравнивания это  выражение с  , видим, что если положить

где , то уравнения системы для будут удовлетворены.

Рекуррентные  соотношения  связывают с . Поэтому, если будут заданы для , то по формулам последовательно можем найти коэффициенты для .

Найдем  для . Из и формулы для непосредственно получим

 

 

 

Далее, подставляя значение при в , получим

Следовательно, будет выполнено, если положить

Используя , можно найти для . Осталось определить , входящие в формулу .

Воспользуемся для этого уравнениями и . Подставляя и при в и сравнивая полученное выражение с , найдем, что и определяются по формулам для . Найдем теперь . Для этого подставим при и в уравнение . Получим

или

где определяется по формуле при .

Объединяя полученные выше формулы, запишем алгоритм правой прогонки для системы в следующем виде:

1). по формулам

 

 

 

 

находятся прогоночные  коэффициенты.

2). Неизвестные  находятся последовательно по формулам

1.3.8. Код программы

Была написана программа, решающая уравнение (8) с использованием пятидиагональной прогонки. Код программы представлен ниже:

Program Progonka

 implicit none

  integer*4 i,n

  real*8 jung,Inercija,dx,mu,h1,h2,delta,u,length,coef

  real*8,allocatable:: a(:),b(:),c(:),d(:),e(:),f(:),T(:),alfa(:)

  real*8,allocatable:: beta(:),gamma(:),x(:),delta_calc(:),y(:),y_analitic(:),a_an(:) 

  character(*),parameter:: f1='E:\My_work_plastina\in.txt', f2='E:\My_work_plastina\out.txt'

 

  open(10,file=f1)  !чтение входных данных

  open(20,file=f2)  !вывод расчетных значений искомых параметров 

 

  read(10,*) jung

  read(10,*) h1

  read(10,*) h2

  read(10,*) delta

  read(10,*) n 

  read(10,*) U

  read(10,*) mu

  read(10,*) length

 

  allocate(a(n+1))

  allocate(b(n+1))

  allocate(c(n+1))

  allocate(d(n+1))

  allocate(e(n+1))

  allocate(f(n+1))

  allocate(T(n+1))

  allocate(alfa(n+1))

  allocate(beta(n+1))

  allocate(gamma(n+1))

  allocate(x(n+1))

  allocate(delta_calc(n+1))

  allocate(y(n+1))

  allocate(y_analitic(n+1))

  allocate(a_an(n+1)) 

 

  Inercija=(delta**3)/12

  x(1)=0

  dx=length/(n-1)

  do i=1,n

   !Координаты узлов сетки

 

 

   if(i/=1) then

    x(i)=x(i-1)+dx   

   endif

    f(i)=12.*mu*((u/h2/h2)-(u/h1/h1))*(length-x(i))

    T(i)=6.*mu*((U/h1)+(u/h2))*(length-x(i))

   enddo 

    !Коэффициенты

    coef=(jung*Inercija)/(dx**4)   

    C(1)=7.*coef+(2.*T(1)/dx/dx)

    D(1)=4.*coef+T(1)/dx/dx;

    E(1)=coef;   

    B(2)=4.*coef+(T(2)/dx/dx);

    C(2)=6.*coef+(2.*T(2)/dx/dx)

    D(2)=4.*coef+(T(2)/dx/dx);

    E(2)=coef;   

    do i=3,N-2   

        A(i)=coef;

        B(i)=4.*coef+(T(i)/dx/dx);

        C(i)=6.*coef+(T(i)*2./dx/dx)

        D(i)=4.*coef+(T(i)/dx/dx);

        E(i)=coef;       

    enddo

    A(N-1)=coef;

    B(N-1)=4.*coef+(T(N-1)/dx/dx)

    C(N-1)=5.*coef+(2.*T(N-1)/dx/dx)

    D(N-1)=2.*coef+(T(N-1)/dx/dx);   

    A(N)=coef;

    B(N)=2.*coef;

    C(N)=coef

 

   !Прямой ход прогонки

   alfa(2)=d(1)/c(1); beta(2)=e(1)/c(1); gamma(2)=f(1)/c(1) 

   alfa(3)=(d(2)-b(2)*beta(2))/(c(2)-b(2)*alfa(2))

   beta(3)=e(2)/(c(2)-b(2)*alfa(2))

   gamma(3)=(f(2)+b(2)*gamma(2))/(c(2)-b(2)*alfa(2))

   do i=3,n  

    delta_calc(i)=c(i)-a(i)*beta(i-1)+alfa(i)*(a(i)*alfa(i-1)-b(i))

    alfa(i+1)=(1d0/delta_calc(i))*(d(i)+beta(i)*(a(i)*alfa(i-1)-b(i))) 

    beta(i+1)=e(i)/delta_calc(i)

    gamma(i+1)=(1d0/delta_calc(i))*(f(i)-a(i)*gamma(i-1)-gamma(i)*(a(i)*alfa(i-1)-b(i)))

   enddo   

   !обратный ход прогонки

      y(1)=0

y(n)=gamma(n+1)

y(n-1)=alfa(n)*gamma(n+1)+gamma(n) 

    do i=n-2,2,-1          

      y(i)=alfa(i+1)*y(i+1)-beta(i+1)*y(i+2)+gamma(i+1)    

    enddo 

 !Вывод результатов

    do i=1,n   

     a_an(i)=12*mu*((u/h2/h2)-(u/h1/h1))/jung/Inercija

     y_analitic(i)=-a_an(i)*x(i)*x(i)*(10*x(i)*length*length-10*(length**3)-5*length*x(i)*x(i)+(x(i)**3))/120

     write(20,'(3(1pg17.10))') x(i),y(i),y_analitic(i)

    enddo

end

 

 

2. Анализ результатов расчетов.

2.1. Течение в каналах

Для анализа  течения в каналах использовался пакет ANSYS CFX. Размерность используемой сетки 100х44х4.

Задача решалась в стационарной постановке. Численная  схема – схема высоко порядка  High Resolution. Использовались следующие граничные условия:

Вход - развитое течение (формула ),

Выход - статическое давление равно нулю,

Боковые поверхности – условие симметрии.

Ширина канала в направлении  поперек течения .

Выражение для  профиля скорости, который ставился в качестве входного условия, имеет вид:

подставляя  сюда градиент давления из формулы  и учитывая, что , получим

Сравним разницу  давления на разных сторонах пластины, с распределением давления, полученным с помощью программного комплекса ANSYS CFX.

Рисунок 5. Разность давлений, действующих на пластину.

Из данного рисунка  видно, что решение, полученное с  помощью пакета Ansys CFX, отлично согласуется с решением, полученным при использовании аналитической формулы. Это говорит о правильности работы пакета Ansys CFX и возможности его использования при проведении дальнейших расчетов.

2.2. Прогиб пластины

2.2.1. Вопросы, связанные с чувствительностью решения к параметрам сетки.

При проведении данного исследования исходная сетка, имеющая размерность 50 узлов вдоль пластины, была дважды последовательно измельчена. Измельченные сетки имели 100 и 200 узлов соответственно. Результаты исследования сеточной сходимости полученного решения, а также сходимости этого решения к аналитическому решению по мере измельчения сетки, представлены на рисунке 6.

Рисунок 6. Сеточная сходимость

Видно, что решения, полученные на сетках, размерностью 100 и 200 узлов в продольном направлении, слабо отличаются друг от друга. Далее будем использовать сетку, размерностью 100 узлов вдоль пластины. 

Кроме решения  дифференциального уравнения (8), решение однонаправленной механической задачи проводилось в программном комплексе ANSYS Mechanical. Была создана конечно-элементная модель пластины с использованием трехмерного оболочечного элемента SHELL93 с толщиной . Количество узлов на пластине N_x=200. На границе x=0 ставится условие жёсткой заделки. К нижней и верхней стороне пластины прикладывается давление, которое было получено при решении задачи о течении в каналах (ANSYS CFX).

2.3. Двунаправленная задача

 

Параметры задачи те же, что были описаны выше.

Для реализации двунаправленного расчета используется взаимодействие между CFX и ANSYS Mechanical. Для этого необходимо выбрать опцию ANSYS MultiField via Prep7 и Mesh Deformation для расчётной области в CFX и опцию MultiField в ANSYS Mechanical.

На рис.7 представлены результаты расчета прогиба пластинки, полученные при расчете однонаправленной (CFX), двунаправленной задач (Mechanical+CFX), численного решения (Fortran) и аналитического решения.