Skip to content
Snippets Groups Projects
Commit 7db786f3 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

up

parent 30d765d1
No related branches found
No related tags found
No related merge requests found
Pipeline #
Hide whitespace changes
Inline Side-by-side
Elasticity/wrench2D.pro
+ 5
− 5
View file @ 7db786f3
......@@ -2,7 +2,7 @@
Tutorial 3 : linear elastic model of a wrench
Features:
- "Grad u" GetDP specific formulation for linear elasticity
- "grad u" GetDP specific formulation for linear elasticity
- first and second order elements
- triangular and quadrangular elements
......@@ -16,15 +16,15 @@
GetDP has a peculiar way to deal with linear elasticity. Instead of a vector
field "u = Vector[ ux, uy, uz ]", the displacement field is regarded as two
(2D case) or 3 (3D case) scalar fields. Unlike conventional formulations
then, GetDP's formulation is written in terms of the gradient "Grad u" of the
then, GetDP's formulation is written in terms of the gradient "grad u" of the
displacement field, which is a non-symmetric tensor, and the needed
symmetrization (to define the strain tensor and relate it to the stress
tensor) is done through the constitutive relationship (Hooke law). The
reason for this unusual formulation is to be able to use also for elastic
problems the powerful geometrical and homological kernel of GetDP, which
relies on the operators Grad, Curl and Div.
relies on the operators grad, curl and div.
The "Grad u" formulation entails a small increase of assembly work but makes
The "grad u" formulation entails a small increase of assembly work but makes
in counterpart lots of geometrical features implemented in GetDP (change of
coordinates, function spaces, etc...) applicable to elastic problems
out-of-the-box, since the scalar fields { ux, uy, uz } have exactly the same
......@@ -142,7 +142,7 @@ Constraint {
u . n = ux Cos [th] + uy Sin [th] = ... ;
are less naturally accounted for within the "Grad u" formulation; but they
are less naturally accounted for within the "grad u" formulation; but they
could be easily implemented with a Lagrange multplier.
Finite element shape (triangles or quadrangles) makes no difference in the
......
Electrostatics/microstrip.pro
+ 12
− 12
View file @ 7db786f3
......@@ -19,15 +19,15 @@
/* In this first tutorial we consider the calculation of the electric field
given a static distribution of electric potential. This corresponds to an
"electrostatic" physical model, obtained by combining the time-invariant
Faraday equation (Curl e = 0, with e the electric field) with Gauss'
law (Div d = rho, with d the displacement field and rho the charge density)
Faraday equation (curl e = 0, with e the electric field) with Gauss'
law (div d = rho, with d the displacement field and rho the charge density)
and the dielectric constitutive law (d = epsilon e, with epsilon the
dielectric permittivity).
Since Curl e = 0, e can be derived from a scalar electric potential v, such
that e = -Grad v. Plugging this potential in Gauss' law and using the
Since curl e = 0, e can be derived from a scalar electric potential v, such
that e = -grad v. Plugging this potential in Gauss' law and using the
constitutive law leads to a scalar Poisson equation in terms of the electric
potential: -Div(epsilon Grad v) = rho.
potential: -div(epsilon grad v) = rho.
We consider here the special case where rho = 0 to model a conducting
microstrip line on top of a dielectric substrate. A Dirichlet boundary
......@@ -52,9 +52,9 @@ Group {
/* We now define abstract regions to be used below in the definition of the
scalar electric potential formulation:
Vol_Ele : volume where -Div(epsilon Grad v) = 0 is solved
Vol_Ele : volume where -div(epsilon grad v) = 0 is solved
Sur_Neu_Ele: surface where non homogeneous Neumann boundary conditions
(on n.d = -n . (epsilon Grad v)) are imposed
(on n.d = -n . (epsilon grad v)) are imposed
Vol_xxx groups contain only volume elements of the mesh (triangles here).
Sur_xxx groups contain only surface elements of the mesh (lines here).
......@@ -172,24 +172,24 @@ Formulation {
deal with it carefully.
A GetDP Formulation encodes a so-called weak formulation of the original
partial differential equation, i.e. of -Div(epsilon Grad v) = 0. This weak
partial differential equation, i.e. of -div(epsilon grad v) = 0. This weak
formulation involves finding v such that
(-Div(epsilon Grad v) , v')_Vol_Ele = 0
(-div(epsilon grad v) , v')_Vol_Ele = 0
holds for all so-called "test-functions" v', where (.,.)_D denotes an inner
product over a domain D. If the test-functions v' are differentiable,
integration by parts using Green's identity leads to finding v such that
(epsilon Grad v, Grad v')_Vol_Ele + (n . (epsilon Grad v), v')_Bnd_Vol_Ele = 0
(epsilon grad v, grad v')_Vol_Ele + (n . (epsilon grad v), v')_Bnd_Vol_Ele = 0
holds for all v', where Bnd_Vol_Ele is the boundary of Vol_Ele. In our
microstrip example this surface term vanishes, as there is either no test
function v' (on the Dirichlet boundary), or "n. (epsilon Grad v)" is zero
function v' (on the Dirichlet boundary), or "n. (epsilon grad v)" is zero
(on the homogeneous Neumann boundary). We are thus eventually looking for
functions v in the function space Hgrad_v_Ele such that
(epsilon Grad v, Grad v')_Vol_Ele = 0
(epsilon grad v, grad v')_Vol_Ele = 0
holds for all v'. Finally, our choice here is to use a Gakerkin method,
where the test functions v' are the same basis functions ("sn_k") as the
......
ElectrostaticsFloating/floating.pro
+ 6
− 6
View file @ 7db786f3
......@@ -56,9 +56,9 @@ Group {
/* Abstract regions:
Vol_Ele : volume where -Div(epsilon Grad v) = 0 is solved
Vol_Ele : volume where -div(epsilon grad v) = 0 is solved
Sur_Neu_Ele : surface where non homogeneous Neumann boundary conditions
(on n.d = -n . (epsilon Grad v)) are imposed
(on n.d = -n . (epsilon grad v)) are imposed
Sur_Electrodes_Ele : electrode regions */
Vol_Ele = Region[ {Air, Diel1} ];
......@@ -140,15 +140,15 @@ FunctionSpace {
carried by that electrode. Indeed, let us consider the electrostatic weak
formulation derived in Tutorial 1: find v in Hgradv_Ele such that
(epsilon Grad v, Grad v')_Vol_Ele + (n . (epsilon Grad v), v')_Bnd_Vol_Ele = 0
(epsilon grad v, grad v')_Vol_Ele + (n . (epsilon grad v), v')_Bnd_Vol_Ele = 0
holds for all test functions v'. When the test-function v' is BF_electrode,
the boundary term reduces to
(n . (epsilon Grad v), BF_electrode)_Sur_Electrodes_Ele.
(n . (epsilon grad v), BF_electrode)_Sur_Electrodes_Ele.
Since BF_electrode == 1 on the electrode, the boundary term is actually
simply equal to the integral of (n . epsilon Grad v) on the electrode,
simply equal to the integral of (n . epsilon grad v) on the electrode,
i.e. the flux of the displacement field, which is by definition the
charge Q_electrode carried by the electrodes.
......@@ -215,7 +215,7 @@ Formulation {
satisfies (just consider the equation corresponding to the test function
BF_electrode):
Q_electrode = (-epsilon Grad v, Grad BF_electrode)_Vol_Ele */
Q_electrode = (-epsilon grad v, grad BF_electrode)_Vol_Ele */
{ Name Electrostatics_v; Type FemEquation;
Quantity {
{ Name v; Type Local; NameOfSpace Hgrad_v_Ele; }
......
Thermics/brick.pro
+ 113
− 145
View file @ 7db786f3
......@@ -15,149 +15,126 @@
Run (button at the bottom of the left panel)
------------------------------------------------------------------- */
/* This model is a rectangular brick with two windows,
where various kinds of thermal constraints can be set.
Dirichlet, Neumann and convection boundary conditions are imposed
on different parts of the surface of the brick.
The model is rather academic but it
demonstrates some useful high-level GetDP features.
/* This model is a rectangular brick with two windows, where various kinds of
thermal constraints can be set. Dirichlet, Neumann and convection boundary
conditions are imposed on different parts of the surface of the brick. The
model is rather academic but it demonstrates some useful high-level GetDP
features.
Governing eqations are
div ( -lambda[] grad T ) = Q in Vol_Lambda_The
-lambda[] grad T . n = qn = 0 on Sur_Neumann_The
div ( -lambda grad T ) = Q in Vol_The
-lambda grad T . n = qn = 0 on Sur_Neu_The
Contact thermal resistance:
First, it is shown how to implement contact thermal resistances
with surface elements. The surface elements are associated with a
thickness and a thermal conductivity (typically much lower than that
of surrounding regions). The implementation takes advantage
of the powerful FunctionSpace definition in GetDP.
With the flag "Flag_Regularization", the contact surface
can be, for the sake of comparison, replaced
by a thin volume conducting region.
First, it is shown how to implement contact thermal resistances with surface
elements. The surface elements are associated with a thickness and a thermal
conductivity (typically much lower than that of surrounding regions). The
implementation takes advantage of the powerful FunctionSpace definition in
GetDP.
With the flag "Flag_Regularization", the contact surface can be, for the sake
of comparison, replaced by a thin volume conducting region.
Thermal "electrode":
The floating potential idea (introduced in tutorial 4)
is reconsidered here in a thermal context to represent a region
with a very large thermal conductivity where, consequently,
the temperature field is uniform (exactly like the electric potential
is uniform on an electrode).
The dual quantity of this uniform temperature "T_electrode" [K]
(which is the "associated global quantity" in GetDP language)
is the heat flux "Q_electrode" [W] injected in the electrode
by the agent that maintains the temperature equal
The floating potential idea (introduced in tutorial 4) is reconsidered here
in a thermal context to represent a region with a very large thermal
conductivity where, consequently, the temperature field is uniform (exactly
like the electric potential is uniform on an electrode). The dual quantity
of this uniform temperature "T_electrode" [K] (which is the "associated
global quantity" in GetDP language) is the heat flux "Q_electrode" [W]
injected in the electrode by the agent that maintains the temperature equal
to the prescribed value.
The value of Q_electrode is a by-product of the system resolution
provided the term
The value of Q_electrode is a by-product of the system resolution provided
the term
GlobalTerm { [-Dof{Q_electrode} , {T_electrode} ] ; In Tfloating_The ; }
is present in the "Resolution" section. This term triggers the writing
in the linear system of a supplementary equation associated
with the global basis function BF{T_electrode}.
All integrations are automatically done by getDP,
and the value of Q_electrode is obtained in postprocessing
with the PostOperation
is present in the "Resolution" section. This term triggers the writing in the
linear system of a supplementary equation associated with the global basis
function BF{T_electrode}. All integrations are automatically done by GetDP,
and the value of Q_electrode is obtained in postprocessing with the
PostOperation
Print[ Q_electrode, OnRegion Tfloating_The, ... ]
Heat flux through surfaces:
The purpose of a thermal simulation usually goes beyond
the mere calculation of a temperature distribution.
One is in general also interested in evaluating
the heat flux q(S) through some specific surface S:
q(S) = ( -lambda[] grad T . n )_S
This quantity cannot be computed from the temperature
distribution available on the surface S only.
As heat flux is related with the gradient of temperature
in the direction normal to the surface, its computation relies on
the temperature distribution in a neighborhood of the surface.
This means that volume elements in contact with the considered surface
need be involved in the computation.
To achieve this with getDP, a good method proceeds
by the definition a smooth auxiliary function g(S),
with g(S)=1 on S, and g(S)=0 outside a finite neighborhood of S.
Typically, g(S) is the sum of the shape functions of the nodes on S.
Let w(S) be the support of g(S),
and let dw(S) denote the boundary of w(S).
We then have, just adding and substracting dw(S)
to the surface of integration S
q(S) = ( -lambda[] grad T . n g(S) )_{ dw(S) - ( dw(S)-S ) }.
dw(S) being a boundary, Stokes theorem can be invoked and,
after an integration by part one ends up with
q(S) = ( -lambda[] grad T . grad g(S) )_w(S)
- ( -lambda[] grad T . n g(S) )_{dw(S)-S}.
The purpose of a thermal simulation usually goes beyond the mere calculation
of a temperature distribution. One is in general also interested in
evaluating the heat flux q(S) through some specific surface S:
q(S) = ( -lambda grad T . n )_S
This quantity cannot be computed from the temperature distribution available
on the surface S only. As heat flux is related with the gradient of
temperature in the direction normal to the surface, its computation relies on
the temperature distribution in a neighborhood of the surface. This means
that volume elements in contact with the considered surface need be involved
in the computation. To achieve this with getDP, a good method proceeds by
the definition a smooth auxiliary function g(S), with g(S)=1 on S, and g(S)=0
outside a finite neighborhood of S. Typically, g(S) is the sum of the shape
functions of the nodes on S. Let w(S) be the support of g(S), and let dw(S)
denote the boundary of w(S). We then have, just adding and substracting
dw(S) to the surface of integration S
q(S) = ( -lambda grad T . n g(S) )_{ dw(S) - ( dw(S)-S ) }.
dw(S) being a boundary, Stokes theorem can be invoked and, after an
integration by part one ends up with
q(S) = ( -lambda grad T . grad g(S) )_w(S)
- ( -lambda grad T . n g(S) )_{dw(S)-S}.
+ ( Q g(S) )_w(S)
Now, g(S) is zero on {dw(S)-S}, except maybe at some surface elements
adjacents to dS, but not in dS.
The second terme vanishes then if either S is closed,
or adjacent to a homogeneous Neumann boundary condition.
The third term also vanishes, except if a region
with a nonzero heatsource Q is in contact with the surface S.
adjacents to dS, but not in dS. The second terme vanishes then if either S
is closed, or adjacent to a homogeneous Neumann boundary condition.
So we have nearly always the following practical formula
to evaluate the heat flux across a surface S,
in terms of a well-chosen auxiliary scalar function g(S).
The third term also vanishes, except if a region with a nonzero heatsource Q
is in contact with the surface S.
q(S) = ( -lambda[] grad T . grad g(S) )_{support of g(S)}
So we have nearly always the following practical formula to evaluate the heat
flux across a surface S, in terms of a well-chosen auxiliary scalar function
g(S).
q(S) = ( -lambda grad T . grad g(S) )_{support of g(S)}
Particular cases:
For the heat flux through the boundary of a thermal electrode,
one uses g(S) = BF(T_electrode).
Note that this heat flux is equal to Q_electrode
in the stationary case.
This is the case for the heat flux through the boundary of Window2
For the heat flux through the boundary of a thermal electrode, one uses g(S)
= BF(T_electrode). Note that this heat flux is equal to Q_electrode in the
stationary case. This is the case for the heat flux through the boundary of
Window2
Auxiliary functions g(S) are also generated
for the surfaces named "Surface_i", i=1,2,3 in the model.
Note that the flux computed through Surface_3 is incorrect
because this surface is not adjacent to surfaces
with homogeneous Neumann boundary conditions.
*/
Auxiliary functions g(S) are also generated for the surfaces named
"Surface_i", i=1,2,3 in the model. Note that the flux computed through
Surface_3 is incorrect because this surface is not adjacent to surfaces with
homogeneous Neumann boundary conditions. */
Include "brick_common.pro";
QWindow1 =
DefineNumber[1e3, Name "Window1/Heat source [W]"];
/* The user is given the choice of setting either the global temperature
or the global heat flux in Window2.
Check how the variable "Flag_ConstraintWin2" is used at different
places to alter the model according to that choice
and to manage the visibility of related input and output data.*/
QWindow2 =
DefineNumber[1e3, Name "Window2/Heat source [W]",
Visible Flag_ConstraintWin2];
TWindow2 =
DefineNumber[50, Name "Window2/Temperature [degC]",
Visible !Flag_ConstraintWin2];
outQWindow2 =
DefineNumber[0, Name "Output/Q window 2 [degC]",
Visible !Flag_ConstraintWin2, Highlight "Ivory"];
outTWindow2 =
DefineNumber[0, Name "Output/T window 2 [degC]",
Visible Flag_ConstraintWin2, Highlight "Ivory"];
ConvectionCoef =
DefineNumber[1000, Name"Surface2/hconv",
Label "Convection coefficient [W/(m^2K)]"];
T_Ambiance =
DefineNumber[20, Name"Surface2/Ambiance temperature [degC]"];
T_Dirichlet =
DefineNumber[20, Name"Surface1/Imposed temperature [degC]"];
QWindow1 = DefineNumber[1e3, Name "Window1/Heat source [W]"];
/* The user is given the choice of setting either the global temperature or the
global heat flux in Window2. Check how the variable "Flag_ConstraintWin2" is
used at different places to alter the model according to that choice and to
manage the visibility of related input and output data.*/
QWindow2 = DefineNumber[1e3, Name "Window2/Heat source [W]",
Visible Flag_ConstraintWin2];
TWindow2 = DefineNumber[50, Name "Window2/Temperature [degC]",
Visible !Flag_ConstraintWin2];
outQWindow2 = DefineNumber[0, Name "Output/Q window 2 [degC]",
Visible !Flag_ConstraintWin2, Highlight "Ivory"];
outTWindow2 = DefineNumber[0, Name "Output/T window 2 [degC]",
Visible Flag_ConstraintWin2, Highlight "Ivory"];
ConvectionCoef = DefineNumber[1000, Name"Surface2/hconv",
Label "Convection coefficient [W/(m^2K)]"];
T_Ambiance = DefineNumber[20, Name"Surface2/Ambiance temperature [degC]"];
T_Dirichlet = DefineNumber[20, Name"Surface1/Imposed temperature [degC]"];
Group {
/* Geometrical regions: */
......@@ -175,17 +152,17 @@ Group {
/* Abstract regions:
Vol_Lambda_Ele : volume regions with a thermal conductivity lambda[]
Sur_Dirichlet_The : Dirichlet bondary condition surface
Sur_Neu_Ele : Neumann bondary condition surface
Vol_The : volume regions with a thermal conductivity
Sur_Dir_The : non-homogeneous Dirichlet boundary condition surface
Sur_Neu_The : homogeneous Neumann boundary condition surface
Sur_Convection_The : convective surface q.n = h ( T-Tinf )
Vol_Qsource_The : volume heat source regions
Tfloating_The : thermal electrodes
*/
Vol_Lambda_The = Region[ {Brick, Window1, Window2, LayerWindow1} ];
Sur_Dirichlet_The = Region[ Surface_1 ];
Sur_Neumann_The = Region[ Surface_3 ];
Vol_The = Region[ {Brick, Window1, Window2, LayerWindow1} ];
Sur_Dir_The = Region[ Surface_1 ];
Sur_Neu_The = Region[ Surface_3 ];
Sur_Convection_The = Region[ Surface_2 ];
Vol_Qsource_The = Region[ Window1 ];
......@@ -200,7 +177,6 @@ Group {
EndIf
}
Function {
lambda_Brick = 50.; // steel
lambda[ Brick ] = lambda_Brick;
......@@ -226,7 +202,7 @@ Function {
Constraint {
{ Name Dirichlet_The ;
Case {
{ Region Sur_Dirichlet_The ; Value T_Dirichlet ; }
{ Region Sur_Dir_The ; Value T_Dirichlet ; }
}
}
{ Name T_Discontinuity ;
......@@ -260,8 +236,7 @@ Constraint {
Integration {
{ Name Int ;
Case {
{
Type Gauss ;
{ Type Gauss ;
Case {
{ GeoElement Triangle ; NumberOfPoints 4 ; }
{ GeoElement Quadrangle ; NumberOfPoints 4 ; }
......@@ -286,12 +261,9 @@ Jacobian {
}
Group {
Dom_Hgrad_T = Region[ {Vol_Lambda_The,
Sur_Convection_The,
Sur_Tdisc_The} ];
Dom_Hgrad_T = Region[ {Vol_The, Sur_Convection_The, Sur_Tdisc_The} ];
DomainWithSurf_TL_The =
ElementsOf[ {Vol_OneSide_The, Sur_Tdisc_The},
OnOneSideOf Sur_Tdisc_The ];
ElementsOf[ {Vol_OneSide_The, Sur_Tdisc_The}, OnOneSideOf Sur_Tdisc_The ];
}
FunctionSpace {
......@@ -336,7 +308,6 @@ FunctionSpace {
}
}
EndFor
}
Formulation {
......@@ -350,36 +321,33 @@ Formulation {
For i In {1:NbSurface}
{ Name un~{i} ; Type Local ; NameOfSpace FluxLayer~{i} ; }
EndFor
}
Equation {
Integral { [ lambda[] * Dof{d T} , {d T} ] ;
In Vol_Lambda_The; Integration Int ; Jacobian Vol ; }
In Vol_The; Integration Int ; Jacobian Vol ; }
Integral { [ ( lambda[]/thickness[] ) * Dof{Tdisc} , {Tdisc} ] ;
In Sur_Tdisc_The; Integration Int ; Jacobian Sur ; }
In Sur_Tdisc_The; Integration Int ; Jacobian Sur ; }
Integral { [ -Qsource[] , {T} ] ;
In Vol_Qsource_The ; Integration Int ; Jacobian Vol ; }
In Vol_Qsource_The ; Integration Int ; Jacobian Vol ; }
Integral { [ h[] * Dof{T} , {T} ] ;
In Sur_Convection_The ; Integration Int ; Jacobian Sur ; }
In Sur_Convection_The ; Integration Int ; Jacobian Sur ; }
Integral { [ -h[] * Tinf[] , {T} ] ;
In Sur_Convection_The ; Integration Int ; Jacobian Sur ; }
In Sur_Convection_The ; Integration Int ; Jacobian Sur ; }
GlobalTerm { [-Dof{Qglob} , {Tglob} ] ; In Tfloating_The ; }
For i In {1:NbSurface}
Integral { [ 0 * Dof{un~{i}} , {un~{i}} ] ;
In Vol_Lambda_The ; Integration Int ; Jacobian Vol ; }
In Vol_The ; Integration Int ; Jacobian Vol ; }
EndFor
}
}
}
Resolution {
{ Name Thermal_T ;
System {
......@@ -411,9 +379,9 @@ PostProcessing {
For i In {1:NbSurface}
{ Name un~{i} ; Value { Local { [ {un~{i}} ] ;
In Vol_Lambda_The ; Jacobian Vol ; } } }
In Vol_The ; Jacobian Vol ; } } }
{ Name IFlux~{i} ; Value { Integral { [ -lambda[]*{d T} * {d un~{i}} ];
In Vol_Lambda_The ; Jacobian Vol ; Integration Int ; } } }
In Vol_The ; Jacobian Vol ; Integration Int ; } } }
EndFor
}
......@@ -422,8 +390,8 @@ PostProcessing {
PostOperation Map_T UsingPost Thermal_T {
If( !Flag_Regularization )
Print[ Tcont, OnElementsOf Vol_Lambda_The, File "Tcont.pos"] ;
Print[ Tdisc, OnElementsOf Vol_Lambda_The, File "Tdisc.pos"] ;
Print[ Tcont, OnElementsOf Vol_The, File "Tcont.pos"] ;
Print[ Tdisc, OnElementsOf Vol_The, File "Tdisc.pos"] ;
EndIf
If(Flag_ConstraintWin2)
......@@ -435,13 +403,13 @@ PostOperation Map_T UsingPost Thermal_T {
EndIf
For i In {1:NbSurface}
Print[un~{i}, OnElementsOf Vol_Lambda_The,
Print[un~{i}, OnElementsOf Vol_The,
File Sprintf("FluxLayer_%g.pos",i)];
Print[ IFlux~{i}[Vol_Lambda_The], OnGlobal,
Print[ IFlux~{i}[Vol_The], OnGlobal,
Format TimeTable, File > "Fluxes.dat", Color "Ivory",
SendToServer Sprintf("Output/Heat flux surface %g [W]", i)];
EndFor
Print[ T, OnElementsOf Vol_Lambda_The, File "T.pos" ] ;
Print[ T, OnElementsOf Vol_The, File "T.pos" ] ;
Echo[ StrCat["l=PostProcessing.NbViews-1;",
"View[l].IntervalsType = 3;",
"View[l].NbIso = 30;",
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment