|
00001 // Copyright (C) 2004 Johan Hoffman. 00002 // Licensed under the GNU GPL Version 2. 00003 00004 #ifndef __NSE_H 00005 #define __NSE_H 00006 00007 #include <dolfin/PDE.h> 00008 00009 namespace dolfin { 00010 00011 class NSE_Momentum : public PDE { 00012 public: 00013 00014 NSE_Momentum(Function::Vector& source, 00015 Function::Vector& uprevious, 00016 Function::Vector& convection, 00017 Function::Vector& pressure) : 00018 PDE(3, 3), f(1), up(3), b(3), p(1) 00019 { 00020 add(f, source); 00021 add(up, uprevious); 00022 add(b, convection); 00023 add(p, pressure); 00024 00025 C1 = 2.0; 00026 C2 = 1.0; 00027 00028 nu = 1.0/1000.0; 00029 00030 bcwt = 1.0e7; 00031 } 00032 00033 real lhs(ShapeFunction::Vector& u, ShapeFunction::Vector& v) 00034 { 00035 00036 unorm = sqrt(sqr(b(0)(cell_->node(0).id()))+sqr(b(1)(cell_->node(0).id()))+sqr(b(2)(cell_->node(0).id()))); 00037 if ( (h/nu) > 1.0 ) d1 = C1 * (0.5 / sqrt( 1.0/sqr(k) + sqr(unorm/h) )); 00038 else d1 = C1 * sqr(h); 00039 00040 if ( (h/nu) > 1.0 ) d2 = C2 * h; 00041 else d2 = C2 * sqr(h); 00042 00043 MASS = (u(0)*v(0) + u(1)*v(1) + u(2)*v(2)) * dx; 00044 STIFF = nu*((grad(u(0)),grad(v(0))) + (grad(u(1)),grad(v(1))) + (grad(u(2)),grad(v(2)))) * dx; 00045 NL = ((b,grad(u(0)))*v(0) + (b,grad(u(1)))*v(1) + (b,grad(u(2)))*v(2)) * dx; 00046 00047 SD = (d1*((b,grad(u(0)))*(b,grad(v(0))) + (b,grad(u(1)))*(b,grad(v(1))) + (b,grad(u(2)))*(b,grad(v(2)))) + 00048 d2*(ddx(u(0))+ddy(u(1))+ddz(u(2)))*(ddx(v(0))+ddy(v(1))+ddz(v(2)))) * dx; 00049 00050 BC = bcwt*(u(0)*v(0) + u(1)*v(1) + u(2)*v(2)) * ds; 00051 00052 return (MASS*(1.0/k) + 0.5*(STIFF + NL + SD) + BC); 00053 00054 } 00055 00056 real rhs(ShapeFunction::Vector& v) 00057 { 00058 00059 unorm = sqrt(sqr(b(0)(cell_->node(0).id()))+sqr(b(1)(cell_->node(0).id()))+sqr(b(2)(cell_->node(0).id()))); 00060 if ( (h/nu) > 1.0 ) d1 = C1 * (0.5 / sqrt( 1.0/sqr(k) + sqr(unorm/h) )); 00061 else d1 = C1 * sqr(h); 00062 00063 if ( (h/nu) > 1.0 ) d2 = C2 * h; 00064 else d2 = C2 * sqr(h); 00065 00066 MASS = ((up(0),v(0)) + (up(1),v(1)) + (up(2),v(2))) * dx; 00067 STIFF = nu*( ddx(up(0))*ddx(v(0)) + ddy(up(0))*ddy(v(0)) + ddz(up(0))*ddz(v(0)) + 00068 ddx(up(1))*ddx(v(1)) + ddy(up(1))*ddy(v(1)) + ddz(up(1))*ddz(v(1)) + 00069 ddx(up(2))*ddx(v(2)) + ddy(up(2))*ddy(v(2)) + ddz(up(2))*ddz(v(2)) ) * dx; 00070 NL = ( (b(0)*ddx(up(0)) + b(1)*ddy(up(0)) + b(2)*ddz(up(0)))*v(0) + 00071 (b(0)*ddx(up(1)) + b(1)*ddy(up(1)) + b(2)*ddz(up(1)))*v(1) + 00072 (b(0)*ddx(up(2)) + b(1)*ddy(up(2)) + b(2)*ddz(up(2)))*v(2) ) * dx; 00073 SD = ( d1*((b(0)*ddx(up(0)) + b(1)*ddy(up(0)) + b(2)*ddz(up(0)))*(b,grad(v(0))) + 00074 (b(0)*ddx(up(1)) + b(1)*ddy(up(1)) + b(2)*ddz(up(1)))*(b,grad(v(0))) + 00075 (b(0)*ddx(up(2)) + b(1)*ddy(up(2)) + b(2)*ddz(up(2)))*(b,grad(v(0)))) + 00076 d2*(ddx(up(0))+ddy(up(1))+ddz(up(2)))*(ddx(v(0))+ddy(v(1))+ddz(v(2))) ) * dx; 00077 SDP = d1*(ddx(p(0))*(b,grad(v(0))) + ddy(p(0))*(b,grad(v(1))) + ddz(p(0))*(b,grad(v(2)))) * dx; 00078 SDF = d1*(f(0)*(b,grad(v(0))) + f(1)*(b,grad(v(1))) + f(2)*(b,grad(v(2)))) * dx; 00079 PDV = p(0)*(ddx(v(0))+ddy(v(1))+ddz(v(2))) * dx; 00080 F = ((f(0),v(0)) + (f(1),v(1)) + (f(2),v(2))) * dx; 00081 00082 ud0 = 0.0; ud1 = 0.0; ud2 = 0.0; 00083 un0 = 0.0; un1 = 0.0; un2 = 0.0; 00084 BC = bcwt*((ud0-un0)*v(0) + (ud1-un1)*v(1) + (ud2-un2)*v(2)) * ds; 00085 00086 return (MASS*(1.0/k) - 0.5*(STIFF + NL + SD) - SDP + SDF + PDV + F + BC); 00087 00088 } 00089 00090 private: 00091 ElementFunction::Vector f; // Source term 00092 ElementFunction::Vector up; // Velocity value at left end-point 00093 ElementFunction::Vector b; // Convection = linearized velocity 00094 ElementFunction::Vector p; // linearized pressure 00095 00096 real nu,d1,d2,C1,C2,unorm,SD,SDP,SDF,MASS,STIFF,NL,PDV,F,BC,ud0,ud1,ud2,un0,un1,un2,bcwt; 00097 }; 00098 00099 //------------------------------------------------------------------------------------------- 00100 00101 class NSE_Continuity : public PDE { 00102 public: 00103 00104 NSE_Continuity(Function::Vector& source, 00105 Function::Vector& convection) : 00106 PDE(1, 1), f(1), b(3) 00107 { 00108 add(f, source); 00109 add(b, convection); 00110 00111 C1 = 2.0; 00112 00113 nu = 1.0/1000.0; 00114 00115 bcwt = 0.0; 00116 } 00117 00118 real lhs(ShapeFunction& u, ShapeFunction& v) 00119 { 00120 00121 STIFF = (grad(u),grad(v)) * dx; 00122 00123 BC = bcwt*(u*v) * ds; 00124 00125 return (STIFF + BC); 00126 00127 } 00128 00129 real rhs(ShapeFunction& v) 00130 { 00131 00132 unorm = sqrt(sqr(b(0)(cell_->node(0).id()))+sqr(b(1)(cell_->node(0).id()))+sqr(b(2)(cell_->node(0).id()))); 00133 if ( (h/nu) > 1.0 ) d1 = C1 * (0.5 / sqrt( 1.0/sqr(k) + sqr(unorm/h) )); 00134 else d1 = C1 * sqr(h); 00135 00136 SD = ((b(0)*ddx(b(0)) + b(1)*ddy(b(0)) + b(2)*ddz(b(0)))*v.ddx() + 00137 (b(0)*ddx(b(1)) + b(1)*ddy(b(1)) + b(2)*ddz(b(1)))*v.ddy() + 00138 (b(0)*ddx(b(2)) + b(1)*ddy(b(2)) + b(2)*ddz(b(2)))*v.ddz()) * dx; 00139 00140 DIV = (ddx(b(0)) + ddy(b(1)) + ddz(b(2)))*v * dx; 00141 00142 ud = 0.0; 00143 un = 0.0; 00144 BC = bcwt*((ud-un)*v) * ds; 00145 00146 return ( (-1.0)*(SD + (1.0/d1)*DIV) + BC); 00147 00148 } 00149 00150 private: 00151 00152 ElementFunction::Vector f; // Source term 00153 ElementFunction::Vector b; // Convection = linearized velocity 00154 00155 real nu,d1,d2,C1,C2,unorm,bcwt,STIFF,SD,DIV,BC,ud,un; 00156 00157 }; 00158 00159 00160 00161 00162 } 00163 00164 #endif
Documentation automatically generated with Doxygen on 10 Sep 2004.