Skip to content

Commit

Permalink
v2.7 (#3)
Browse files Browse the repository at this point in the history
  • Loading branch information
TLCFEM committed Dec 21, 2022
1 parent 0a6750c commit 39580a6
Show file tree
Hide file tree
Showing 28 changed files with 4,230 additions and 350 deletions.
15 changes: 14 additions & 1 deletion docs/Collection/Process/set.md
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,19 @@

The `set` command is used to set some properties of the analysis.

## Linear Elastic Simplification

In a general non-linear context, each substep requires at least two system solving, one for initial guess, one for convergence test.
There is no general way to detect linear system automatically, with which we know the initial guess leads to convergence so that
the further system solving is not necessary. To speed up simulation, one can use the following command to indicate the system is
linear.

```
set linear_system true
```

With this enabled, only one system solving will be performed in each substep.

## Substepping Control

To use a fixed substep size, users can define
Expand Down Expand Up @@ -36,7 +49,7 @@ set min_step_size (1)

By default, an **asymmetric banded** matrix storage scheme is used. For 1D analysis, the global stiffness matrix is
always symmetric. However, when it comes to 2D and 3D analyses, the global stiffness may be structurally symmetric in
terms of its layout, but there may not have the same value on the symmetric entries. In math language, if $$K(i,j)
terms of its layout, but there may not have the same value on the symmetric entries. In maths language, if $$K(i,j)
\neq0$$ then $$K(j,i)\neq0$$, but $$K(i,j)\neq{}K(j,i)$$. Hence, an asymmetric banded storage is the safest.
The `_gbsv()` LAPACK subroutine is used for matrix solving.

Expand Down
Binary file added docs/Example/Solid/wave-explicit.gif
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added docs/Example/Solid/wave-implicit.gif
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
157 changes: 157 additions & 0 deletions docs/Example/Solid/wave-propagation.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,157 @@
# [★★★☆☆] Propagation of A Planar Wave

In this example, we demonstrate an example of wave propagation in a 2D solid.

Essentially, it is a 2D dynamics problem. The [CP4](../../Library/Element/Membrane/Plane/CP4.md) element is used to
model the solid.
Various time integration methods will be used to compare their performance regarding numerical energy dispersion.

The model script can be downloaded [here](wave-propagation.zip).

**This model contains 16641 nodes and 16384 elements. The memory usage is about 1.2 GB.**

**Due to small step size and need to export visualization data, the full analysis takes around 20 minutes to complete on
an average PC platform (6 physical cores @ 3 GHz).**

## Model

A square solid of size $$3200\times3200$$ is used. The structured mesh can be generated by using whatever mesh generator
available. It is not difficult to generate an array of squares using scripting languages such as Python or Matlab.
Here, [Gmsh](https://gmsh.info/) is used.

The left and right boundary are constrained along the horizontal direction. The bottom boundary is constrained along the
vertical direction. An initial velocity is assigned to the node at the centre of the top boundary.

### Material

Whether plane stress or plane strain assumption is adopted is not the focus of this example, we simply use a plane
stress element with a unit thickness.

```text
material Elastic2D 1 1E7 .2 1 1
```

### Visualisation

For visualisation, we define a [`Visualisation`](../../Library/Recorder/Recorder.md) recorder.
We record von Mises stress to represent the propagation of stress field.

```text
hdf5recorder 1 Visualisation MISES width 5
```

### IBC

The boundaries can be extracted by [generating](../../Collection/Define/generate.md) node groups.

```text
generatebyrule nodegroup 1 1 1. 0. # left
generatebyrule nodegroup 2 1 1. 3200. # right
generatebyrule nodegroup 3 2 1. 0. # bottom
```

Then BCs can be applied via [`groupmultiplierbc`](../../Collection/Define/bc.md).

```text
groupmultiplierbc 1 1 1 2
groupmultiplierbc 2 2 3
```

The initial condition can be applied using [`initial`](../../Collection/Define/initial.md) command.

```text
# node 322 is the centre of the top boundary (1600,3200)
initial velocity -1 2 322
```

### Time Integration

We use both implicit and explicit time integration methods.

#### Implicit

The implicit time integration methods are the default.
If no integrator is defined, a default [Newmark](../../Library/Integrator/Newmark/Newmark.md) integrator will be used.

An implicit integrator shall be used with an implicit step.

```text
step ImplicitDynamic 1 4
# or just
# step Dynamic 1 4
```

#### Explicit

Most explicit methods use acceleration as the primary variable, the equations of motion are often expressed as a
function of acceleration.
This differs from implicit methods that often use displacement as the primary variable.
In order to adopt such a difference, one needs to define an `ExplicitDynamic` step, similar to the setting in ABAQUS.

```text
step ExplicitDynamic 1 4
```

Please be aware that most displacement-based constraints **cannot** be used in explicit analysis.

#### Other Settings

Since we are using a linear elastic material with the CP4 elements.
The elemental stiffness is symmetric.
As there are no other non-linear constraints defined in the model, the global stiffness/mass matrix is also symmetric.
It is possible to then turn on symmetric banded storage to save memory space.

Also, since the system is linear, the global stiffness/mass matrix does not change once assembled.
It is possible to indicate the solver to skip iterations.

```text
set symm_mat 1
set band_mat 1
set linear_system 1
```

## Results

For a not-so-rigorous comparison, different spectral radii are used for different methods, mainly for the purpose of
showcasing different methods.

Also, the chosen model parameters are quite arbitrary.
Sufficently accurate results often require an accurate estimation of the highest frequency of the model, which governs the time step size.

### Implicit

The [Bathe](../../Library/Integrator/BatheTwoStep.md) two-step method appears to have the best numerical dispersion
performance among implicit methods.

The [GSSSS](../../Library/Integrator/GSSSS.md) optimal scheme is also fine if the spectral radius is chosen properly.

The second-order, unconditionally stable [Newmark](../../Library/Integrator/Newmark/Newmark.md) method has significant
high-frequency noise.
This explains why it is mainly used for structural dynamics in which the low-frequency response is of interest.

The chosen time step size is 0.01 s for all three cases. The following radii are used:

- Bathe Implicit: 0.9
- GSSSS Optimal: 0.8
- Newmark: 1.0

![implicit methods](wave-implicit.gif)

### Explicit

The explicit methods show better numerical dispersion.

The [Tchamwa](../../Library/Integrator/Tchamwa.md) method is first-order accurate, and does not require corrector.

The [Noh-Bathe](../../Library/Integrator/BatheExplicit.md) two-step explicit method, as discussed in the original
reference, shows superior performance.
However, it requires a corrector step, which requires an additional element-wise computation for each substep.
This increases the computational cost.

The chosen time step size is 0.001 s. The following radii are used:

- Tchamwa: 0.9
- Tchamwa: 0.6
- Bathe Explicit: 0.9

![explicit methods](wave-explicit.gif)
Binary file added docs/Example/Solid/wave-propagation.zip
Binary file not shown.
161 changes: 161 additions & 0 deletions docs/Example/Structural/Dynamics/mass-spring-dashpot-system.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,161 @@
# [★★☆☆☆] Mass-Spring-Dashpot System

This example is taken from the ABAQUS benchmark manual, see section 2.6.1.
The original problem is also reported in the research paper,
see [10.1002/nme.1620170902](https://doi.org/10.1002/nme.1620170902).

The model can be [downloaded](mass-spring-dashpot-system.zip).

## The System

The configuration of the model is shown below.

![mass-spring-dashpot system](mass-spring-dashpot-system.svg)

## Numerical Model

### Nodes

There are two DoFs in the system, two masses are connected to fixed points via nonlinear elastic springs.
In order to do so, we define four nodes.

```text
node 1 0 0
node 2 1 0
node 3 2 0
node 4 3 0
```

Here we use unit distance between two adjacent nodes.
If the elements used are based on strain and strain rate, the unit distance is the only correct choice.
If the elements used are based on displacement and velocity, the unit distance is not a must.

### Materials

The left spring uses a $$\tanh$$ function.
To model it, we use [Tanh1D](../../../Library/Material/Material1D/Elastic/Tanh1D.md) material.

```text
material Tanh1D 1 1000
```

The right spring uses a $$\sinh$$ function.
To model it, we use [Sinh1D](../../../Library/Material/Material1D/Elastic/Sinh1D.md) material.

```text
material Sinh1D 3 1000
```

The middle spring is a linear spring, we simply
use [Elastic1D](../../../Library/Material/Material1D/Elastic/Elastic1D.md) material.

```text
material Elastic1D 2 100
```

The dashpot is linear. We use $$\alpha=1$$
in [Viscosity01](../../../Library/Material/Material1D/Viscosity/Viscosity01.md) material.
The viscosity coefficient is $$5$$.

```text
material Viscosity01 4 1 5
```

### Elements

The springs can be modelled by using either [T2D2](../../../Library/Element/Truss/T2D2.md), which uses strain and strain
rate as the basic quantities, or [Spring01](../../../Library/Element/Special/Spring01.md), which users displacement and
velocity as the basic quantities.

```text
element Spring01 1 1 2 1
element Spring01 2 2 3 2
element Spring01 3 3 4 3
```

For a linear dashpot, we use [Damper01](../../../Library/Element/Special/Damper01.md).

```text
element Damper01 4 2 3 4
```

In addition to the above, it is necessary to define two mass elements.

```text
mass 5 2 1 1
mass 6 3 1 1
```

### IBC

The vertical DoFs of all nodes shall be fixed. The horizontal DoFs of the first and last nodes shall be fixed.

```text
fix2 1 1 1 4
fix2 2 2 1 2 3 4
```

The initial condition can be applied to node 2 via

```text
initial velocity 100 1 2
```

### Load

To apply a step load, one shall use tabular amplitude to define the load curve. The following table can be stored in
file `h`.

```text
0 0
0.499999999999 0
0.5 1
0.55 1
0.550000000001 0
100 0
```

Then the load can be applied such that

```text
amplitude Tabular 1 h
cload 1 1 3000 1 3
```

### Response

To record response, we define two recorders, one for displacement, one for velocity.

```text
hdf5recorder 1 Node U1 2 3
hdf5recorder 2 Node V1 2 3
```

### Step and Analysis

The remaining settings are pretty standard. We use a dynamic step with a default integration
scheme ([Newmark](../../../Library/Integrator/Newmark/Newmark.md)).
If one wishes, other integration schemes can be used.

```text
step dynamic 1 1
set ini_step_size 1E-3
set fixed_step_size 1
converger RelIncreDisp 1 1E-11 10 1
analyze
save recorder 1 2
exit
```

## Results

One could compare the results with the original results in the paper.

![displacement](mass-spring-dashpot-u.svg)

![velocity](mass-spring-dashpot-v.svg)
Loading

0 comments on commit 39580a6

Please sign in to comment.