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

Merge branch 'master' of https://gitlab.onelab.info/doc/tutorials

parents f6a03d0d c9549fb0
Branches
No related tags found
No related merge requests found
Pipeline #2015 passed
...@@ -2,7 +2,7 @@ ...@@ -2,7 +2,7 @@
Include "R_circuit_common.pro"; Include "R_circuit_common.pro";
lc = 0.25; // Cjaracteristic length lc = 0.25; // Characteristic length
Point(1) = {0, 0, 0, lc}; Point(1) = {0, 0, 0, lc};
Extrude {1, 0, 0} { Extrude {1, 0, 0} {
......
/* ------------------------------------------------------------------- /* -------------------------------------------------------------------
Tutorial 3 : linear elastic model of a wrench Tutorial 3 : linear elastic model of a wrench
Features: Features:
- "grad u" GetDP specific formulation for linear elasticity - "grad u" GetDP specific formulation for linear elasticity
- first and second order elements - first and second order elements
- triangular and quadrangular elements - triangular and quadrangular elements
To compute the solution interactively from the Gmsh GUI: To compute the solution in a terminal:
getdp wrench2D.pro -solve Elast_u -pos pos
To compute the solution interactively from the Gmsh GUI:
File > Open > wrench.pro File > Open > wrench.pro
Run (button at the bottom of the left panel) Run (button at the bottom of the left panel)
------------------------------------------------------------------- */ ------------------------------------------------------------------- */
/* Linear elasticity with GetDP: /*
Particularities of linear elasticity in GetDP:
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 Instead of a vector field "u = Vector[ ux, uy, uz ]", the displacement field
(2D case) or 3 (3D case) scalar fields. Unlike conventional formulations is regarded as two (2D case) or three (3D case) scalar fields.
then, GetDP's formulation is written in terms of the gradient "grad u" of the Unlike conventional formulations, GetDP formulation is then written in terms
displacement field, which is a non-symmetric tensor, and the needed of the gradient "grad u" of the displacement field, which is a non-symmetric
symmetrization (to define the strain tensor and relate it to the stress tensor, and the needed symmetrization (to define the strain tensor and relate
tensor) is done through the constitutive relationship (Hooke law). The it to the stress tensor) is done via the constitutive relationship (Hooke law).
reason for this unusual formulation is to be able to use also for elastic This unusual formulation allows to take advantage of the powerful geometrical
problems the powerful geometrical and homological kernel of GetDP, which and homological GetDP kernel, 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
in counterpart lots of geometrical features implemented in GetDP (change of coordinates, function spaces, etc...) applicable to elastic problems
coordinates, function spaces, etc...) applicable to elastic problems out-of-the-box, since the scalar fields { ux, uy, uz } have exactly the same
out-of-the-box, since the scalar fields { ux, uy, uz } have exactly the same geometrical properties as, e.g. an electric scalar potential or a temperature
geometrical properties as, e.g. a scalar electric potential or a temperature field.
field. */ */
Include "wrench2D_common.pro"; Include "wrench2D_common.pro";
...@@ -54,24 +57,30 @@ Group { ...@@ -54,24 +57,30 @@ Group {
Sur_Clamp_Mec = Region[ Grip ]; Sur_Clamp_Mec = Region[ Grip ];
Sur_Force_Mec = Region[ Force ]; Sur_Force_Mec = Region[ Force ];
/* Signification of the abstract regions: /* Meaning of abstract regions:
- Vol_Mec : Elastic domain - Vol_Mec : Elastic domain
- Vol_Force_Mec : Region with imposed volumic force - Vol_Force_Mec : Region with imposed volume force
- Sur_Force_Mec : Surface with imposed surface traction - Sur_Force_Mec : Surface with imposed surface traction
- Sur_Clamp_Mec : Surface with imposed zero displacements (all components) */ - Sur_Clamp_Mec : Surface with imposed zero displacements (all components) */
} }
Function { Function {
/* Material coefficients. /*
No need to define them regionwise here ( E[{Wrench}] = ... ; ) Material coefficients
as there is only one region in this model. */ No need of regionwise definition ( E[{Wrench}] = ... ; )
as this model comprises only one region.
If there is more than one region and coefficients are not particularised,
the same value holds for all of them.
*/
E[] = Young; E[] = Young;
nu[] = Poisson; nu[] = Poisson;
/* Components of the volumic force applied to the region "Vol_Force_Mec" /*
Gravity could be defined here ( force_y[] = 7000*9.81; ) ; */ Volume force components applied to the region "Vol_Force_Mec"
Gravity could be defined here as well ( force_y[] = 7000*9.81; ) ;
*/
force_x[] = 0; force_x[] = 0;
force_y[] = 0; force_y[] = 0;
/* Components of the surface traction force applied to the region "Sur_Force_Mec" */ /* Surface traction force components applied to the region "Sur_Force_Mec" */
pressure_x[] = 0; pressure_x[] = 0;
pressure_y[] = -AppliedForce/(SurfaceArea[]*Thickness); // downward vertical force pressure_y[] = -AppliedForce/(SurfaceArea[]*Thickness); // downward vertical force
} }
...@@ -82,7 +91,7 @@ Function { ...@@ -82,7 +91,7 @@ Function {
sigma_ij = C_ijkl epsilon_ij sigma_ij = C_ijkl epsilon_ij
is represented in 2D by 4 2x2 tensors C_ij[], i,j=1,2, depending on the Lamé is represented in 2D by four 2x2 tensors C_ij[], i,j=1,2, depending on the Lamé
coefficients of the isotropic linear material, coefficients of the isotropic linear material,
lambda = E[]*nu[]/(1.+nu[])/(1.-2.*nu[]) lambda = E[]*nu[]/(1.+nu[])/(1.-2.*nu[])
...@@ -94,7 +103,8 @@ Function { ...@@ -94,7 +103,8 @@ Function {
EPD: a[] = lambda + 2 mu b[] = mu c[] = lambda EPD: a[] = lambda + 2 mu b[] = mu c[] = lambda
3D: a[] = lambda + 2 mu b[] = mu c[] = lambda 3D: a[] = lambda + 2 mu b[] = mu c[] = lambda
respectively for the 2D plane strain (EPD), 2D plane stress (EPS) and 3D cases. */ respectively for the 2D plane strain (EPD), 2D plane stress (EPS) and 3D cases.
*/
Function { Function {
Flag_EPC = 1; Flag_EPC = 1;
...@@ -128,7 +138,7 @@ Constraint { ...@@ -128,7 +138,7 @@ Constraint {
} }
} }
/* As explained above, the displacement field is discretized as two scalar /* As mentioned above, the displacement field is discretized as two scalar
fields "ux" and "uy", which are the spatial components of the vector field fields "ux" and "uy", which are the spatial components of the vector field
"u" in a fixed Cartesian coordinate system. "u" in a fixed Cartesian coordinate system.
...@@ -143,17 +153,17 @@ Constraint { ...@@ -143,17 +153,17 @@ Constraint {
u . n = ux Cos [th] + uy Sin [th] = ... ; 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. could be easily implemented with e.g. a Lagrange multiplier.
Finite element shape (triangles or quadrangles) makes no difference in the The finite element shape (e.g. triangles or quadrangles in 2D) has no influence
definition of the FunctionSpaces. The appropriate shape functions to be used in the definition of the FunctionSpaces. The appropriate shape functions
are determined by GetDP at a much lower level on basis of the information are determined by GetDP at a much lower level on basis of the information
contained in the *.msh file. contained in the *.msh file.
Second order elements, on the other hand, are implemented in the hierarchical Second order elements are hierarchically implemented by adding to the first
fashion by adding to the first order node-based shape functions a set of order node-based shape functions a set of second order edge-based functions
second order edge-based functions to complete a basis for 2d order to complete a basis for 2D order polynomials on the reference element.
polynomials on the reference element. */ */
// Domain of definition of the "ux" and "uy" FunctionSpaces // Domain of definition of the "ux" and "uy" FunctionSpaces
Group { Group {
...@@ -217,8 +227,8 @@ Jacobian { ...@@ -217,8 +227,8 @@ Jacobian {
} }
} }
/* Adapt the number of Gauss points to the polynomial degree of the finite elements /* Adapting the number of Gauss points to the polynomial degree of the finite elements
is as simple as this: */ is simple: */
Integration { Integration {
{ Name Gauss_v; { Name Gauss_v;
Case { Case {
......
...@@ -36,12 +36,12 @@ ...@@ -36,12 +36,12 @@
boundary condition (zero flux of the displacement field, i.e. n.d = 0) is boundary condition (zero flux of the displacement field, i.e. n.d = 0) is
imposed on the left boundary of the domain to account for the symmetry of the imposed on the left boundary of the domain to account for the symmetry of the
problem, as well as on the top and right boundaries that truncate the problem, as well as on the top and right boundaries that truncate the
modelling domain. */ simulation domain. */
Group { Group {
/* One starts by giving explicit meaningful names to the Physical regions /* One starts by giving explicit meaningful names to the Physical regions
defined in the "microstrip.msh" mesh file. We only use 2 volume regions and defined in the "microstrip.msh" mesh file. This model comprises only
2 surface regions in this model. */ 2 volume regions and 2 surface regions. */
Air = Region[101]; Air = Region[101];
Diel1 = Region[111]; Diel1 = Region[111];
...@@ -61,6 +61,10 @@ Group { ...@@ -61,6 +61,10 @@ Group {
Since there are no non-homogeneous Neumann conditions in this particular Since there are no non-homogeneous Neumann conditions in this particular
example, Sur_Neu_Ele is defined as empty. example, Sur_Neu_Ele is defined as empty.
Note that volume elements are those that correspond to the higher dimension
of the model at hand (2D elements here), surface elements correspond to the
higher dimension of the model minus one (1D elements here).
*/ */
Vol_Ele = Region[ {Air, Diel1} ]; Vol_Ele = Region[ {Air, Diel1} ];
...@@ -78,8 +82,8 @@ Function { ...@@ -78,8 +82,8 @@ Function {
} }
Constraint { Constraint {
/* As for material laws, the Dirichlet boundary condition is defined /* The Dirichlet boundary condition is also defined piecewise.
piecewise. The constraint "Dirichlet_Ele" is invoked in the FunctionSpace The constraint "Dirichlet_Ele" is invoked in the FunctionSpace
below. */ below. */
{ Name Dirichlet_Ele; Type Assign; { Name Dirichlet_Ele; Type Assign;
...@@ -97,13 +101,13 @@ Group{ ...@@ -97,13 +101,13 @@ Group{
Contrary to Vol_xxx and Sur_xxx regions, which contain only volume or Contrary to Vol_xxx and Sur_xxx regions, which contain only volume or
surface regions, resp., domains of definition Dom_xxx regions may contain surface regions, resp., domains of definition Dom_xxx regions may contain
both volume and surface regions. Hence the use of the prefixes Vol_, Sur_ both volume and surface regions. Hence the use of the prefixes Vol_, Sur_
and Dom_ to avoid confusions.*/ and Dom_ to avoid confusion.*/
Dom_Hgrad_v_Ele = Region[ {Vol_Ele, Sur_Neu_Ele} ]; Dom_Hgrad_v_Ele = Region[ {Vol_Ele, Sur_Neu_Ele} ];
} }
FunctionSpace { FunctionSpace {
/* The function space in which we shall pick the electric scalar potential "v" /* The function space in which we pick the electric scalar potential "v"
solution is defined by solution is defined by
- a domain of definition (the "Support": "Dom_Hgrad_v_Ele") - a domain of definition (the "Support": "Dom_Hgrad_v_Ele")
- a type ("Form0" means scalar field) - a type ("Form0" means scalar field)
...@@ -127,7 +131,7 @@ FunctionSpace { ...@@ -127,7 +131,7 @@ FunctionSpace {
{ Name sn; NameOfCoef vn; Function BF_Node; { Name sn; NameOfCoef vn; Function BF_Node;
Support Dom_Hgrad_v_Ele; Entity NodesOf[ All ]; } Support Dom_Hgrad_v_Ele; Entity NodesOf[ All ]; }
// using "NodesOf[All]" instead of "NodesOf[Dom_Hgrad_v_Ele]" is an // using "NodesOf[All]" instead of "NodesOf[Dom_Hgrad_v_Ele]" is an
// optimization, which allows GetDP to not explicitly build the list of // optimization, which avoids explicitly building the list of
// all the nodes // all the nodes
} }
Constraint { Constraint {
...@@ -138,7 +142,7 @@ FunctionSpace { ...@@ -138,7 +142,7 @@ FunctionSpace {
Jacobian { Jacobian {
/* Jacobians are used to specify the mapping between elements in the mesh and /* Jacobians are used to specify the mapping between elements in the mesh and
the reference elements (defined in standardized unit cells) over which the reference elements (defined in standardized/reference unit cells) over which
integration is performed. "Vol" represents the classical 1-to-1 mapping integration is performed. "Vol" represents the classical 1-to-1 mapping
between identical spatial dimensions, i.e. in this case a reference between identical spatial dimensions, i.e. in this case a reference
triangle/quadrangle onto triangles/quadrangles in the z=0 plane (2D <-> triangle/quadrangle onto triangles/quadrangles in the z=0 plane (2D <->
...@@ -173,7 +177,7 @@ Integration { ...@@ -173,7 +177,7 @@ Integration {
} }
Formulation { Formulation {
/* The syntax of the Formulation section is a harder nut to crack. So let's /* The syntax of the Formulation section is a harder nut to crack. So let us
deal with it carefully. deal with it carefully.
A GetDP Formulation encodes a so-called weak formulation of the original A GetDP Formulation encodes a so-called weak formulation of the original
...@@ -201,11 +205,11 @@ Formulation { ...@@ -201,11 +205,11 @@ Formulation {
ones used to interpolate the unknown potential v. ones used to interpolate the unknown potential v.
The "Integral" statement in the Formulation is a symbolic representation of The "Integral" statement in the Formulation is a symbolic representation of
this weak formulation. It has got 4 semicolon separated arguments: this weak formulation. It has 4 semicolon separated arguments:
* the density [.,.] to be integrated (note the square brackets instead of * the density [.,.] to be integrated (note the square brackets instead of
the parentheses), with the test-functions (always) after the comma the parentheses), with the test-functions (always) after the comma;
* the domain of integration, * the domain of integration;
* the Jacobian of the transformation reference element <-> real element, * the Jacobian of the transformation reference element <-> real element;
* the integration method to be used. * the integration method to be used.
In the density, braces around a quantity (such as "{v}") indicate that this In the density, braces around a quantity (such as "{v}") indicate that this
...@@ -223,33 +227,33 @@ Formulation { ...@@ -223,33 +227,33 @@ Formulation {
be something like this [ ... , basis_functions_of {d v} ]. However, since be something like this [ ... , basis_functions_of {d v} ]. However, since
the second term is always devoted to test functions, the operator the second term is always devoted to test functions, the operator
"basis_functions_of" would always be there. It can therefore be made "basis_functions_of" would always be there. It can therefore be made
implicit, and, according to the GetDP syntax, it is omitted. So, one writes implicit and, according to the GetDP syntax, omitted. So, one simply writes
simply [ ... , {d v} ]. [ ... , {d v} ].
The first term, on the other hand, can contain a much wider variety of On the other hand, the first term can contain a much wider variety of
expressions than the second one. In our case it should be expressed in expressions than the second one. In our case it should be expressed in
terms of the FE expansion of "v" at the present system solution, i.e. when terms of the FE expansion of "v" at the present system solution, i.e. when
the coefficients vn_k in the expansion of "v = Sum_k vn_k sn_k" are the coefficients vn_k in the expansion of "v = Sum_k vn_k sn_k" are
unknown. This is indicated by prefixing the braces with "Dof", which leads unknown. This is indicated by prefixing the braces with "Dof" (degrees of
to the following density: freedom), which leads to the following density:
[ epsilon[] * Dof{d v} , {d v} ], [ epsilon[] * Dof{d v} , {d v} ],
a so-called bilinear term that will contribute to the stiffness matrix of a so-called bilinear term that will contribute to the stiffness matrix of
the electrostatic problem at hand. the electrostatic problem at hand.
Another option, which would not work here, is to evaluate the first Another option, which does not work here, is to evaluate the first
argument with the last available already computed solution, i.e. simply argument with the last available computed solution, i.e. simply
perform the interpolation with known coefficients vn_k. For this the Dof perform the interpolation with known coefficients vn_k. The Dof
prefix would be omitted and we would have: prefix is then omitted and we have:
[ epsilon[] * {d v} , {d v} ], [ epsilon[] * {d v} , {d v} ],
a so-called linear term that will contribute to the right-hand side of the a so-called linear term that contributes to the right-hand side of the
linear system. linear system.
Both choices are commonly used in different contexts, and we shall come Both choices are commonly used in different contexts, and we shall often
back on this often in subsequent tutorials. */ come back in subsequent tutorials. */
{ Name Electrostatics_v; Type FemEquation; { Name Electrostatics_v; Type FemEquation;
Quantity { Quantity {
...@@ -259,20 +263,21 @@ Formulation { ...@@ -259,20 +263,21 @@ Formulation {
Integral { [ epsilon[] * Dof{d v} , {d v} ]; Integral { [ epsilon[] * Dof{d v} , {d v} ];
In Vol_Ele; Jacobian Vol; Integration Int; } In Vol_Ele; Jacobian Vol; Integration Int; }
/* Additional Integral terms could be added here. For example, the /* Additional Integral terms can be added here. For example, the
following term would account for non-homogeneous Neumann boundary following term may account for non-homogeneous Neumann boundary
conditions, provided that the function nd[] is defined: conditions, provided that the function nd[] is defined:
Integral { [ nd[] , {v} ]; Integral { [ nd[] , {v} ];
In Sur_Neu_Ele; Jacobian Sur; Integration Int; } In Sur_Neu_Ele; Jacobian Sur; Integration Int; }
All the terms are added, and an implicit "= 0" is considered at the end. */ All the terms in the Equation environment are added,
and an implicit "= 0" is considered at the end. */
} }
} }
} }
/* What to do with a weak formulation is specified in the Resolution: here we /* In the Resolution environment we specify what to do with a weak formulation:
simply generale a linear system, solve it and save the solution (.res file) here we simply generate a linear system, solve it and save the solution (.res file)
to disk. */ to disk. */
Resolution { Resolution {
...@@ -291,18 +296,18 @@ Resolution { ...@@ -291,18 +296,18 @@ Resolution {
The first part defines, in terms of the Formulation, which itself refers to The first part defines, in terms of the Formulation, which itself refers to
the FunctionSpace, a number of quantities that can be evaluated at the the FunctionSpace, a number of quantities that can be evaluated at the
postprocessing stage. The three quantities defined here are: postprocessing stage. The three quantities defined here are:
- the scalar electric potential, - the scalar electric potential;
- the electric field, - the electric field;
- the electric displacement. - the electric displacement.
The second part consists in defining groups of post-processing operations, The second part consists in defining groups of post-processing operations,
which can be invoked separately. The first group is invoked by default when which can be invoked separately. The first group is invoked by default when
Gmsh is run interactively. Each Operation specifies Gmsh is run interactively. Each Operation specifies
- a quantity to be eveluated, - a quantity to evaluate;
- the region on which the evaluation is done, - the region on which the evaluation is done;
- the name of the output file. - the name of the output file.
The generated post-processing files are automatically displayed by Gmsh if The generated post-processing files are automatically displayed by Gmsh if
the "Merge result automatically" option is enabled (which is the default). */ the "Merge result automatically" option is enabled (the default). */
PostProcessing { PostProcessing {
{ Name EleSta_v; NameOfFormulation Electrostatics_v; { Name EleSta_v; NameOfFormulation Electrostatics_v;
......
...@@ -12,7 +12,7 @@ ...@@ -12,7 +12,7 @@
------------------------------------------------------------------- */ ------------------------------------------------------------------- */
/* A thing GetDP is pretty good at is the management of global (non-local) basis /* GetDP is pretty good at the management of global (non-local) basis
functions. Finite element expansions typically associate basis functions to functions. Finite element expansions typically associate basis functions to
individual nodes or edges in the mesh. But consider the situation where a individual nodes or edges in the mesh. But consider the situation where a
scalar field is set to be uniform over a region of the problem (say a scalar field is set to be uniform over a region of the problem (say a
...@@ -20,12 +20,12 @@ ...@@ -20,12 +20,12 @@
identical nodal value "v_electrode", a global (non-local) basis function identical nodal value "v_electrode", a global (non-local) basis function
"BF_electrode" is obtained as factor which is the sum of the shape functions "BF_electrode" is obtained as factor which is the sum of the shape functions
of all the nodes in the electrode region. This basis function "BF_electrode" of all the nodes in the electrode region. This basis function "BF_electrode"
- is a continuous function, scalar in this case, - is a continuous function, scalar in this case;
- is equal to 1 at the nodes of the electrode, and to 0 at all other nodes, - is equal to 1 at the nodes of the electrode, and to 0 at all other nodes;
- decreases from 1 to 0 over the one-element-thick layer of elements sharing - decreases from 1 to 0 over the one-element-thick layer of elements sharing
at least one node with the electrode region. at least one node with the electrode region.
One such glabal basis function can be associated with each electrode in the One such global basis function can be associated with each electrode in the
system, so that the finite element expansion of the electric scalar potential system, so that the finite element expansion of the electric scalar potential
reads: reads:
...@@ -34,13 +34,13 @@ ...@@ -34,13 +34,13 @@
with the the sum_k running over all nodes except those of the electrode with the the sum_k running over all nodes except those of the electrode
regions. regions.
We show in this tutorial how GetDP takes advantage of global quantities and This tutorial shows how GetDP takes advantage of global quantities and
the associated global basis functions the associated global basis functions
- to reduce the number of unknowns - to reduce the number of unknowns;
- to compute efficiently the electrode charges "Q_electrode", which are - to efficiently compute the electrode charges "Q_electrode", which are
precisely the energy duals of the global "v_electrode" quantities precisely the energy duals of the global "v_electrode" quantities;
- to deal with floating potentials, which are the computed electrode - to deal with floating potentials, which are the computed electrode
potential when the electrode charge is imposed potentials when the electrode charge is imposed;
- to provide output quantities (charges, armature voltages, capacitances, - to provide output quantities (charges, armature voltages, capacitances,
...) that can be immediately used in a external circuit. */ ...) that can be immediately used in a external circuit. */
...@@ -57,8 +57,8 @@ Group { ...@@ -57,8 +57,8 @@ Group {
/* Abstract regions: /* 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 Sur_Neu_Ele : surface with imposed non homogeneous Neumann boundary conditions
(on n.d = -n . (epsilon grad v)) are imposed (on n.d = -n . (epsilon grad v))
Sur_Electrodes_Ele : electrode regions */ Sur_Electrodes_Ele : electrode regions */
Vol_Ele = Region[ {Air, Diel1} ]; Vol_Ele = Region[ {Air, Diel1} ];
...@@ -66,7 +66,7 @@ Group { ...@@ -66,7 +66,7 @@ Group {
Sur_Electrodes_Ele = Region [ {Ground, Microstrip} ]; Sur_Electrodes_Ele = Region [ {Ground, Microstrip} ];
} }
/* A number of ONELAB parameters are defined to define model parameters or model /* A number of ONELAB parameters are defined to provide model parameters or model
options interactively. */ options interactively. */
MicrostripTypeBC = DefineNumber[0, Name "1Microstrip excitation/Type", MicrostripTypeBC = DefineNumber[0, Name "1Microstrip excitation/Type",
...@@ -86,7 +86,7 @@ Function { ...@@ -86,7 +86,7 @@ Function {
Constraint { Constraint {
/* The Dirichlet boundary condition on the local electric potential is no /* The Dirichlet boundary condition on the local electric potential is no
longer used. The microstrip and the ground are now treated as electrodes, longer used. The microstrip and the ground are herein treated as electrodes,
whose voltage is imposed with the "SetGlobalPotential" constraint below. */ whose voltage is imposed with the "SetGlobalPotential" constraint below. */
{ Name Dirichlet_Ele; Type Assign; { Name Dirichlet_Ele; Type Assign;
Case { Case {
...@@ -95,12 +95,12 @@ Constraint { ...@@ -95,12 +95,12 @@ Constraint {
{ Name SetGlobalPotential; Type Assign; { Name SetGlobalPotential; Type Assign;
Case { Case {
/* Define the imposed potential regionwise on the different parts of /* Impose the potential regionwise on the different parts of
"Sur_Electrodes_Ele". No voltage is imposed to the Microstrip electrode "Sur_Electrodes_Ele". No voltage is imposed to the Microstrip electrode
when the "Fixed charge" option is enabled (if MicrostripTypeBC != 0). */ when the "Fixed charge" option is enabled (if MicrostripTypeBC != 0). */
{ Region Ground; Value 0; } { Region Ground; Value 0; }
If(!MicrostripTypeBC) If(!MicrostripTypeBC)
{ Region Microstrip; Value MicrostripValueBC; } { Region Microstrip; Value MicrostripValueBC; }
EndIf EndIf
} }
} }
...@@ -108,26 +108,26 @@ Constraint { ...@@ -108,26 +108,26 @@ Constraint {
Case { Case {
/* Impose the charge if MicrostripTypeBC != 0 */ /* Impose the charge if MicrostripTypeBC != 0 */
If(MicrostripTypeBC) If(MicrostripTypeBC)
{ Region Microstrip; Value MicrostripValueBC; } { Region Microstrip; Value MicrostripValueBC; }
EndIf EndIf
} }
} }
} }
Group{ Group{
/* The domain of definition lists all regions on which the field "v" is /* The domain of definition comprises all regions on which the field "v" is
defined.*/ defined.*/
Dom_Hgrad_v_Ele = Region[ {Vol_Ele, Sur_Neu_Ele, Sur_Electrodes_Ele} ]; Dom_Hgrad_v_Ele = Region[ {Vol_Ele, Sur_Neu_Ele, Sur_Electrodes_Ele} ];
} }
FunctionSpace { FunctionSpace {
/* The magic in the treatment of global quantitities by GetDP is in the fact /* The magic in the treatment of global quantitities by GetDP lies in the fact
that nearly all the work is done at the level of the FunctionSpace that nearly all the work is done at the level of the FunctionSpace
definition. The finite element expansion of "v" is definition. The finite element expansion of "v" is
v = Sum_k sn_k vn_k + Sum_electrode v_electrode BF_electrode v = Sum_k sn_k vn_k + Sum_electrode v_electrode BF_electrode
with the the sum_k running over all nodes except those of the electrode with the sum_k running over all nodes except those of the electrode
regions. This is exactly what one finds in the FunctionSpace definition regions. This is exactly what one finds in the FunctionSpace definition
below with "sf" standing for "BF_electrode" and "vf" for "v_electrode". below with "sf" standing for "BF_electrode" and "vf" for "v_electrode".
...@@ -135,10 +135,10 @@ FunctionSpace { ...@@ -135,10 +135,10 @@ FunctionSpace {
"GlobalQuantity" section; these names are used in the corresponding "GlobalQuantity" section; these names are used in the corresponding
"GlobalTerm" in the Formulation. Such global terms are the equivalent of a "GlobalTerm" in the Formulation. Such global terms are the equivalent of a
"Integral" term, but where no integration needs to be performed. The "Integral" term, but where no integration needs to be performed. The
"AssociatedWith" statement manifests the fact that the global potential of "AssociatedWith" statement refers to the fact that the global potential of
an electrode is the (electrostatic) energy dual of the electric charge an electrode is the (electrostatic) energy dual of the electric charge
carried by that electrode. Indeed, let us consider the electrostatic weak carried by that electrode. Indeed, let us consider the electrostatic weak
formulation derived in Tutorial 1: find v in Hgradv_Ele such that formulation derived in Tutorial 1: find v in Hgrad_v_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
...@@ -153,14 +153,14 @@ FunctionSpace { ...@@ -153,14 +153,14 @@ FunctionSpace {
charge Q_electrode carried by the electrodes. charge Q_electrode carried by the electrodes.
By checking the "Display global basis functions" checkbox and running the By checking the "Display global basis functions" checkbox and running the
model, you can take a look on how the two "BF_electrode" basis functions in model, you can take a look at how the two "BF_electrode" basis functions in
this model look like. Constraints can then be set on either component of this model look like. Constraints can then be set on either component of
the FunctionSpace. Besides the usual Dirichlet boundary condition on the the FunctionSpace. Besides the usual Dirichlet boundary condition on the
local field, which is left here for the sake of completeness but is not local field (left here for the sake of completeness but not
used in this model, there is the possibility to fix either the used in this model), one may fix either the
GlobalPotential or the ArmatureCharge of each indidual electrode (not both, GlobalPotential or the ArmatureCharge of each indidual electrode (never both,
of course). When the ArmatureCharge is fixed, the computed GlobalPotential of course). When the ArmatureCharge is fixed, the computed GlobalPotential
computed for that electrode is the so-called floating potential. */ for that electrode is the so-called floating potential. */
{ Name Hgrad_v_Ele; Type Form0; { Name Hgrad_v_Ele; Type Form0;
BasisFunction { BasisFunction {
...@@ -207,8 +207,8 @@ Integration { ...@@ -207,8 +207,8 @@ Integration {
} }
Formulation { Formulation {
/* The formulation only contains minor changes compared to formulation from /* The formulation contains only minor changes compared to formulation from
the first tutorial. The global quantities are declared as "Global" in the the first tutorial. The global quantities are declared as "Global" in the
"Quantity" section, and a "GlobalTerm" is added that triggers the assembly "Quantity" section, and a "GlobalTerm" is added that triggers the assembly
of the additional equation per electrode (the "pre-integrated" boundary of the additional equation per electrode (the "pre-integrated" boundary
term) in the system to compute the charge Q_electrode, which term) in the system to compute the charge Q_electrode, which
...@@ -290,8 +290,8 @@ PostProcessing { ...@@ -290,8 +290,8 @@ PostProcessing {
} }
} }
/* Various output results are generated, which are both displayed in the /* Several output results are generated, which are both displayed in the
graphical user interface, and stored in disk files. In particular, global graphical user interface, and stored in disk files. In particular, global
quantities related results are stored in the "output.txt" file. A user option quantities related results are stored in the "output.txt" file. A user option
allows to chose to not overwrite the "output.txt" file when running a new allows to chose to not overwrite the "output.txt" file when running a new
simulation. */ simulation. */
......
...@@ -117,7 +117,7 @@ p6 = newp; Point (p6) = {xCnt+xInt, yCnt-yInt, zCnt+zInt, lc1inf}; ...@@ -117,7 +117,7 @@ p6 = newp; Point (p6) = {xCnt+xInt, yCnt-yInt, zCnt+zInt, lc1inf};
p7 = newp; Point (p7) = {xCnt+xInt, yCnt+yInt, zCnt+zInt, lc1inf}; p7 = newp; Point (p7) = {xCnt+xInt, yCnt+yInt, zCnt+zInt, lc1inf};
p8 = newp; Point (p8) = {xCnt-xInt, yCnt+yInt, zCnt+zInt, lc1inf}; p8 = newp; Point (p8) = {xCnt-xInt, yCnt+yInt, zCnt+zInt, lc1inf};
VolInf[]={}; // Define empty array in case Flag_InfiniteBox is not active
If(Flag_InfiniteBox) If(Flag_InfiniteBox)
xExt = xInt * ratioInf; xExt = xInt * ratioInf;
yExt = yInt * ratioInf; yExt = yInt * ratioInf;
......
...@@ -8,35 +8,39 @@ ...@@ -8,35 +8,39 @@
- Maxwell stress tensor and rigid-body magnetic forces - Maxwell stress tensor and rigid-body magnetic forces
To compute the solution in a terminal: To compute the solution in a terminal:
getdp magnets.pro First generate the (3D) mesh and then run getdp with the chosen resolution
gmsh magnets.geo -3
getdp magnets.pro -solve MagSta_a
OR
getdp magnets.pro -solve MagSta_phi
To compute the solution interactively from the Gmsh GUI: To compute the solution interactively from the Gmsh GUI:
File > Open > magnets.pro File > Open > magnets.pro
Resolution can be chosen from the menu on the left:
MagSta_a (default) or MagSta_phi
Run (button at the bottom of the left panel) Run (button at the bottom of the left panel)
------------------------------------------------------------------- */ ------------------------------------------------------------------- */
/* /*
This rather didactic tutorial solves the electromagnetic field This tutorial solves the electromagnetic field
and the rigid-body forces acting on a set of magnetic pieces and the rigid-body forces acting on a set of magnetic pieces
of parallelepipedic or cylindrical shape. of either parallelepipedic or cylindrical shape.
Besides position and dimension, each piece is attributed Besides position and dimension, each piece is attributed
a (constant) magnetic permability and/or a remanence field. a (constant) magnetic permeability and/or a remanence field.
The pieces are all (a bit abusively) generically called "Magnet" Hereafter, the pieces are all, simply though imprecisely, referred to as "Magnet",
in the problem decription below, irresective of whether they are irresective of whether they are truly permanent magnets or ferromagnetic barrels.
truly permanent magnets or ferromagnetic barrels.
The tutorial model proposes two dual 3D magnetostatic formulations:
The tutorial model proposes both dual 3D magnetostatic formulations: - the magnetic vector potential formulation with spanning-tree gauging;
the magnetic vector potential formulation with spanning-tree gauging, - the scalar magnetic potential formulation.
and the scalar magnetic potential formulation. As there are no conductors, the later is rather simple. The source field "hs" is
The later is rather simple in this case since, as there are no conductors, directly the the known coercive field hc[]:
the known coercive field hc[] is the only source field "hs" one needs
to represent the magnetic field:
h = hs - grad phi , hs = -hc. h = hs - grad phi , hs = -hc.
As in tutorial 2 (magnetostatic field of an electromagnet), a shell If the "Add infinite box" box is ticked, a transformation to infinity shell is
of so-called infinite elements is used here to impose the exact used to impose the exact zero-field boundary condition at infinity.
zero-field boundary condition at infinity. See also Tutorial 2: magnetostatic field of an electromagnet.
The shell is generated automatically by including "InfiniteBox.geo" The shell is generated automatically by including "InfiniteBox.geo"
at the end of the geometrical description of the model. at the end of the geometrical description of the model.
It can be placed rather close of the magnets without loss of accuracy. It can be placed rather close of the magnets without loss of accuracy.
...@@ -44,24 +48,23 @@ ...@@ -44,24 +48,23 @@
The preferred way to compute electromagnetic forces in GetDP The preferred way to compute electromagnetic forces in GetDP
is as an explicit by-product of the Maxwell stress tensor "TM[{b}]", is as an explicit by-product of the Maxwell stress tensor "TM[{b}]",
which is a material dependent function of the magnetic induction "b" field. which is a material dependent function of the magnetic induction "b" field.
Exactly like we computed the heat flux "q(S)" through a surface "S" The magnetic force acting on a rigid body in empty space can be evaluated
using a special auxiliary function "g(S)" associated with that surface as the flux of the Maxwell stress tensor through a surface "S" (surrounding the body).
in "Tutorial 5 : thermal problem with contact resistances", A special auxiliary function "g(S)" linked "S" is defined for each magnet, i.e.
the magnetic force acting on a rigid body in empty space can be evaluated "g(SkinMagnet~{i}) = un~{i}".
as the flux of the Maxwell stress tensor through that surface. The resultant magnetic force acting on "Magnet~{i}" is given by the integral:
There is one auxiliary function "g(SkinMagnet~{i}) = un~{i}"
for each magnet, and the resultant magnetic force acting on "Magnet~{i}"
is given by the integral:
f~{i} = Integral [ TM[{b}] * {-grad un~{i}} ] ; f~{i} = Integral [ TM[{b}] * {-grad un~{i}} ] ;
It should be insisted on the fact that the Maxwell stress tensor This approach is analogous to the computation of heat flux "q(S)" through a
is always discontinuous on material discontinuities, surface "S" described in "Tutorial 5: thermal problem with contact resistances".
Note that the Maxwell stress tensor is always discontinuous on material discontinuities,
and that magnetic forces acting on rigid bodies and that magnetic forces acting on rigid bodies
only depend on the Maxwell stress tensor of empty space, depend only on the Maxwell stress tensor in empty space,
and on the "b" and "h" field distribution, and on the "b" and "h" field distribution,
on the external side of "SkinMagnet~{i}" on the external side of "SkinMagnet~{i}"
(the side of the surface in contact with air). (side of the surface in contact with air).
"{-grad un~{i}}" in the above formula can be regarded "{-grad un~{i}}" in the above formula can be regarded
as the normal vector to "SkinMagnet~{i}" as the normal vector to "SkinMagnet~{i}"
...@@ -79,7 +82,7 @@ ...@@ -79,7 +82,7 @@
Include "magnets_common.pro" Include "magnets_common.pro"
DefineConstant[ DefineConstant[
// preset all getdp options and make them invisible // preset all getdp options and make them (in)visible
R_ = {"MagSta_a", Name "GetDP/1ResolutionChoices", Visible 1, R_ = {"MagSta_a", Name "GetDP/1ResolutionChoices", Visible 1,
Choices {"MagSta_a", "MagSta_phi"}}, Choices {"MagSta_a", "MagSta_phi"}},
C_ = {"-solve -v 5 -v2 -bin", Name "GetDP/9ComputeCommand", Visible 0} C_ = {"-solve -v 5 -v2 -bin", Name "GetDP/9ComputeCommand", Visible 0}
...@@ -198,7 +201,7 @@ Constraint { ...@@ -198,7 +201,7 @@ Constraint {
} }
FunctionSpace { FunctionSpace {
{ Name Hgrad_phi ; Type Form0 ; // scalar magnetic potential { Name Hgrad_phi ; Type Form0 ; // magnetic scalar potential
BasisFunction { BasisFunction {
{ Name sn ; NameOfCoef phin ; Function BF_Node ; { Name sn ; NameOfCoef phin ; Function BF_Node ;
Support Dom_Hgrad_phi ; Entity NodesOf[ All ] ; } Support Dom_Hgrad_phi ; Entity NodesOf[ All ] ; }
...@@ -207,7 +210,7 @@ FunctionSpace { ...@@ -207,7 +210,7 @@ FunctionSpace {
{ NameOfCoef phin ; EntityType NodesOf ; NameOfConstraint phi ; } { NameOfCoef phin ; EntityType NodesOf ; NameOfConstraint phi ; }
} }
} }
{ Name Hcurl_a; Type Form1; // vector magnetic potential { Name Hcurl_a; Type Form1; // magnetic vector potential
BasisFunction { BasisFunction {
{ Name se; NameOfCoef ae; Function BF_Edge; { Name se; NameOfCoef ae; Function BF_Edge;
Support Dom_Hcurl_a ;Entity EdgesOf[ All ]; } Support Dom_Hcurl_a ;Entity EdgesOf[ All ]; }
......
...@@ -14,6 +14,8 @@ ...@@ -14,6 +14,8 @@
To compute the solution interactively from the Gmsh GUI: To compute the solution interactively from the Gmsh GUI:
File > Open > electromagnet.pro File > Open > electromagnet.pro
You may choose the Resolution in the left panel:
Magnetodynamics2D_av (default) or Magnetostatics2D_a
Run (button at the bottom of the left panel) Run (button at the bottom of the left panel)
------------------------------------------------------------------- */ ------------------------------------------------------------------- */
...@@ -32,8 +34,8 @@ Group { ...@@ -32,8 +34,8 @@ Group {
Vol_Mag = Region[{Air, Core, Ind, AirInf}]; // full magnetic domain Vol_Mag = Region[{Air, Core, Ind, AirInf}]; // full magnetic domain
Vol_C_Mag = Region[Core]; // massive conductors Vol_C_Mag = Region[Core]; // massive conductors
Vol_S_Mag = Region[Ind]; // stranded conductors (coils) Vol_S_Mag = Region[Ind]; // stranded conductors (coils)
Vol_Inf_Mag = Region[AirInf]; // annulus for infinite shell transformation Vol_Inf_Mag = Region[AirInf]; // ring-shaped shell for infinite transformation
Val_Rint = rInt; Val_Rext = rExt; // interior and exterior radii of annulus Val_Rint = rInt; Val_Rext = rExt; // interior and exterior radii of ring
} }
Function { Function {
...@@ -50,7 +52,7 @@ Function { ...@@ -50,7 +52,7 @@ Function {
sigma[ Ind ] = 5e7; sigma[ Ind ] = 5e7;
Ns[ Ind ] = 1000 ; // number of turns in coil Ns[ Ind ] = 1000 ; // number of turns in coil
Sc[ Ind ] = SurfaceArea[] ; // surface (cross section) of coil Sc[ Ind ] = SurfaceArea[] ; // area of coil cross section
// Current density in each coil portion for a unit current (will be multiplied // Current density in each coil portion for a unit current (will be multiplied
// by the actual total current in the coil) // by the actual total current in the coil)
js0[ Ind ] = Ns[]/Sc[] * Vector[0,0,-1]; js0[ Ind ] = Ns[]/Sc[] * Vector[0,0,-1];
...@@ -66,7 +68,7 @@ Constraint { ...@@ -66,7 +68,7 @@ Constraint {
} }
{ Name Current_2D; { Name Current_2D;
Case { Case {
// represents the phasor amplitude for a dynamic analysis // represents the phasor amplitude (peak to peak value) for a dynamic analysis
{ Region Ind; Value Current; } { Region Ind; Value Current; }
} }
} }
......
...@@ -29,20 +29,19 @@ ...@@ -29,20 +29,19 @@
the reluctivity. the reluctivity.
Electromagnetic fields expand to infinity. The corresponding boundary Electromagnetic fields expand to infinity. The corresponding boundary
condition can be imposed rigorously by means of a gometrical transformation condition can be imposed rigorously by means of a geometrical transformation
that maps a ring (or shell) of finite elements to the complementary of its that maps a ring (or shell) of finite elements to the complementary of its
interior. As this is a mere geometric transformation, it is enough in the interior. As this is a mere geometric transformation, it is enough in the
model description to attribute a special Jacobian to the ring region model description to attribute a special Jacobian to the ring region
("AirInf") - see the "Jacobian" section below. With this information, GetDP ("AirInf") - see the "Jacobian" section below. With this information, GetDP
is able to deal with the correct transformation of all quantities involved in is able to correctly transform all quantities involved in the model.
the model.
The special Jacobian "VolSphShell" takes 2 parameters in this case, The special Jacobian "VolSphShell" takes 2 parameters in this case,
"Val_Rint" and "Val_Rext", which represent the inner and outer radii of the "Val_Rint" and "Val_Rext", which represent the inner and outer radii of the
transformed ring region and whose value must match those used in the transformed ring region and the value of which must match those used in the
geometrical description of the model (.geo file). This is a typical case geometrical description of the model (.geo file). This is a typical case
where Gmsh and GetDP must consistently share parameter values. To ensure where Gmsh and GetDP must consistently share parameter values. To ensure
consistency in all cases, common parameters are defined is a specific file consistency, common parameters are defined is a specific file
"electromagnet_common.pro", which is included in both the .geo and .pro file "electromagnet_common.pro", which is included in both the .geo and .pro file
of the model. of the model.
...@@ -72,8 +71,8 @@ Group { ...@@ -72,8 +71,8 @@ Group {
The abstract regions in this model have the following interpretation: The abstract regions in this model have the following interpretation:
- Vol_Mag : full volume domain - Vol_Mag : full volume domain
- Vol_S_Mag : region where the current source js is defined - Vol_S_Mag : region with imposed current source js
- Vol_Inf_Mag : region where the infinite ring geometric transformation is applied - Vol_Inf_Mag : region where the infinite ring geometric transformation applies
- Sur_Dir_Mag : part of the boundary with homogenous Dirichlet conditions - Sur_Dir_Mag : part of the boundary with homogenous Dirichlet conditions
- Sur_Neu_Mag : part of the boundary with non-homogeneous Neumann conditions - Sur_Neu_Mag : part of the boundary with non-homogeneous Neumann conditions
*/ */
...@@ -85,7 +84,7 @@ Group { ...@@ -85,7 +84,7 @@ Group {
} }
/* The weak formulation for this problem is derived in a similar way to the /* The weak formulation for this problem is derived in a similar way to the
electrostatic weak formulation from Tutorial 1. The main difference is that electrostatic weak formulation from Tutorial 1. The main differences are that
the fields are now vector-valued, and that we have one linear (source) term the fields are now vector-valued, and that we have one linear (source) term
in addition to the bilinear term. The weak formulation reads: find a such in addition to the bilinear term. The weak formulation reads: find a such
that that
...@@ -109,13 +108,13 @@ Group { ...@@ -109,13 +108,13 @@ Group {
Function { Function {
mu0 = 4.e-7 * Pi; mu0 = 4.e-7 * Pi;
/* New ONELAB variables can then be defined using defineNumber, e.g.: */ /* New ONELAB variables can then be defined using DefineNumber, e.g.: */
murCore = DefineNumber[100, Name "Model parameters/Mur core", murCore = DefineNumber[100, Name "Model parameters/Mur core",
Help "Magnetic relative permeability of Core"]; Help "Magnetic relative permeability of Core"];
/* When the script is run, if the parameter named "Model parameters/Mur core" /* When the script is run, if the parameter named "Model parameters/Mur core"
has not been previously defined, it takes the value (100) provided in has not been previously defined, it takes the value (100) provided in
defineNumber and is sent to the ONELAB server. The "/" character in the DefineNumber and is sent to the ONELAB server. The "/" character in the
variable name is interpreted as a path separator, and results in the variable name is interpreted as a path separator, and results in the
creation of a sub-tree in the graphical user interface. If the script is creation of a sub-tree in the graphical user interface. If the script is
re-run later, the value will be updated using the value from the server re-run later, the value will be updated using the value from the server
...@@ -150,7 +149,7 @@ Constraint { ...@@ -150,7 +149,7 @@ Constraint {
} }
/* In the 2D approximation, the magnetic vector potential a and the current /* In the 2D approximation, the magnetic vector potential a and the current
density js are vectors with a z-component only, i.e.: density js are vectors with only the z-component, i.e.:
a = Vector [ 0, 0, az(x,y) ] a = Vector [ 0, 0, az(x,y) ]
js = Vector [ 0, 0, jsz(x,y) ] js = Vector [ 0, 0, jsz(x,y) ]
...@@ -162,7 +161,7 @@ Constraint { ...@@ -162,7 +161,7 @@ Constraint {
with "BF_PerpendicularEdge" basis functions associated to nodes of the mesh) with "BF_PerpendicularEdge" basis functions associated to nodes of the mesh)
and in the "Hregion_j_Mag_2D" FunctionSpace for js ("Vector" with and in the "Hregion_j_Mag_2D" FunctionSpace for js ("Vector" with
"BF_RegionZ" basis functions associated with the region Vol_S_Mag). Without "BF_RegionZ" basis functions associated with the region Vol_S_Mag). Without
giving all the details, a is thus node-based, whereas js is region-wise providing explicit details, a is thus node-based, whereas js is region-wise
constant. */ constant. */
Group { Group {
...@@ -232,7 +231,7 @@ Integration { ...@@ -232,7 +231,7 @@ Integration {
Integral { [ - Dof{js} , {a} ]; Integral { [ - Dof{js} , {a} ];
In Vol_S_Mag; Jacobian Vol; Integration Int; } In Vol_S_Mag; Jacobian Vol; Integration Int; }
The chosen implementation below is however more effeicient as it avoids The chosen implementation below is however more efficient as it avoids
evaluating repeatedly the function js_fct[] during assembly. evaluating repeatedly the function js_fct[] during assembly.
*/ */
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment