[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[ff3d-users] Problem with elliptic hollow cylinder section
From: |
Tobias Nähring |
Subject: |
[ff3d-users] Problem with elliptic hollow cylinder section |
Date: |
Thu, 22 Apr 2004 21:44:13 +0200 |
User-agent: |
Mozilla/5.0 (X11; U; Linux i686; en-US; rv:1.6) Gecko/20040401 Debian/1.6-4 |
Note:
A short summary of the problem with the corresponding pictures and
files is available at: http://www.tn-home.de/Soft/Fem3d/index.html.
I want to solve the Laplacian equation in the domain
G := {(x,y,z)| 1 < (x/rx)^2 + (y/ry)^2 < a^2; 0 < z < L }
which resembles a section of a hollow cylinder with elliptical
cross-section. The constants are: rx=1, ry=0.5, a=2, L=5.
The boundary condition is
u(x,y,z) = (|x,y|-1)*(a-|x,y|)*(L-z)/L
with
|x,y|:=sqrt( (x/rx)^2 + (y/ry)^2 ).
That means effectively that u(x,y,z)=0 for z=L, and for (x,y,z) on the
cylinder mantles. Only at z=0 one has some nice non-zero boundary
condition.
Because of the symmetry of the boundary condition the
computational domain can be confined to the quarter
Q:={(x,y,z) in G| x>0, y>0}
of G with natural boundary conditions for u on the boundaries at x=0
and y=0.
I tried to solve this problem with FreeFEM3D without much luck.
The results of my attempts are shown at
http://www.tn-home.de/Tobias/Soft/Fem3d/index.html
ff3d accepts the input file and the corresponding pov file (included below)
and the solver options are set to the save values
cg(epsilon=1e-8) and cg(maxiter=2000)
but the output is rather erroneous. I have played
already with other solver settings, such as bicg(...) and
krylov(type=bicg) without any luck. The grid size is as large as
possible on my computer with 512MB RAM. (If I enlarge the grid size
then ff3d freezes during "Marching cube ...".
With the current grid size I get the error message:
.
warning : meshing object plane {
(0, 1, 0),2
}
a problem occured ...
Check that:
- the object is not to small compare to the background mesh
- the object intersects the background mesh
I don't know how to handle that.
The ff3d computed solution of the problem are graphically represented
in the figures:
http://www.tn-home.de/Tobias/Soft/Fem3d/
Is there any chance to get a better solution?
One more remark: There would be a way to specify the mantle surface of
a cylinder section: the `open' directive of the povray object
cylinder (as described in the povray manual).
The input file for ff3d:
//////////////////////////////////////////////////////////////////////
// Model of a cable with ellipsoidal cross-section.
vector Grid=(70,60,30);
double rx = 1.0; // radius of the inner conductor in x-direction
double ry = 0.5; // radius of the inner conductor in y-direction
double a=2.0; // The cross-section boundary of the inner conductor
// scaled by a gives the cross-section boundary of the
// shielding.
double CableW = a*rx;
double CableH = a*ry;
double CableL = 5.0; // length of the cable
double Border = 0.1; // Border added around the modelled piece of the
// cable to build a fictitious domain
scene myScene = pov("cable-ellipse-annulus-num.pov");
mesh myMesh = structured( Grid, (-Border,-Border,-2.0*Border),
(a*rx+Border,a*ry+Border,CableL+2.0*Border));
// Comment: in z-direcion we provided a bit more space for the
fictitious domain (just guessing).
domain myDomain = domain(myScene, inside(<1,0,0>));
mesh myCylMesh = surface(myDomain,myMesh);
function norm2d=sqrt((x/rx)^2+(y/ry)^2);
// Realizing the boundary conditions with method-type penalty doesnot
// work (no convergence).
solve(u) in myDomain by myMesh
cg(epsilon=1e-8),
cg(maxiter=2000)
{
pde(u) div(grad(u))=0;
u=(norm2d-1)*(a-norm2d)*(CableL-z)/CableL on myCylMesh;
}
//////////////////////////////////////////////////////////////////////
// Results output for postprocessing via gnuplot:
double myx; double myy; double myz;
double mydx = CableW/x(Grid);
double mydy = CableH/y(Grid);
double mydz = CableL/z(Grid);
double myxmax=CableW+0.5*mydx;
double myymax=CableH+0.5*mydy;
double myzmax=CableL+0.5*mydz;
for( myz=0.0; myz <= myzmax; myz = myz + mydz ) {
for( myx=0.0; myx <= myxmax; myx = myx + mydx ) {
for( myy=0.0; myy <= myymax; myy = myy + mydy ) {
cout << "@" << myx << ", " << myy << ", " << myz << ", " <<
u(myx,myy,myz) << "\n";
}
cout << "@\n";
}
cout << "@\n";
}
//////////////////////////////////////////////////////////////////////
The pov file (named cable-ellipse-annulus-num.pov):
//////////////////////////////////////////////////////////////////////
// Model of a cable with ellipsoidal cross-section.
// rx=1.0 radius of the inner conductor in x-direction
// ry=0.5 radius of the inner conductor in y-direction
// a=2.0 the inner conductor scaled by a gives the shielding
// CableL=5 length of the cable.
// Since we impose symmetric boundary conditions we need only
// to model a quarter of the cable. (0 < x < a*rx, 0 < y < a*ry)
intersection {
difference {
cylinder {
<0,0,0>,
<0,0,1>,
2.0 // a==2.0
}
cylinder{
<0,0,-0.01>,
<0,0,1.01>,
1.0
}
plane{
<0,1,0>, 0
}
plane{
<1,0,0>, 0
}
}
box{
<0,0,0>,
<2,2,1>
}
scale <1.0,0.5,5.0> // <rx,ry,CableL>
pigment { color rgb <1,0,0> }
}
- [ff3d-users] Problem with elliptic hollow cylinder section,
Tobias Nähring <=