• Data Preparation
• Uniform Flow in a Rectangular Domain

FINITE ELEMENTS

SYMBOLIC PROGRAMMING IN MAPLE

by Artur Portela

Potential-Flow Example 1: Uniform Flow in a Rectangular Domain

> restart:interface(verboseproc=3):printlevel:=3:

> libname := "C:/mylib/fem",libname:

> with(Plotter):

> with(Cgt_fem):

> with(G_cgt_fem):

>

Data Preparation

The best way to input data for Cgt_fem is to use the procedure read_save_data , which reads a file with the following structure:

* control * [title,point sources,element sources,boundary velocities]
title Potential Flow in a Rectangular Domain
point sources
element sources
boundary velocities

* materials * [material,x-permeability,y-permeability,angle of the local x-direction]
1 2 2 45
2 8 8 45

* nodes * [node,x,y]
1 0 0
2 0 1
3 0 2
4 1 0
5 1 1
6 1 2
7 2 0
8 2 1
9 2 2
10 3 0
11 3 1
12 3 2
13 4 0
14 4 1
15 4 2

* elements * [element,node1,node2,node3,material]
1 1 4 5 1
2 1 5 2 1
3 2 5 3 1
4 3 5 6 1
5 4 7 8 1
6 5 4 8 1
7 5 8 6 1
8 6 8 9 1
9 7 10 11 1
10 7 11 8 1
11 8 11 9 1
12 11 12 9 1
13 10 13 14 1
14 10 14 11 1
15 11 14 12 1
16 12 14 15 1

* constraints * [node,potential]
1 10
2 10
3 10
13 0
14 0
15 0

* boundary velocities * [node1,node2,normal velocity]
1 2 0
2 3 0
5 4 0

* point sources * [node,Q]
1 -50

* element sources * [element,q]
2 30

* end *

The data blocks, with the respective keyword on the top, can be defined in any order.

Alternatively, data can be given manually through the definition of the variables: tcase , control , nods , elems , mat_props , bdr_conds , b_velts , p_srcs and e_srcs . See bellow the structure of these variables.

>

Uniform Flow in a Rectangular Domain

As the first problem, consider the case of a uniform flow in a rectangular domain, 4 m long and 2 m wide. Define flow potentials, uniformly distributed on the two opposite sides, perpendicular to the x direction. Consider the hydraulic conductivity specified as Kx = Ky =2 m/h , in the first half of the domain, and Kx = Ky =8 m/h in the second half, with the angle of the local x -direction, measured from the global x -direction, equal to 45 degrees.

> read_save_data();

read data from a file (y/n) ? y;

file name: "dat_test1.txt";

save data into a file (y/n) ? y;

file name: "check.txt";

The first thing to do is to check the problem data

> plot_problem_data();

> plot_mesh();

After checking the data, run the finite element procedure

> #cgt_fem("anything here will print intermediate results"):

The results obtained are stored in the following variables:

total_kinetic_energy

> n_potls;e_grads;e_velts;n_velts;n_fluxes;n_kined;total_kinetic_energy;

The distribution of the flow potential, along a set of arbitrarilly-specified lines, across the finite element mesh, can be displayed through the following interactive procedure:

> plot_line_potentials();

enter line end-points as X1,Y1,X2,Y2 or 0 to finish: -.5,.5,4.5,.5;

line label: "line A";

color as R,G,B: 1,0,0;

enter line end-points as X1,Y1,X2,Y2 or 0 to finish: .5,2.5,.5,-.5;

line label: "line B";

color as R,G,B: 0,1,0;

enter line end-points as X1,Y1,X2,Y2 or 0 to finish: 0,0,2.,2;

line label: "line C";

color as R,G,B: 0,0,1;

enter line end-points as X1,Y1,X2,Y2 or 0 to finish: 4,0,0,2;

line label: "line D";

color as R,G,B: 1,.8,0;

enter line end-points as X1,Y1,X2,Y2 or 0 to finish: 2,0,2,2;

line label: "line E";

color as R,G,B: 0,1,1;

enter line end-points as X1,Y1,X2,Y2 or 0 to finish: 0;

plot title: "Testing Line Graphics";

3D animated views of the flow potential can be displayed.

> animate_3D_rot(n_potls);

plot title: "Nodal Flow Potential";

> animate_3D_def(n_potls);

plot title: "Nodal Flow Potential";

> plot_contour_lines(n_potls);

plot title: "Nodal Flow Potential";

> plot_velocities();

> animate_velocities();

To display the fluxes through the sides of all the finite elements of the mesh use the following procedure, with no arguments:

> plot_fluxes();

> n_fluxes;

To display the fluxes through the sides of a set of finite elements of the mesh use the same procedure, with the list of the specified elements in argument. Note that, as a consequence, the vector n_fluxes contains the displayed results.

> plot_fluxes([1,2,3,4]);

> n_fluxes;

To restore it, give the commands:

> compute_fluxes():n_fluxes;

> animate_3D_rot(n_kined);

plot title: "Nodal kinetic-Energy Density";

> plot_contour_lines(n_kined);

plot title: "Nodal kinetic-Energy Density";

>